function [fams_ints, fams_intsp, fams_ints_raw, usable, fams_counts, gr_vols] = gtCrystComputeExperimentalFamilyIntensities_v2(phase_id, p, det_ind, varargin)
% FUNCTION GTCRYSTCOMPUTEEXPERIMENTALFAMILYINTENSITIES_V2
% [fams_ints, fams_intsp, fams_ints_raw, usable, fams_counts, gr_vols] = gtCrystComputeExperimentalFamilyIntensities_v2(phase_id, p, det_ind, varargin)
% INPUTS:
% phase_id = <int> Phase ID. The default is 1.
% p = <struct> The parameter structure.
% det_ind = <int> Detector index. The default is 1.
% OPTIONAL INPUTS:
% fams_ints_norm = <int> Estimate the intensity of each hkl
% family by [1]: median spot intensity,
% or 2: average spot intensity.
% only_selected = <bool> [true]: only use selected spots;
% false: use all the included spots.
% use_spot2grain = <bool> [true]: deselect the shared spots
% indicated by spot2grain.mat;
% mu_att = <double> Attenuation coefficient. If empty,
% absorption volume will be used.
% flat_image = Flat image used for direct beam profile.
% If empty, it will be generated automatically.
% gr_ids = <int> ID of used grains.
% usedfam = <bool> Used hkl families.
% fams_counts_thres = <int> For each grain, only if the spots are
% more than this threshold, the hkl
% family will be used. Default: 4
% abs_vol = Absorption volume.
%
% 06-09-2021 by Zheheng. (based on Nicola's function gtCrystComputeExperimentalFamilyIntensities)
if (~exist('phase_id', 'var') || isempty(phase_id))
phase_id = 1;
end
if (~exist('p', 'var') || isempty(p))
p = gtLoadParameters();
end
if (~exist('det_ind', 'var') || isempty(det_ind))
det_ind = 1;
end
conf = struct(...
'fams_ints_norm', 1, ... % 0 and 1 are median, 2 is average
'only_selected', true, ... % To use included spots or only use selected spots
'use_spot2grain', true, ... % deselect the shared spots indicated by spot2grain.mat
'mu_att', [], ... % 32.2686 cm^-1 43.57keV Ni
'flat_image', [], ...
'gr_ids', [], ...
'usedfam', [], ...
'fams_counts_thres', 4, ...
'abs_vol', []);
conf = parse_pv_pairs(conf, varargin);
conf.mu_att = conf.mu_att * p.acq(det_ind).pixelsize * 0.1; % transform unit from per cm to per pixel: mu(cm^-1)*pixelsize(mm/pixel)*0.1(cm/mm)=mu_att_pixel(pixel^-1)
sample = GtSample.loadFromFile();
tot_grains = sample.phases{phase_id}.getNumberOfGrains();
if isempty(conf.gr_ids)
num_grains = max(tot_grains - 100, round(tot_grains * 0.9));
conf.gr_ids = 1:num_grains;
else
num_grains = numel(conf.gr_ids);
end
fams_ints = nan(num_grains, numel(p.cryst(phase_id).int));
fams_counts = zeros(num_grains, numel(p.cryst(phase_id).int));
if det_ind > 1
conf.use_spot2grain = false;
end
if conf.use_spot2grain
spot2grain = load('4_grains/spot2grain.mat');
spot2grain = spot2grain.spot2grain;
end
gauge = GtGauge(num_grains, 'Loading info: ');
c = tic();
for ii_g = 1:num_grains
gr_id = conf.gr_ids(ii_g);
if sample.phases{phase_id}.selectedGrains(gr_id)
gr = gtLoadGrain(phase_id, gr_id);
gr = gtCalculateGrain(gr, p);
if conf.only_selected
selected = gr.proj(det_ind).selected;
else
selected = true(size(gr.proj(det_ind).included));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate the average blob intensity of each family of each
% grain: fams_ints
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gauge.incrementAndDisplay()
% Working around the fact that we initially implemented a wrong
% formula from Henning's book
if (isfield(gr.proj, 'ondet'))
included = gr.proj(det_ind).ondet(gr.proj(det_ind).included);
else
included = gr.ondet(grs(ii_g).included);
end
bl_ints = arrayfun(@(x) x.intensity, gr.proj(det_ind).bl);
bl_ints = reshape(bl_ints, size(included));
not_nan = ~isnan(bl_ints);
if conf.use_spot2grain
% not using shared spots
to_use = cellfun(@(x) numel(x) == 1, spot2grain(gr.difspotID));
to_use = to_use & selected & not_nan;
else
to_use = selected & not_nan;
end
inc_good = included(to_use);
L_fac = gr.allblobs(det_ind).lorentzfac(inc_good) ./ sind(2 * gr.allblobs(det_ind).theta(inc_good)); % Lorentz factor
P_fac = 1 - (sind(2 * gr.allblobs(det_ind).theta(inc_good)) .* sind(gr.allblobs(det_ind).eta(inc_good))) .^ 2; % polarization factor
theta_types = gr.allblobs(det_ind).thetatype(inc_good);
if isempty(conf.abs_vol)
[atts_tot, ~, conf.abs_vol] = gtGrainComputeBeamAttenuation(gr, p, det_ind, []);
else
atts_tot = gtGrainComputeBeamAttenuation(gr, p, det_ind, conf.abs_vol);
end
if isempty(conf.flat_image)
[beam_ints, conf.flat_image] = gtGrainComputeIncomingBeamIntensity(gr, p, det_ind);
else
beam_ints = gtGrainComputeIncomingBeamIntensity(gr, p, det_ind, conf.flat_image);
end
ints = bl_ints(to_use) ./ L_fac ./ atts_tot(to_use) ./ beam_ints(to_use) ./ P_fac;
fams_temp_counts = accumarray(theta_types, ones(size(theta_types)), [numel(p.cryst(phase_id).int), 1]);
if conf.fams_ints_norm == 2
fams_temp_ints = accumarray(theta_types, ints, [numel(p.cryst(phase_id).int), 1]);
fams_temp_ints = fams_temp_ints ./ (fams_temp_counts + (fams_temp_counts == 0));
else
fams_temp_ints = sfCalcMedianInts(theta_types, ints, [numel(p.cryst(phase_id).int), 1]);
end
fams_ints(ii_g, fams_temp_counts > 0) = fams_temp_ints(fams_temp_counts > 0);
fams_counts(ii_g, :) = fams_temp_counts;
end
end
fprintf('Done in %f seconds.\n', toc(c))
usable = fams_counts > 0;
fams_ints_raw = fams_ints;
if ~isempty(conf.usedfam)
usedfam = false(size(p.cryst.int));
usedfam(conf.usedfam) = true;
fams_ints(:, ~usedfam) = nan;
end
conf.fams_counts_thres = min(conf.fams_counts_thres, max(fams_counts, [], 1));
usedfamsints = bsxfun(@ge, fams_counts, conf.fams_counts_thres);
fams_ints(~usedfamsints) = nan;
[fams_ints, gr_vols] = gtMathsEstimateCoefVecFromIncompDataMatByNullSpace(fams_ints');
% fams_ints contain 1/sin(2theta), which is part of Lorentz factor.
fams_ints = fams_ints' ./ sind(2 * p.cryst(phase_id).theta(1:numel(fams_ints)));
fams_intsp = fams_ints(p.cryst(phase_id).thetatypesp);
fprintf('Theoretical Intensities:\n');
disp(p.cryst(phase_id).int / norm(p.cryst(phase_id).int));
fprintf('Calculated Intensities:\n');
disp(fams_ints);
end
function fams_temp_ints = sfCalcMedianInts(thetatype, ints, fams_size)
fams_temp_ints = zeros(fams_size);
for ii_th = 1:numel(fams_temp_ints)
theta_inds = (thetatype == ii_th);
if any(theta_inds)
fams_temp_ints(ii_th) = median(ints(theta_inds));
end
end
end