gt6DPostProcessOrientationSpread.m

function post_result = gt6DPostProcessOrientationSpread(test_data, result, disable_gtdisorientation)
    if (~exist('disable_gtdisorientation', 'var'))
        disable_gtdisorientation = false;
    end

    c = tic();
    fprintf('Cropping volume..')
    bigSize = [size(result.solution{1}, 1), size(result.solution{1}, 2), size(result.solution{1}, 3)];
    finalSize = [size(result.deviation{1}, 1), size(result.deviation{1}, 2), size(result.deviation{1}, 3)];
    lims = round((bigSize - finalSize) /2) + 1;
    lims = [lims, (lims + finalSize - 1)];

    post_result.recon = cellfun(@(x){x(lims(1):lims(4), lims(2):lims(5), lims(3):lims(6))}, result.solution);

    fprintf('\b\b (%f s), computing Average R-Vectors..', toc(c))
    c = tic();
    post_result.avg_orientations = getAverageOrientations(result.grains, post_result.recon, size(test_data.dmvol(:, :, :, 1)));
    fprintf('\b\b (%f s), merging Theoretical Deviations..', toc(c))
    c = tic();
    post_result.domains_theo = mergeBlockTheoDeviations(test_data);
    fprintf('\b\b (%f s), merging Reconstructed Deviations..', toc(c))
    c = tic();
    post_result.domains_recon = mergeBlockReconstructions(test_data, result.grains, post_result.recon);

    fprintf('\b\b (%f s), finding component-wise error Distance (quick)..', toc(c))
    c = tic();
    post_result.distance_comp_deg = 2*atand(getComponentsDistances(test_data, result.grains, post_result.recon));
    fprintf('\b\b (%f s), finding error Distance (quick)..', toc(c))
    c = tic();
    post_result.distance_deg_alt = 2*atand(getDistance(test_data, result.grains, post_result.recon));
    fprintf('\b\b (%f s), finding error Distance (exact)..', toc(c))
    c = tic();
    if (~disable_gtdisorientation)
        post_result.distance_deg = getDisorientation(test_data, result.grains, post_result.recon);
    end
    fprintf('\b\b (%f s), Done.\n', toc(c))

    % And then things like:
%     [ post_result.distance ] = gtCSGetOrientationDistance( recon, testData );
end

function merged = mergeDeviations(deviations, isResult)
    merged = zeros(size(deviations{1}));
    if (isResult)
        values = zeros(size(merged));
    else
        values = ones(size(merged));
    end
    for ii = 1:numel(deviations)
        if (isResult)
            indx = deviations{ii} > values;
        else
            indx = deviations{ii} < values;
        end
        values(indx) = deviations{ii}(indx);
        merged(indx) = ii;
    end
end

function merged = mergeBlockTheoDeviations(testData)
    orientations = testData.gv.dm';
    finalVolSize = size(testData.dmvol(:, :, :, 1));
    tiles = tileOrientations(orientations, 5);

    deviations = getDeviations(orientations, tiles, finalVolSize);

    merged = mergeDeviations(deviations, false);
end

function merged = mergeBlockReconstructions(testData, grains, recon)
    orientations = testData.gv.dm';
    finalVolSize = size(testData.dmvol(:, :, :, 1));
    tiles = tileOrientations(orientations, 5);

    reconRvals = getAverageOrientations(grains, recon, finalVolSize);

    deviations = getDeviations(reconRvals, tiles, finalVolSize);

    merged = mergeDeviations(deviations, false);
end

function distance = getDistance(testData, grains, recon)
    finalVolSize = size(testData.dmvol(:, :, :, 1));
    reconRvals = getAverageOrientations(grains, recon, finalVolSize);

    theoRvals = testData.gv.dm;

    distance = theoRvals - reconRvals';
    distance = sqrt(sum(distance .^ 2, 1));
    distance = reshape(distance, finalVolSize);
end

