gtDefComputeKernelAverageMisorientation.m

function [kam, gam] = gtDefComputeKernelAverageMisorientation(dmvol, intvol, varargin)
    conf = struct('threshold', 5, 'order', 1, 'symm_ops', []);
    conf = parse_pv_pairs(conf, varargin);

    [vol_size(1), vol_size(2), vol_size(3), vol_size(4)] = size(dmvol);

    cum_dists = zeros(vol_size(1:3), 'like', dmvol);
    cum_weights = zeros(vol_size(1:3), 'like', dmvol);

    valid_vol = intvol > 0;

    steps = min(conf.order, vol_size(1:3) - 1);

    for incrx = -steps(1):steps(1)
        for incry = -steps(2):steps(2)
            for incrz = -steps(3):steps(3)
                incr = [incrx, incry, incrz];
                if (all(incr == 0))
                    continue;
                end
                dists = compute_ang_dists_map(dmvol, valid_vol, incrx, incry, incrz, conf.symm_ops);
                dists_weights = cast(dists < conf.threshold, 'like', dmvol);

                shift_pre = incr .* (incr >= 0);
                shift_post = - incr .* (incr < 0);

                pre_weigths = gtPlaceSubVolume(zeros(vol_size(1:3), 'like', dmvol), dists_weights, shift_pre);
                post_weigths = gtPlaceSubVolume(zeros(vol_size(1:3), 'like', dmvol), dists_weights, shift_post);

                pre_dists = gtPlaceSubVolume(zeros(vol_size(1:3), 'like', dmvol), dists, shift_pre);
                post_dists = gtPlaceSubVolume(zeros(vol_size(1:3), 'like', dmvol), dists, shift_post);

                cum_dists = cum_dists + pre_dists .* pre_weigths + post_dists .* post_weigths;
                cum_weights = cum_weights + pre_weigths + post_weigths;
            end
        end
    end

    cum_weights = cum_weights + (cum_weights == 0);
    kam = cum_dists ./ cum_weights;
    gam = sum(kam(:) .* intvol(:)) ./ sum(intvol(:));
end

function dists = compute_ang_dists_map(dmvol, valid_vol, incrx, incry, incrz, symm_ops)
    if (incrx >= 0)
        inds_x_1 = (1+incrx):size(valid_vol, 1);
        inds_x_2 = 1:(size(valid_vol, 1)-incrx);
    else
        inds_x_1 = 1:(size(valid_vol, 1)+incrx);
        inds_x_2 = (1-incrx):size(valid_vol, 1);
    end

    if (incry >= 0)
        inds_y_1 = (1+incry):size(valid_vol, 2);
        inds_y_2 = 1:(size(valid_vol, 2)-incry);
    else
        inds_y_1 = 1:(size(valid_vol, 2)+incry);
        inds_y_2 = (1-incry):size(valid_vol, 2);
    end

    if (incrz >= 0)
        inds_z_1 = (1+incrz):size(valid_vol, 3);
        inds_z_2 = 1:(size(valid_vol, 3)-incrz);
    else
        inds_z_1 = 1:(size(valid_vol, 3)+incrz);
        inds_z_2 = (1-incrz):size(valid_vol, 3);
    end

    valid_diffs = ...
        valid_vol(inds_x_1, inds_y_1, inds_z_1) ...
        .* valid_vol(inds_x_2, inds_y_2, inds_z_2);
    if (~isempty(symm_ops))
        [gvdm_1, dmvol_size] = gtDefDmvol2Gvdm(dmvol(inds_x_1, inds_y_1, inds_z_1, :));
        [gvdm_2, ~] = gtDefDmvol2Gvdm(dmvol(inds_x_2, inds_y_2, inds_z_2, :));

        dists = gtDisorientation(gvdm_1, gvdm_2, symm_ops, 'mode', 'active');
        dists = reshape(dists, dmvol_size(1:3)) .* valid_diffs;
    else
        [gvdm_1, dmvol_size] = gtDefDmvol2Gvdm(dmvol(inds_x_1, inds_y_1, inds_z_1, :));
        [gvdm_2, ~] = gtDefDmvol2Gvdm(dmvol(inds_x_2, inds_y_2, inds_z_2, :));

        diffs = gtMathsRodSum(gvdm_1, -gvdm_2);
        dists = sqrt(sum(diffs .^ 2, 1));
        dists = reshape(dists, dmvol_size(1:3)) .* valid_diffs;
        dists = 2 * atand(dists);
    end
end