  function  vfrac = mmd_dect_im_priority(x, mtype, energies, varargin)
%|function  vfac = mmd_dect_im_priority(x, mtype, energies, [options])
%|
%| image domain multi-material decomposition with a priority list using
%| dual-energy CT images
%| allows only three materials for each pixel
%| 
%| reference paper:
%| mendonca:13:afm TMI
%| A flexible method for multi-material decomposition of dual-energy CT images 
%|
%| in
%|  x	[np, 2]		LAC at two distinct energies
%|  mtype    [L] materials types for multi-material decomposition
%|                  'air' has to be one 
%|  energies(E1, E2)   two distinct energies for image domain decomposition
%|
%| options
%|	'units'		cm | mm	default: cm
%|  'triplets'   library of possible material triplets
%|  'pixmin'    minimum pixel value
%|  'pixmax'   maximum pixel value
%|
%| out 
%|  vfrac  [np L] volume fractions coefficients
%|
%|
%| Yong Long, 2012-12-25
%|


if ~nargin, help(mfilename), error(mfilename), return, end
% todo: add test routine


if isnumeric(mtype) && length(mtype) > 1
	mtype = num2cell(mtype);
end

if ~iscell(mtype)
	mtype = {mtype};
end

L = length(mtype);

% defaults
arg.units = 'cm';
arg.triplets = nchoosek(mtype, 3);
arg.pixmin = 0;
arg.pixmax = 1;
arg = vararg_pair(arg, varargin);


if ~streq(arg.units, 'cm')
	error 'bad units'
end

% image domain multi-material decomposition pixel by pixel
vfrac = mmd_dect_im_pixel(x, mtype, arg.triplets, energies, arg.pixmin, arg.pixmax);

end
  
  
  
%
% mmdect_pixel()
% multi-material decomposition pixel by pixel
%
function vfrac = mmd_dect_im_pixel(x, mtype, triplets, energies, pixmin, pixmax)
% function vfrac = mmd_dect_im_pixel(ig.maskit(fbp), mtype, energies) 
% in
% x	[np, 2]		density pairs
% mtype  L  materials types for multi-material decomposition
%                  'air' has to be one 
% triplets  library of possible material triplets
% energy(E1, E2)   two distinct energies for image domain decomposition
% pixmin, pixmax   constraints    pixmin <= xl <= pixmax 
% out
%  vfrac  [np L] volume fractions coefficients

[np M] = size(x);


if M~= 2
    error 'LAC at two distinct energies'
end

L = length(mtype);

ntri = size(triplets, 1);
mflg = zeros(ntri, L);
vfrac = zeros(np, L);

% lac
for t = 1 : ntri %loop over triplet materials
    ttype = triplets(t, :);        
    lin = xray_read_lac(ttype);
    lac = lin.lac(energies); %[ 2 3]
    % NaN is due to air
    if sum(isnan(lac(:))) >= 1
        lac( find(isnan(lac))) = 0;
    end
    lac_tri{t} = lac;       
    
    for i = 1 : 3 %triplet 
        for j = 1 : L %mtype
            if isequal(mtype{j}, ttype{i})
                mflg(t, j) = 1;
            end
        end
    end
end



%------------------------------------
% image domain decomposition pixel by pixel       
%------------------------------------

for j = 1 : np%50858
    
    printm('%g pixel of %g total', j, np) 
    
    mu = x(j, :); % LAC pair [1 2]
    nopt = 0; % # optimal solutions to volume fractions
    
    for t = 1 : ntri %loop over triplet materials
        
     %   ttype = triplets(t, :);        
        lac = lac_tri{t};
         
        %------------------------------------
        % solve the system equations
        % mu1(E1) a1 + mu2(E1) a2 + mu3(E1) a3 = mu(E1)
        % mu1(E2) a1 + mu2(E2) a2 + mu3(E2) a3 = mu(E2)
        % a1               + a2               + a3              = 1
        % subject to 0 <= a1, a2, a3 <= 1
        %------------------------------------
        A = [lac; ones(1, 3)];
        b = [mu 1]';
        a = A\b;
        
        
        %------------------------------------
        % check whether volume fraction coefficients satisfy 
        % the box constraints
        %------------------------------------   
        if sum(pixmin <= a & a <= pixmax) == 3
          %  a = max(a, 0); let them stay negative?
          %  a = min(a, 1);
          %  nopt = nopt + 1;
           % vfopt(nopt, :) = a; %[nopt, 3]
           % typeopt(nopt) = t; 
            nopt = 1;
            jtype =  t;   
            jvfc = a;
            break
        end
               
    end
    
     
    if nopt == 0%no optimal solution
         %----------------------------------------
         % find feasible solution
         %----------------------------------------
         hf_dis = zeros(ntri, 1); %housdorff distance
         
         for t = 1 : ntri %loop over triplet materials
            ttype = triplets(t, :);        
            lac = lac_tri{t};
            a = wls_simplex(lac, mu', [], [], 'how', 'gp');    
            hf_dis(t) = sum((lac * a - mu') .^ 2);
            vfopt(t, :) = a; %[nopt, 3]          
         end
         
         [dummy nind] = min(hf_dis);
         jtype = nind;
         jvfc = vfopt(nind, :);
         
    end
    
    
        
    vfrac(j, find(mflg(jtype,:) == 1) ) = jvfc;
    
     clear vfopt typeopt
end
  
        

end

























