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