% mmd_true_object.m
% 
% generate true object from multi-material decomposition
% 1. generate 1024 * 1024 object, and then generate sinogram
% and measurements based on it
% 2. downsample it to 512 * 512 to compare with 
% the reconstructed results which are 512 * 512
%
%
% modified NCAT phantom with
% 1. fat
% 2. blood
% 3. iodine contrast agent (6% iodine + 94%water)
% 4. cortical bone
% 5. air
% the order in xftrue
%
%
% Yong Long, 2012-02-14


%save the true fractions
if ~isvar('f.dir'), printm 'data directory'
    f.dir = './true/';
      
    f.home = which('mmd_true_object'); % find setup.m
	[f.home dummy] = fileparts(f.home);
	clear dummy
    
 %   if streq(f.home, '../mmd')
        f.dir = './true/';
 %   elseif streq(f.home, '../mmd')
 %       f.dir = './true/';
 %   else
 %       error('wrong directory')
 %   end
end

% NCAT phantom containing water and bone
% 1024 * 1024
if ~isvar('x'), printm 'ncat phantom'
    f.nx_over = 1024;
    f.ny_over = 1024;
    f.dx_over = 50 / f.nx_over; %cm
   
    ig = image_geom('nx', f.nx_over, 'ny', f.ny_over, 'dx',f.dx_over); %cm
    ig.mask = ig.circ > 0;
    
    x = read_ncat('nx', f.nx_over, 'ny', f.ny_over, 'marrow', true, 'marrow_erode', 6, ...
        'mu', [0 0 1 2 2 1]);
end

% change water to fat
% fat, bone and air volume fractions
if ~isvar('fbone'), printm 'fat, bone and air fractions'
    % 'either water+air=1 or water+bone=1'
    fair = (x <= 1 ) .* (1 - x);
    fbone = (x > 1)  .* (x - 1);
    fwater =1 - fair - fbone;
    
    ffat = fwater; %let water be fat
    clear fwater
    
    cden = [ 0 2];
    cfrac = [ 0 1];
    
    if 0
    im pl 2 2
    im(1, x, 'density', cden), cbar(cden)
    im(2, fair, 'air fraction', cfrac), cbar(cfrac)
    im(3, ffat, 'fat fraction', cfrac), cbar(cfrac)
    im(4, fbone, 'bone fraction', cfrac), cbar(cfrac)
    end
end

% add blood in the heart
if ~isvar('fblood'), printm 'blood fraction'
     ell_blood = [1.5 2.5 3.5 4 0 1
        ];
     fblood = ellipse_im(ig, ell_blood, 'oversample', 3);
    
     ffat = ffat -  fblood;
     ffat = max(0, ffat);
     
     if 0 
     im(ffat, cfrac);
     end
end

% add iodine contrast agent in the lung
if ~isvar('fcont'), printm 'contrast agent fraction'
      ell_cont = [1 1 1 1 0 0.2
         3 4 1 1 0 0.2
        ];
     fcont = ellipse_im(ig, ell_cont, 'oversample', 3);
     fblood = fblood - fcont;
     fblood = max(0, fblood);
          
     
     if 0
     im(fair, cfrac);
     end
     
     % the volume fractions should sum to one
end

% generate fractions
if ~isvar('xftrue'), printm 'xftrue'
      xftrue(:,:,1) = ffat;
      xftrue(:,:,2) = fblood;
      xftrue(:,:,3) = fcont;
      xftrue(:,:,4) = fbone;
      xftrue(:,:,5) = fair;
      
        
      f.name = [f.dir 'true_fractions_1024.fld'];
      fld_write (f.name, xftrue)
end

if  ~isvar('xden'), printm 'true density'
     mtype = {'adipose', 'blood', 'omnipaque300', 'bone', 'air'};
     density = xray_read_dens_yl(mtype);
     LL = length(mtype);
     xden = zeros(ig.nx, ig.ny);
      
     for ll = 1 : LL -1
         xden = xden +  xftrue(:,:,ll) * density(ll);
     end
     
     %save
      fld_write ([f.dir 'true_density_1024.fld'], xden)
      
     if 0
         im(xden, 'true density', cden), cbar(cden)
     end
end

% downsample to 512*512 for compare with reconstructed images
% downsample by averaging
if  ~isvar('xf512'), printm 'downsampled true fractions by averaging'
    % downsample by 4
    for ll = 1 : LL
        xtmp = xftrue(:, :, ll);
        xtmp = ( xtmp(1:2:end, :) + xtmp(2:2:end, :) ) / 2;
        xtmp = ( xtmp(:, 1:2:end) + xtmp(:, 2:2:end) ) / 2;
        xf512(:, :, ll) = xtmp;
    end
    
    % check volume fractions sum to 1
    % minmax( sum(xf512, 3) )
  
    f.name = [f.dir 'true_fractions_down512.fld'];
    fld_write (f.name, xf512)
    
     if 0
         im(xf512, cfrac), cbar(cfrac)
     end

end



if  ~isvar('xden512'), printm 'downsampled density'
     mtype = {'adipose', 'blood', 'omnipaque300', 'bone', 'air'};
     density = xray_read_dens_yl(mtype);
     LL = length(mtype);
     xden512 = zeros(512, 512);
      
     for ll = 1 : LL -1
         xden512 = xden512 +  xf512(:,:,ll) * density(ll);
     end
     
     %save
      fld_write ([f.dir 'true_density_down512.fld'], xden)
      
     if 0
         im(xden512, 'true density downsampled', cden), cbar(cden)
     end
end

%-----------------------------------------------------------
% monoenergy image [HU]
%------------------------------------------------------------
if  ~isvar('xhu'), printm 'monochromatic image'
     
        kev = 70;
        mtype = { 'adipose', 'blood', 'omnipaque300', 'bone', 'air'};
        LL = length(mtype);
        lin = xray_read_lac(mtype(1:LL-1));
        lac = lin.lac(kev); %energy
                

        xhu = zeros(ig.nx, ig.ny);
        
        for ll = 1 : (LL - 1)     
            xhu =xhu + lac(ll) * xftrue(:, :, ll);
        end
        
        %virtual nonenhanced images
        xvir = xhu - lac(3) * xftrue(:, :, 3) +  lac(2) * xftrue(:, :, 3) ;
              
        lin_water = xray_read_lac('water');
        lac_water = lin_water.lac(kev);
       
       
        xhu = xhu * 1000 / lac_water;
        xvir = xvir * 1000 / lac_water;

      %  im_toggle(xhu, xvir), cbar
      
end
