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

    fprintf('Post-processing theoretical reconstruction:\n')
    c = tic();
    c_in = c;
    fprintf('- Cropping volume: ')
    rec_size = [size(result.solution{1}, 1), size(result.solution{1}, 2), size(result.solution{1}, 3)];
    result_size = [size(result.expected_dmvol, 1), size(result.expected_dmvol, 2), size(result.expected_dmvol, 3)];
    lims = floor((rec_size - result_size) / 2) + 1;
    lims = [lims, (lims + result_size - 1)];

    post_result.recon = cell(size(result.solution));
    num_vols = numel(post_result.recon);
    sample_rate = 1000;
    num_chars = fprintf('%04d/%04d', 0, num_vols);
    for ii = 1:num_vols
        if (mod(ii, sample_rate) == 0)
            fprintf(repmat('\b', [1, num_chars]));
            fprintf(repmat(' ', [1, num_chars]));
            fprintf(repmat('\b', [1, num_chars]));
            num_chars = fprintf('%04d/%04d', ii, num_vols);
        end

        post_result.recon{ii} = result.solution{ii}(lims(1):lims(4), lims(2):lims(5), lims(3):lims(6));
    end
    fprintf(repmat('\b', [1, num_chars]));
    fprintf(repmat(' ', [1, num_chars]));
    fprintf(repmat('\b', [1, num_chars]));
    post_result.avg_orientations = result.ODF6D.voxels_avg_R_vectors(lims(1):lims(4), lims(2):lims(5), lims(3):lims(6), :);

    tiles = get_orientation_tiles(result.expected_dmvol, 5);

    fprintf('\b\b: Done %d in %g s.\n- Merging Theoretical Deviations..', num_vols, toc(c))
    c = tic();
    post_result.domains_theo = tile_dmvol(result.expected_dmvol, tiles);
    fprintf('\b\b: Done in %g s.\n- Merging Reconstructed Deviations..', toc(c))
    c = tic();
    post_result.domains_recon = tile_dmvol(post_result.avg_orientations, tiles);

    fprintf('\b\b: Done in %g s.\n- Finding component-wise error Distance (quick)..', toc(c))
    c = tic();
    post_result.distance_comp_deg = 2 * atand(abs(result.expected_dmvol - post_result.avg_orientations));
    fprintf('\b\b: Done in %g s.\n- Finding error Distance (quick)..', toc(c))
    c = tic();
    post_result.distance_deg_alt = 2 * atand(sqrt(sum((result.expected_dmvol - post_result.avg_orientations) .^ 2, 4)));
    fprintf('\b\b: Done in %g s.\n- Finding error Distance (exact)..', toc(c))
    c = tic();
    if (~disable_gtdisorientation)
        post_result.distance_deg = get_disorientation(test_data, result, post_result.avg_orientations);
    end
    fprintf('\b\b: Done in %g s.\nDone in %g.\n', toc(c), toc(c_in))
end

function tiles = get_orientation_tiles(dmvol, num_dim_tiles)
    gvdm = gtDefDmvol2Gvdm(dmvol);
    extremes = [min(gvdm, [], 2); max(gvdm, [], 2)];

    tilesX = linspace(extremes(1), extremes(4), num_dim_tiles);
    tilesY = linspace(extremes(2), extremes(5), num_dim_tiles);
    tilesZ = linspace(extremes(3), extremes(6), num_dim_tiles);

    [tilesX, tilesY, tilesZ] = ndgrid(tilesX, tilesY, tilesZ);

    tiles = [tilesX(:), tilesY(:), tilesZ(:)];
end

function tiling = tile_dmvol(dmvol, tiles)
    distances = cell(size(tiles, 1), 1);
    for ii = 1:numel(distances)
        tile_r_vec = reshape(tiles(ii, :), 1, 1, 1, 3);
        dm_diff_vol = bsxfun(@minus, dmvol, tile_r_vec);
        distances{ii} = sqrt(sum((dm_diff_vol .^ 2), 4));
    end

    tiling = zeros(size(distances{1}));
    curr_tile_dist = Inf(size(tiling));
    for ii = 1:numel(distances)
        indx = distances{ii} < curr_tile_dist;
        curr_tile_dist(indx) = distances{ii}(indx);
        tiling(indx) = ii;
    end
end

function distances = get_disorientation(test_data, result, voxels_avg_R_vectors)
    vol_size = size(result.expected_dmvol(:, :, :, 1));

    theo_gvdm = gtDefDmvol2Gvdm(result.expected_dmvol);
    reco_gvdm = gtDefDmvol2Gvdm(voxels_avg_R_vectors);
    cryst = test_data.parameters.cryst;
    symm = gtCrystGetSymmetryOperators(cryst.crystal_system, cryst.spacegroup);

    fprintf('\b\b: ')
    num_vecs = size(theo_gvdm, 2);
    chunk_size = 1000;
    distances = zeros(num_vecs, 1);
    for ii = 1:chunk_size:num_vecs
        num_chars = fprintf('%04d/%04d', ii, num_vecs);

        last_chunk_ind = min(num_vecs, ii + chunk_size - 1);

        distances(ii:last_chunk_ind) = gtDisorientation( ...
            theo_gvdm(:, ii:last_chunk_ind), ...
            reco_gvdm(:, ii:last_chunk_ind), ...
            symm, 'input', 'orimat', 'mode', 'passive');

        fprintf(repmat('\b', [1 num_chars]));
        fprintf(repmat(' ', [1 num_chars]));
        fprintf(repmat('\b', [1 num_chars]));
    end
    distances = reshape(distances, vol_size);
end