 function lin = xray_read_lac(mtype, varargin)
%function lin = xray_read_lac(mtype, [options])
%|
%| Read linear attenuation coefficients for a given material type.
%| Returns a strum that can interpolate onto desired energies.
%| 
%| in
%|	mtype		'aluminum', 'copper', 2, '2', '02-helium', ...
%|			See xray_material_file_name.m
%|			(Optionally can be a cell array of several materials.)
%|
%| options
%|	'units'		cm | mm	default: cm
%|	'interp'	{}	interpolator type.  default {'cubic', 'extrap'}
%|	'shortfile'	0|1	return short file name instead of full path
%|
%| out
%|	lin	strum
%| data:
%|	lin.type	{L}	material type
%|	lin.file	{L}	file for each material
%|	lin.units	'char'	cm | mm
%| methods:
%|	lin.plot		show it
%|	lin.lac(kev)	[N L]	linear attenuation coefficients [cm^-1],
%|				given vector of energies "kev"
%|	max.mean(kev, Ide) [M L] spectrally-weighted "mean" lac
%| private:
%|	lin.lac_raw	{L}	raw linear attenuation coefficients [cm^-1]
%|	lin.lac_kev	{L}	raw energies from tables
%|
%| Copyright 2008-6-15, Jeff Fessler, The University of Michigan
%|
%| modified version of xray_read_mac.m
%| Yong Long, 2011-10-06
%| mas --> lin
%| mac --> lac
%| add omnipaque2: 2% iodine homogenously dissolved in 98% water

if nargin < 1, help(mfilename), error(mfilename), end
if nargin == 1 && streq(mtype, 'test'), xray_read_lac_test, return, end
if nargin == 1 && streq(mtype, 'test-water'), xray_read_lac_water, return, end

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

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

% defaults
lin.units = 'cm';
%arg.interp = {'linear', 'extrap'};
%arg.interp = {'spline', 'extrap'};
lin.interp = {'cubic', 'extrap'};
lin.shortfile = false;

lin = vararg_pair(lin, varargin);
lin.type = mtype;

is_mm = 0;
if streq(lin.units, 'mm')
	is_mm = 1;
elseif ~streq(lin.units, 'cm')
	error 'bad units'
end

if streq(mtype{1}, 'given') % special case {'given', mtype, lac, kev}
	lin.type = mtype{2};
	for ll=1:length(lin.type)
		lin.lac_raw{ll} = mtype{3}(:,ll);
		kev = mtype{4};
		if ncol(kev) > 1
			lin.kev_raw{ll} = kev(:,ll);
		else
			lin.kev_raw{ll} = kev(:,1);
		end
		clear kev
	end

elseif streq(mtype{1}, 'pca_do')
	fail 'not done'
%lac.type = {'pca1', 'pca2'}; % hardwired to first two components
%tmp = type{2}; % arguments for de_component2()
%tmp = de_component2('kev', en, tmp{:});
%lac.lac = tmp.basis; % [ne,2]

else
	for ll=1:length(lin.type)
        mtype = lin.type{ll};
        if streq(mtype, 'omnipaque', 9)
            [lin.lac_raw{ll} lin.kev_raw{ll} lin.file{ll}] = ...
				xray_read_lac_contrast_raw(mtype, is_mm);
        else
            [lin.lac_raw{ll} lin.kev_raw{ll} lin.file{ll}] = ...
                    xray_read_lac_raw(mtype, is_mm);
        end
		if lin.shortfile
			lin.file{ll} = regexprep(lin.file{ll}, '.*\/', '');
		end
	end
end

meth = {'plot', @xray_read_lac_plot, '()';
	'lac', @xray_read_lac_interp, '([kev])';
	'mean', @xray_read_lac_mean, '(kev, Ide)';
	};
lin = strum(lin, meth);

%
% xray_read_lac_contrast_raw()
% lac of constrast agent
% only omnipaque works 
% 2011-11-23
%
function [lac kev file] = xray_read_lac_contrast_raw(mtype, is_mm)

btype = {'iodine', 'water'};
for ll = 1 : length(btype)
    bfile{ll} = xray_material_file_name(btype{ll});
    tmp = load_ascii_skip_header(bfile{ll}); % read uncommented lines
    bkev{ll} = tmp(:,1) * 1000;	% keV
    mac = tmp(:,2);	% mass atten. coeff. [cm^2/g]
    density{ll} = xray_read_dens_yl(btype{ll});
    blac{ll} = mac * density{ll}; % linear atten. coeff. [1/cm]
    if is_mm
        blac{ll} = blac{ll} / 10; 
    end
 end
% interpolate water onto the energy samples of what iodine has
% trick: allow for the k-edge jumps!
interp =  {'cubic', 'extrap'};
tmp = interp1_jump(bkev{2}, log(blac{2}), bkev{1}, ...
		'monodown', interp{:});% interpolate on a log scale
