% CompMethods() % % This script compares various methods for the forward problem % using an infinite medium since the exact solution for a sphere is % only accurate for an infinite medium. % % Methods include: Exact, 1st Born, Nth Born, Full Born, Ext Born, % and Nth Ext Born % % Feel free to change the parameters to explore different optical % contrast in the object. The code is set for an absorption % perturbation. If you include scattering then don't forget to % switch the ds.Fwd.Method.ObjVec_musp flag to 1. % % Note that at present the Full Born and Ext Born methods do not % work for a scattering perturbation, but Nth Born and Nth Ext Born % do. % % Note that small differences between the exact solution and the % discretized solutions are expected because of small differences % in the volume of the object. % clear ds ds.Debug = 0; %Show Debug Information % % OBJECT POSITION AND RADIUS % xo = 0; yo = 0; zo = 3; ro = 0.4; Thickness = 6; % % SET THE FORWARD MODEL STRUCTURE % ds.Fwd.Lambda = [780]; ds.Fwd.idxRefr = [1.37]; ds.Fwd.Mu_s = [100]; ds.Fwd.g = [0.9]; ds.Fwd.Mu_sp = ds.Fwd.Mu_s .* (1 - ds.Fwd.g); ds.Fwd.v = 2.99e10 ./ ds.Fwd.idxRefr; ds.Fwd.Mu_a = [0.05]; x = [0:0.4:6]; ds.Fwd.Src = SetOptode( 0, 0, 0, 1e-3 ); ds.Fwd.Det = SetOptode( x, 0, Thickness, 1e-3 ); ds.Fwd.ModFreq = 100; ds.Fwd.Boundary.Geometry = 'infinite'; % ONLY INCLUDE THE COMPUTATIONAL VOLUME AROUND THE SPHERE step = ro/5; ds.Fwd.CompVol.Type = 'uniform'; ds.Fwd.CompVol.ZStep = step; ds.Fwd.CompVol.Z=[zo-ro+step/2:step:zo+ro-step/2]; ds.Fwd.CompVol.XStep = step; ds.Fwd.CompVol.X=[xo-ro+step/2:step:xo+ro-step/2]; ds.Fwd.CompVol.YStep = step; ds.Fwd.CompVol.Y=[yo-ro+step/2:step:yo+ro-step/2]; ds.Fwd.MeasList = genMeasList('all', ds.Fwd); ds.Inv = ds.Fwd; % % SET THE OBJECT PARAMETERS % ds.Object = cell(1); ds.Object{1}.Type = 'Sphere'; ds.Object{1}.Pos = [xo yo zo]; ds.Object{1}.Radius = ro; ds.Object{1}.Mu_a = [0.3]; ds.Object{1}.Mu_s = [100]; ds.Object{1}.g = [0.9]; ds.Object{1}.Mu_sp = ds.Object{1}.Mu_s .* (1 - ds.Object{1}.g); % % Perform the Forward Calculations % ds.Fwd.Method.ObjVec_mua = 1; ds.Fwd.Method.ObjVec_musp = 0; ds.Fwd.Method.ObjVec_sd = 0; ds.Fwd.Method.Type = 'BornN'; ds.Fwd.Method.Born_Order = 10; ds.Fwd = genFwdMat(ds.Fwd); ds.Fwd = genMeasData(ds.Fwd, ds.Object, 0); rbn = divide_PhiByPhi( ds.Fwd, ds.Fwd.P.PhiTotal, ... ds.Fwd.P.PhiInc, 1 ); ds.Fwd.Method.Type = 'FullBorn'; ds.Fwd = genMeasData(ds.Fwd, ds.Object, 0); rbf = divide_PhiByPhi( ds.Fwd, ds.Fwd.P.PhiTotal, ... ds.Fwd.P.PhiInc, 1 ); ds.Fwd.Method.Type = 'Born'; ds.Fwd = genFwdMat(ds.Fwd); ds.Fwd = genMeasData(ds.Fwd, ds.Object, 0); rb1 = divide_PhiByPhi( ds.Fwd, ds.Fwd.P.PhiTotal, ... ds.Fwd.P.PhiInc, 1 ); ds.Fwd.Method.Type = 'ExtBorn'; ds.Fwd.Method.ExtBorn_Radius = 0.9; ds.Fwd = genFwdMat(ds.Fwd); ds.Fwd = genMeasData(ds.Fwd, ds.Object, 0); reba = divide_PhiByPhi( ds.Fwd, ds.Fwd.P.PhiTotal, ... ds.Fwd.P.PhiInc, 1 ); ds.Fwd.Method.Type = 'ExtBornN'; ds.Fwd.Method.ExtBorn_Radius = 0.9; ds.Fwd = genFwdMat(ds.Fwd); ds.Fwd = genMeasData(ds.Fwd, ds.Object, 0); rebna = divide_PhiByPhi( ds.Fwd, ds.Fwd.P.PhiTotal, ... ds.Fwd.P.PhiInc, 1 ); ds.Fwd.Method.Type = 'Spherical'; ds.Fwd = genMeasData(ds.Fwd, ds.Object, 0); re = divide_PhiByPhi( ds.Fwd, ds.Fwd.P.PhiTotal, ... ds.Fwd.P.PhiInc, 1 ); figure plot(x,abs(re),x,abs(rb1),x,abs(rbn),x,abs(reba),x,abs(rbf),x,abs(rebna)) legend('Exact','Born1','BornN','ExtBorn','FullBorn','ExtBornN') title('Amplitude Perturbation') xlabel('X Position of Detector (cm)') ylabel('Relative change') if ds.Fwd.ModFreq > 0 figure plot(x,angle(re)*180/pi, ... x,angle(rb1)*180/pi, ... x,angle(rbn)*180/pi, ... x,angle(reba)*180/pi, ... x,angle(rbf)*180/pi, ... x,angle(rebna)*180/pi ) legend('Exact','Born1','BornN','ExtBorn','FullBorn','ExtBornN') title('Phase Perturbation') xlabel('X Position of Detector (cm)') ylabel('Phase change (degrees)') end |