% mmd_setup.m
%
% setup for multi-material decomposition using dect measurements
%
%
% 5 materials:
% 1. fat
% 2. blood
% 3. iodine contrast agent (6% iodine + 94%water)
% 4. cortical bone
% 5. air
%
% Note
% 1. true object is 1024 * 1024 object, and then 
% 2. generate sinogram measurements based on the true object
%
% Yong Long, 2012-02-17

%savefig on


if ~isvar('xftrue'), printm 'load true volume fractions'
    
    f.home = which('mmd_setup'); % find setup.m
	[f.home dummy] = fileparts(f.home);
	clear dummy
    
    if streq(f.home, pwd)
        f.dir = './true/';
    else
        error('wrong directory')
    end
    

    f.name = [f.dir 'true_fractions_1024.fld'];
    xftrue = fld_read(f.name);
          
    xden =  fld_read([f.dir 'true_density_1024.fld']);
    mtype = {'adipose', 'blood', 'omnipaque300', 'bone', 'air'};
        
    if 0
      im clf, im pl 2 3, cfrac = [0 1]; cden= [ 0 2];
      im(1, xftrue(:,:,1), cfrac, 'fat'), cbar
      im(2, xftrue(:,:,2), cfrac, 'blood'), cbar 
      im(3, xftrue(:,:,3), cfrac, 'omnipaque300'), cbar
      im(4, xftrue(:,:,4), cfrac, 'bone'), cbar
      im(5, xftrue(:,:,5), cfrac, 'air'), cbar
      im(6, xden, cden, 'density'), cbar
      
      
      savefig('eps_c', f.dir, 'true_frac_1024')
      end
    
end

    
%
% system matrix, DD
%
if ~isvar('G'), printm 'G'
    [f.nx_over f.ny_over LL] = size(xftrue);
    f.dx_over = 50/f.nx_over; %cm
    
    f.down = 1;
    ig_over = image_geom('nx', f.nx_over, 'ny', f.ny_over, 'dx',f.dx_over, 'down', f.down); %cm
    
    sg = sino_geom('ge1', 'strip_width', 'd', 'units', 'cm', 'down', f.down);
      
    G_over = Gtomo2_dscmex(sg, ig_over);
end


if ~isvar('ssino'), printm 'ssino'
    LL = size(xftrue, 3);
    for ll = 1 : LL-1
        ssino(:,:,ll) = G_over * xftrue(:,:,ll);
    end


end


%
% generate x-ray spectra 
% 
if ~isvar('xrs'), printm  'generate x-ray spectra'  
        
     xrs = xray_read_spectra('poly1,140,80', ...
		'filters', {{'aluminum', 0.25, 'copper', 0.05}});  
    if 0
        figure
        xrs.plot
        title 'ploy1,140,80' 
        savefig('eps_c', f.fig, 'sp_140_80')  
    end
   
end


if ~isvar('ftab'), printm 'ftab - generate x-ray spectra'
   
    sls = de_ftab_sls('max', [50 20 20 20], 'min', [0 0 0 0], 'n', [101 51 51 51]);   
    lin = xray_read_lac(mtype(1:end-1));
    ftab = de_ftab_lac_yl(xrs, lin, 'sls', sls,...
        'ctype', 'newt', ...
		'ctype', 'pre10');   
    
end


if ~isvar('fsino'), printm 'fsino'
	fsino = ftab.fit.fmfun(ssino);

    if 0
	im(sg.s, sg.ad, fsino, 'fsino'), cbar
    savefig('eps_c', f.fig, 'sino')
    end
end


if ~isvar('yi'), printm 'yi'
    
    if 0 %noiseless
        hfi = fsino; % no noise  
        yi = exp(- fsino);
        bi = ones(size(fsino));
    end
    
    if 1%noisy data
        yi = zeros(size(fsino));
        hfi = zeros(size(fsino));
        bi = zeros(size(fsino));
    
        f.I0 = 2 * 10^5;
        f.I0 = f.I0 * xrs.I / xrs.I(1); % xrs.I [1 M]
        f.I0 = f.I0 * f.down;
       
        for mm=1: size(fsino, 3)
            yim =  poisson(f.I0(mm) * exp(-fsino(:,:,mm)), 7) / f.I0(mm);
            
            %Borrow neighbor(s) for any log(0) values
            off = 1; %DECT
        
            while any( yim(:)== 0)
                    ii = find(yim(:) == 0);
                    printm('fixing %d zeros in yi with %d', length(ii), off)
                    iim = max(1, ii - off);
                    iip = min(nd, ii + off);
                    yim(ii) = max(yim(iip), yim(iim));
                    off = off + 1;
            end, clear off ii
            
            yi(:,:,mm) = yim;
            hfi(:,:,mm) = -log(yim);
        end
  
    end

end

%
if ~isvar('xf512'), printm 'downsampled true fractiosn for comparison'
    
     f.name = [f.dir 'true_fractions_down512.fld'];
     xf512 = fld_read(f.name);
end



if ~isvar('ig'), printm 'image geometry for reconstruction'
    f.nx= f.nx_over / 2;
    f.ny= f.ny_over / 2;
    f.dx = 50/f.nx; %cm
    
    ig = image_geom('nx', f.nx, 'ny', f.ny, 'dx',f.dx, 'down', f.down); %cm 
    ig.mask = ig.circ > 0;
    
    G = Gtomo2_dscmex(sg, ig);
end