blac{2} = exp(tmp);

% calculate LAC according to concentration(mg/ml)
mgI = str2num(mtype(10:end));
vfrac(1) = mgI / density{1} * 10^-3;
vfrac(2) = 1 - vfrac(1);

lac = zeros(length(blac{1}), 1);
for ll = 1 : 2
    lac = lac + vfrac(ll) * blac{ll};
end

kev = bkev{1};
file = bfile{1};

%
% xray_read_lac_raw()
%
function [lac kev file] = xray_read_lac_raw(mtype, is_mm)

if streq(mtype, 'blood', 5)
    %file = [path_find_dir([filesep 'dect_yl']) filesep 'xray-mass-atten'];
	%file = [file filesep 'compound' filesep mtype];
    file='./xray-mass-atten/blood';
else
    file = xray_material_file_name(mtype);
end
tmp = load_ascii_skip_header(file); % read uncommented lines
kev = tmp(:,1) * 1000;	% keV
mac = tmp(:,2);	% mass atten. coeff. [cm^2/g]
density = xray_read_dens_yl(mtype);
lac = mac * density; % linear atten. coeff. [1/cm]

if is_mm
	%lac = lac * 100;
    lac = lac / 10; 
end

if 0 % look at k-edges
	jump = find(diff(kev) == 0);
	if ~isempty(jump)
		file = regexprep(file, '.*\/', '');
		jump = num2str(kev(jump)', '%5.1f');
		printf('%s %s', file, jump)
	end
end


%
% xray_read_lac_interp()
% interpolate onto desired energies
%
function lac = xray_read_lac_interp(lin, kev)

if ~isvar('kev') | isempty(kev)
	kev = lin.kev{1};
	if length(lin.type) ~= 1
		warn('using energies of 1st material type')
	end
end

LL = length(lin.type);
lac = zeros(length(kev), LL); % [N,L]

% trick: allow for the k-edge jumps!
interp = lin.interp;
for ll = 1:LL
	tmp = log(lin.lac_raw{ll}); % interpolate on a log scale
	lac(:,ll) = interp1_jump(lin.kev_raw{ll}, tmp, kev, ...
		'monodown', interp{:});
end
lac = exp(lac);


%
% xray_read_lac_mean()
% spectrally-weighted "mean" lac
%
function bar = xray_read_lac_mean(lin, kev, Ide)
lac = lin.lac(kev);
bar = diag(1 ./ sum(Ide)) * (Ide' * lac);


%
% xray_read_lac_plot()
%
function xray_read_lac_plot(lin, varargin)

arg.kev = [10:510]';
arg = vararg_pair(arg, varargin);
kev = arg.kev;

lac = lin.lac(kev);

argp = {};
for ll=1:length(lin.type)
	kev_raw = lin.kev_raw{ll};
	ie = min(kev) <= kev_raw & kev_raw <= max(kev);
	argp = {argp{:}, kev_raw(ie), lin.lac_raw{ll}(ie), 'o'};
end

semilogy(argp{:})
mtypes = lin.type;
legend(mtypes{:})
hold on, semilogy(kev, lin.lac(kev), '-'), hold off
xlabel 'KeV', ylabel 'lins atten. coef.', axis tight


%
% xray_read_lac_test()
% example usage
%
function xray_read_lac_test
mtypes = {'lead', 'aluminum', 'water', 'lexan'};
mtypes = {'water', 'bone', 'iodine', 'adipose'};
%mtypes = {'water', 'bone', 'calcium'};
%mtypes = {'04-beryllium', 'water'}
%mtypes = {'iodine', 'cesium', 'csi'};
%mtypes = {'cadmium', 'tellurium', 'cdte'};
%mtypes = {'06', '06-carbon-graphite', 'carbon-graphite'};
lin = xray_read_lac(mtypes);
lin.plot('kev', 20:10:2000)

xrs = xray_read_spectra('ps1');
lin.mean(xrs.en, xrs.Ide);

if 0 % examine CdTe
	% Cd atomic mass about 112
	% Te atomic mass about 128
	en = 5:500;
	tmp = lin.lac(en);
	tmp = [tmp(:,1) * 112 / (112 + 128) + tmp(:,2) * 128 / (112 + 128) ...
		tmp(:,3)];
	plot(en, tmp(:,1) - tmp(:,2))
end


% xray_read_lac_water()
% test lac for water
%
function xray_read_lac_water
wlin = xray_read_lac({'water'});
elin = xray_read_lac({'hydrogen', 'oxygen'});
kev = 40:1:200;
%elin.plot('kev', kev)
tmp = elin.lac(kev);
pred = 2/18 * tmp(:,1) + 16/18 * tmp(:,2); % approximate lac
plot(kev, pred, '-', kev, wlin.lac(kev), '--')
legend('predict', 'nist')
%plot(kev, pred-wlin.lac(kev), '--')
