% mmd_pl_decom.m
%
% penalized-likelihood (PL) method for multi-material decomposition
% including kappa
%
% run mmd/mmd_setup.m first
%
% 5 materials:
% 1. fat
% 2. blood
% 3. iodine contrast agent (6% iodine + 94%water)
% 4. cortical bone
% 5. air
%
% constraints:
% -a_l <= x_l <= b_l  ( a_l <= b_l)
% x1 + x2 + x3 = 1
%
% Yong Long, 2012-04-04


%savefig off % enable figure saving
savefig on % disable figure saving
     
if ~isvar('yi'), printm 'run setup first'
    mmd_setup
    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
end


%mmd_im_decom.m
if ~isvar('xinit'), printm 'initialize with image-domain method results'
     file_con = [f.dir 'im/im_priority_34215_results.fld'];
     xinit = fld_read(file_con);
     
end

if ~isvar('f.outdir'), printm 'output directory'
    f.outdir = './pl'; % output dir
    f.fig = [f.outdir, 'fig'];
end


if ~isvar('triplets'), printm 'triplets: library of possible material triplets' 
    ntype = {'omnipaque300', 'bone'; 'omnipaque300', 'adipose'};
    triplets = make_material_triplets(mtype, 'ntype', ntype);
    plist = [3 4 2 1 5];
    for i = 1 : size(triplets, 1)
       priority_triplets(i, :) = triplets(plist(i), :);
    end
    
end


%regularization s
if 1 || ~isvar('R'), printm 'R'
    f.beta2l = [8 12 12 8 6] ;  % much better
    f.beta = 2 .^ f.beta2l;
    f.delta = [0.01 0.01 0.005 0.01 0.1];    
  
    W = mean(yi, 3);
    kappa = div0(G' * W, G' * sg.ones);
    
    for ll = 1 : LL
        R{ll} = Reg1(kappa, 'type_denom', 'matlab', 'beta', f.beta(ll), ...
            'pot_arg', {'hyper3', f.delta(ll)}, ...
            'offsets', '2d:hvd'); 
    end
end


if ~isvar('Gb'), printm 'do Gb'
	f.nsubset = 41;
	Gb = Gblock(G, f.nsubset);
end

%
% iterative recon: os
% 
if ~isvar('xh'), printm 'xh recon iterative - may be slow'
    I0 = ones(size(yi));
	data = {yi, I0, ftab};
    dercurv = @de_pl_dercurv;
    
	f.niter = 200;
    f.isave = [0:1:f.niter];
    f.riter = 1;
     
    epsilon = -[3 3 3 3 3]';
    f.epsilon = 10 .^ epsilon;
    f.pixmin = -f.epsilon; % minimum volume fraction value
    f.pixmax = 1 + f.epsilon; % maximum volume fraction value
    
   % PL PWSQS (pixel) + GSMO (symmetric matrix) 
   % materials are coupled
   % June 2013
   
   if 1 
    tic;
   [xh, f.info] = mmd_pl_nsqs_os_fly(ig.maskit(xinit), Gb, ...
                data, R, ...
                'dercurv', dercurv, ...
                'riter', f.riter, ...
                'niter', f.niter, 'isave', f.isave,... 
                'pixmin', f.pixmin, 'pixmax', f.pixmax, ...
                'mtype', mtype, ...
                'triplets', priority_triplets,...
                'chat', 0, ...
                'how', 'mgsmo',...%'cgsmo',...
                'ig', ig,...
                'coupled_material',  1, 'xtrue', xtrue);
            toc
   end
         
end

 f.name = ['id-34215-medfilt-epsilon-' sprintf('%g', epsilon),...
        'size-' sprintf('%g', ig.nx),...
        '-nsub' sprintf('%g', f.nsubset),...
        '-niter' sprintf('%g', f.niter) '-riter' sprintf('%g', f.riter), ...
        '-beta2l-' sprintf('%g', f.beta2l), ...
        '-delta-' sprintf('%3.3f', f.delta)];
    

if ~isvar('cost1'), printm 'PL cost'
	cost1 = de_pl_cost(xh, G, data, ...
        R, ig.mask);
	im clf, plot(f.isave, cost1, '-o')
    xlabel '#iteration', ylabel 'cost'
end

if 1, printm 'rms errors'
    
    rms_err = zeros(length(f.isave), LL);
    indm = ig.mask == 1;
	
    
    for ii = 1 : length(f.isave)
     
        for ll = 1 : LL
            err = col(xh(:, :, ll, ii) -  xf512(:, :, ll));
            rms_err(ii,ll) = rms(err(indm));
        end
        
   
    end
end


if 1, printm 'max err along iterations'

    iter_max_err = zeros(f.niter, LL);
	
   
    for ii = 2 : length(f.isave)
        
        for ll = 1: LL
            iter_max_err(ii-1, ll) = max( col(abs( xh(:, :, ll, ii)- xh(:, :, ll, ii - 1))) );
        end
        
   
        
     end
end


if 1 %save
    save([f.outdir f.name '.mat'], '-v7.3')
end
