gt6DGetMatchingReflections.m

function [matches_ref_to_twin, matches_twin_to_ref] = gt6DGetMatchingReflections(ref_or, ors, angular_tols, det_ind)
    % matches_twin_to_ref contains the indices of shared reflections in the
    % twin_gr.allblobs list: i.e. [0 10; ...] would tell us that
    % ref.allblobs(1) shares a spot with twin_or(2).allblobs(10)
    %
    % matches_ref_to_twin contains the indices of shared reflections in the
    % ref_gr.allblobs list: i.e. [10 0; ...] would tell us that
    % twin_or(1).allblobs(1) shares a spot with ref_gr.allblobs(10)

    if (~exist('angular_tols', 'var') || isempty(angular_tols))
        angular_tols = 1; % degrees
    end
    angular_tols = reshape(angular_tols, 1, []);
    if (~exist('det_ind', 'var') || isempty(det_ind))
        det_ind = 1;
    end

    num_ors = numel(ors);

    if (numel(angular_tols) == 1 && num_ors > 1)
        angular_tols = angular_tols(1, ones(num_ors, 1));
    end

    ref_ws = ref_or.allblobs(det_ind).omega;
    ref_ns = ref_or.allblobs(det_ind).eta;
    ref_tt = ref_or.allblobs(det_ind).thetatype;
    ref_Lfactors = ref_or.allblobs(det_ind).lorentzfac;

    num_spots = numel(ref_ws);
    angular_tols = angular_tols(ones(num_spots, 1), :);

    w_tol = angular_tols .* ref_Lfactors(:, ones(num_ors, 1));
    n_tol = angular_tols;

    ors_allblobs = [ors(:).allblobs];
    ors_ws = cat(2, ors_allblobs(:).omega);
    ors_ns = cat(2, ors_allblobs(:).eta);
    ors_tt = cat(2, ors_allblobs(:).thetatype);

    ones_array = ones(num_spots, num_ors);
    matches_ref_to_twin = zeros(num_spots, num_ors);
    matches_twin_to_ref = zeros(num_spots, num_ors);

    for ii_refl = 1:num_spots
        ref_ws_array = ref_ws(ii_refl) * ones_array;
        ref_ns_array = ref_ns(ii_refl) * ones_array;
        ref_tt_array = ref_tt(ii_refl) * ones_array;
        tt_equal = ref_tt_array == ors_tt;
        w_diff = abs(mod(ref_ws_array - ors_ws + 180, 360) - 180);
        n_diff = abs(mod(ref_ns_array - ors_ns + 180, 360) - 180);
        [twin_refl_inds, twin_inds] = find((w_diff < w_tol) & (n_diff < n_tol) & tt_equal);

        for jj = 1:numel(twin_refl_inds)
            matches_twin_to_ref(twin_refl_inds(jj), twin_inds(jj)) = ii_refl;
            matches_ref_to_twin(ii_refl, twin_inds(jj)) = twin_refl_inds(jj);
        end
    end
end