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

Note that the factorization method is restricted to plane waves and far field data.

Output Arguments

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

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 $T$ on a Hilbert space over $\bf{C}$ are the self-adjoint operators, see [3],

$$ \mathrm{Re}(T) := \frac{1}{2}(T + T^\ast), \quad \mathrm{Im}(T) := \frac{1}{2\mathrm{i}}(T-T^\ast). $$

These definitions are in accordance with the corresponding definition for complex numbers, see [1]. That does mean that

$$T = \mathrm{Re} T + \mathrm{i}\, \mathrm{Im} T. $$

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 $T$.

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.

$$ z \mapsto \left(\sum_{j = 1}^\infty \frac{\lambda_j}{(\lambda_j + \varepsilon)^2} \left|\langle \Phi_\infty(\cdot,y),\psi_j\rangle_{L^2(\S)}\right|^2\right)^{-1}. $$

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