% ds = ImgRecon( lambda, delta_Mua, delta_Musp, Freq, F_flag, R_flag); % % This function does a simulation of performing the forward % calculation in a slab geometry and then reconstructing an image. % % % lambda - Script handles any number of wavelengths % between 600 and 900nm. Must be a 1xn vector. % % delta_Mua - Change from mu_a = 0.1 [1/cm] for object. % Must be a 1xn vector where n is the number of % wavelengths. % % delta_Musp - Change from mu_sp = 10 [1/cm] for object. % Must be a 1xn vector where n is the number of % wavelengths. % % Freq - Modulation frequency in MHz. This can % specify multiple frequencies in a 1xm vector, % in which case all frequencies will be used % for the reconstruction. % % F_flag, R_flag - vector of length 3, % F_flag specifies the forward problem % R_flag specifies the forward problem % flag(1) - Calculation of Mu_a of Object % (1-Yes, 0-No) % flag(2) - Calculation of Mu_s of Object % (1-Yes, 0-No) % flag(3) - Calculation of SD Pairs (1-Yes, 0-No) % % % function [ds]=ImgRecon( lambda, delta_Mua, delta_Musp, Freq, F_flag, R_flag); % % Local Variables for the script % zObject1 = 2.5; zObject2 = 4.0; Thickness = 6; % Thickness of slab ds.Debug = 0; % Don't Show Debug Information % % SET THE FORWARD MODEL STRUCTURE % nWavelengths = size(lambda,2); ds.Fwd.Lambda = lambda; ds.Fwd.idxRefr = 1.37 * ones(1,nWavelengths); ds.Fwd.Mu_sp = 10 * ones(1,nWavelengths); ds.Fwd.v = 2.99e10 ./ ds.Fwd.idxRefr; ds.Fwd.ModFreq = Freq; %Set Freq ds.Fwd.Mu_a = 0.1 * ones(1,nWavelengths); ds.Fwd.Src = SetOptode( [0:2:6], [0:2:6], 0, 1e-3*ones(1,nWavelengths) ); ds.Fwd.Det = SetOptode( [0:2:6], [0:2:6], Thickness,... 1e-3*ones(1,nWavelengths) ); ds.Fwd.Boundary.Geometry = 'slab'; ds.Fwd.Boundary.Thickness = Thickness; ds.Fwd.CompVol.Type = 'uniform'; ds.Fwd.CompVol.XStep = 0.5; ds.Fwd.CompVol.YStep = 0.5; ds.Fwd.CompVol.ZStep = 0.5; ds.Fwd.CompVol.X = [0 : ds.Fwd.CompVol.XStep : 6]; ds.Fwd.CompVol.Y = [0 : ds.Fwd.CompVol.YStep : 6]; ds.Fwd.CompVol.Z = [1 : ds.Fwd.CompVol.ZStep : Thickness-1]; ds.Fwd.MeasList = genMeasList('all',ds.Fwd); ds.Inv = ds.Fwd; ds.Fwd.Method.Type = 'Born'; % % Specify the object parameters % ds.Object = cell(2); ds.Object{1}.Mu_a = ds.Fwd.Mu_a + delta_Mua; ds.Object{1}.Mu_sp = ds.Fwd.Mu_sp; ds.Object{1}.Type = 'Sphere'; ds.Object{1}.Pos = [1 , 2, zObject1]; ds.Object{1}.Radius = 0.5; ds.Object{2}.Mu_a = ds.Fwd.Mu_a; ds.Object{2}.Mu_sp = ds.Fwd.Mu_sp + delta_Musp; ds.Object{2}.Type = 'Sphere'; ds.Object{2}.Pos = [4, 4, zObject2]; ds.Object{2}.Radius = 0.5; % % Generate Forward Data % ds.Fwd.Method.ObjVec_mua = F_flag(1); ds.Fwd.Method.ObjVec_musp = F_flag(2); ds.Fwd.Method.ObjVec_sd = F_flag(3); ds.Fwd = genFwdMat(ds.Fwd); ds.Fwd = genMeasData(ds.Fwd, ds.Object, ds.Debug); % % Generate Noise % ds.Noise.ShotSNRflag = 0; ds.Noise.ShotSNR = 2000; ds.Noise.ElectronicSNRflag = 0; ds.Noise.ElectronicSNR = 1e-12; ds = genNoise(ds); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Done With the Simulation of the Data Now Time for the Reconstruction % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ds.Inv.Method.Type = 'Rytov'; ds.Inv.Method.ObjVec_mua = R_flag(1); ds.Inv.Method.ObjVec_musp = R_flag(2); ds.Inv.Method.ObjVec_sd = R_flag(3); ds.Inv = genFwdMat(ds.Inv); ds.Recon.ReconAlg = 'Tik'; ds.Recon.TIK_alpha = 1e-2; ds = genRecon(ds); showImage(ds); |