function distances = getComponentsDistances(testData, grains, recon)
    finalVolSize = size(testData.dmvol(:, :, :, 1));
    reconRvals = getAverageOrientations(grains, recon, finalVolSize);

    theoRvals = testData.gv.dm;

    distances = permute(abs(theoRvals - reconRvals'), [2 3 4 1]);
    distances = reshape(distances, [finalVolSize 3]);
end

function distance = getDisorientation(testData, grains, recon)
    finalVolSize = size(testData.dmvol(:, :, :, 1));
    reconRvals = getAverageOrientations(grains, recon, finalVolSize)';
%     reconRvals = getMaxOrientations(grains, recon, finalVolSize)';

    theoRvals = testData.gv.dm;
    symm = gtCrystGetSymmetryOperators('cubic');

    fprintf('\b\b: ')
    num_vecs = size(theoRvals, 2);
    g_theo = gtMathsRod2OriMat(theoRvals);
    g_recon = gtMathsRod2OriMat(reconRvals);
    for ii = num_vecs:-1:1
        if (mod(ii, 100) == 0)
            num_chars = fprintf('%4d/%04d', ii, num_vecs);
        end
        distance(ii) = gtDisorientation(g_theo(:, :, ii), g_recon(:, :, ii), symm, 'input', 'orimat', 'mode', 'active');
        if (mod(ii, 100) == 0)
            fprintf(repmat('\b', [1 num_chars]));
        end
    end
    distance = reshape(distance, finalVolSize);
end

function tiles = tileOrientations(orientations, numTilesEdge)
    extremes = [min(orientations, [], 1), max(orientations, [], 1)];
    tilesX = linspace(extremes(1), extremes(4), numTilesEdge);
    tilesY = linspace(extremes(2), extremes(5), numTilesEdge);
    tilesZ = linspace(extremes(3), extremes(6), numTilesEdge);
    [tilesX, tilesY, tilesZ] = meshgrid(tilesX, tilesY, tilesZ);
    tiles = [tilesX(:), tilesY(:), tilesZ(:)];
end

function deviation = getDeviations(orientations, tiles, volSize)
    deviation = cell(size(tiles, 1), 1);

    for ii = 1:numel(deviation)
        ddm = orientations - repmat(tiles(ii, :), [size(orientations, 1), 1]);
        ddm = sqrt(sum((ddm .^ 2), 2));
        deviation{ii} = reshape(ddm, volSize);
    end
end

function reconRvals = getMaxOrientations(gr, recon, finalVolSize)
    reconRvals = zeros([3, finalVolSize]);
    mergedVols = mergeDeviations(recon, true);
    mergedVols = reshape(mergedVols, [1, finalVolSize]);
    mergedVols = repmat(mergedVols, [3, 1, 1, 1]);

    for ii = 1:numel(recon)
        indices = mergedVols == ii;
        num_reps = sum(sum(sum(indices))) / 3;
        reconRvals(indices) = repmat(gr{ii}.R_vector', [num_reps, 1]);
    end

    reconRvals = reshape(reconRvals, [3, prod(finalVolSize)])';
end

function [avg_R_vecs, avg_R_vecs_int] = getAverageOrientations(gr, recon_vols, final_vol_size)
    avg_R_vecs = zeros([3, final_vol_size]);
    avg_R_vecs_int = zeros([3, final_vol_size]);

    for ii = 1:numel(recon_vols)
        weights = reshape(recon_vols{ii}, [1, final_vol_size]);
        weights = repmat(weights, [3, 1, 1, 1]);
        avg_R_vecs_int = avg_R_vecs_int + weights;
        avg_R_vecs = avg_R_vecs + weights .* ...
            repmat(gr{ii}.R_vector', [1, final_vol_size]);
    end

    avg_R_vecs = avg_R_vecs ./ (avg_R_vecs_int + (avg_R_vecs_int == 0));
    avg_R_vecs = reshape(avg_R_vecs, [3, prod(final_vol_size)])';
end