function [fams_ints, fams_ints_raw, usable, fams_counts, gr_vols] = gtCrystComputeExperimentalFamilyIntensities_v2(phase_id, p, det_ind, varargin)
% FUNCTION GTCRYSTCOMPUTEEXPERIMENTALFAMILYINTENSITIES_V2
% [fams_ints, fams_intsp, usable, fams_counts] = 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:
%         only_selected  = <bool>     If true, only selected spots will be
%                                     considered. (The default is false)
%         deselect_annealing_twin_shared_spots
%                        = <bool>     If true, the spots shared by
%                                     annealing twins will be deselected.
%
% 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
        'use_nullspace', true, ...
        '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
        'verbose', false, ...
        '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

%     if conf.deselect_annealing_twin_shared_spots
%         clusters = sfTwinClusters(conf.gr_ids, tot_grains, sample.phases{phase_id}.clusters);
%     else
%         clusters = num2cell(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)));

    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