
% 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); 



