factMeth
Computes an image of the shape of the scatterer from far field pattern using the factorization method.
Contents
Syntax
z = factMeth(seti);
Description
z = factMeth(seti) computes a real matrix z.
Input Arguments
- seti : structure array
- seti.incNb : number of transmitters (number of incident fields)
- seti.delta : Relative (artificial) noise level of data (default: 0.01)
- seti.FmeasDelta : scattered field at receivers positions (i.e. 'the data') (uScaRX = uTotRX - uIncRX, i.e. total field minus incident field) (complex matrix of size seti.measNb x seti.incNb)
Note that the factorization method is restricted to plane waves and far field data.
Output Arguments
- z : Result of the numercial implementation of the factorization method. (real matrix of size seti.nROI x seti.nROI)
More About
Note that the factorization method's result separates the obstacle from the background but does not distinguish whether the contrast is different from zero in the real or imaginary part. This is important for the interpretation of the result.
References
- [1] Armin Lechleiter. Factorization Methods for Photonics and Rough Surfaces. PhD thesis, Universität Karlsruhe, 2008. URL: https://doi.org/10.5445/KSP/1000008735.
- [2] Willy Dörfler, Armin Lechleiter, Michael Plum, Guido Schneider, and Christian Wieners. Photonic Crystals: Mathematical Analysis and Numerical Approximation, volume 42 of Oberwolfach Seminars. Birkhäuser, Basel, 2011. URL: https://doi.org/10.1007/978-3-0348-0113-3.
- [3] Andreas Kirsch and Natalia Grinberg. The Factorization Method for Inverse Problems, volume 36 of Oxford Lecture Series in Mathematics and its Applications. Oxford University Press, 2008. URL: https://doi.org/10.1093/acprof:oso/9780199213535.001.0001.
See Also
Code: function: factMeth
function z = factMeth(seti)
Load data F (farField)
Notation: F: data (far field)
farField = seti.FmeasDelta; % requires far field data (as factorization method is restricted to far field)
Code: function: factMeth: Construct F_# and Compute the SVD
Construct F_# = Re(F) + Im(F) and compute its SVD.
Construct F_#: 1st step: Re(F) and Im(F)
Note that the real and imaginary parts of a bounded linear operator on a Hilbert space over
are the self-adjoint operators, see [3],
These definitions are in accordance with the corresponding definition for complex numbers, see [1]. That does mean that
This does not mean that in the numerical implementation the real part only contains real numbers and the imaginary part only imaginary ones for a matrix .
realPart = (farField + farField')/2; % Re(F) imagPart = (farField - farField')/(2*1i); % Im(F)
Construct F_#: 2nd step: Re(F)
[V,D] = eig(realPart);
absRealPart = V*abs(D)*V'; % |Re(F)|
Construct: F_#: 3rd step: F_#
farFieldSharp = absRealPart + imagPart; % F_# = |Re(F)| + Im(F)
Compute SVD of F_#
[V,D,~] = svd(farFieldSharp); % V will be used for \psi
Code: function: factMeth: Numerical Implementation of Factorization Method
The factorization method decides whether a point belongs to the image or not. In theory, we have (4.32) in Th. 4.2.6 in [2]. As we have to deal with a finite number of receivers and perturbed data we follow (4.33) with a regularization due to noise, i.e.
Preparations
k = seti.k; % wave number s = zeros(size(seti.G(seti.qROIexact))); % ROI as matrix directions = seti.incPnts; % The directions are the coordinates of transmitters as well as receivers. % Actually it does not matter because they have the same positions. % Anyway, the same number of transmitters and receivers is needed (incNb = measNb). % Use grid points from ROI x1ROI = seti.gridROI(2,1:seti.nROI); gridPnts = x1ROI;
Numercial Implementation of Series of Factorization Method
reg = seti.facRegFac*seti.delta; % regularization parameter Neig = floor(seti.measNb*seti.facNeigFac); % numbers of eigenvalues used in series before truncated if Neig > seti.measNb error('Neig is greater than seti.measNb. Set seti.facNeigFac not bigger than 1.'); end for l = 1:length(gridPnts) for m = 1:length(gridPnts) testfunc = exp(-1i*k* (directions(1,:)*gridPnts(l)+directions(2,:)*gridPnts(m))); % \Phi_\infty % requires: incNb = measNb (e.g. 64) % testfunc: (e.g. size(testfunc) is 1 x 64) % V: eigenvectors in cols (e.g. size(V) is 64 x 64) % Use less eigenvalues and eigenvectors: % Reduce number of eigenvalues and eigenvectors. D = D(1:Neig,1:Neig); V = V(:,1:Neig); % Series of factorization method: s(l,m) = 1./sum(diag(D)./(diag(D)+reg).^2.*transpose(abs(testfunc*V)).^2); % (4.33) in [2] end end z = s;
Rotate the result Rotate the result such that it fits (probably because of the grid)
z = transpose(z); % z = rot90(rot90((z))); % rotate 180 degrees
end