function post_result = gt6DPostProcessOrientationSpread(test_data, result, disable_gtdisorientation, remove_axes)
    if (~exist('disable_gtdisorientation', 'var') || isempty(disable_gtdisorientation))
        disable_gtdisorientation = false;
    end
    if (~exist('remove_axes', 'var'))
        remove_axes = [];
    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 = gt6DComputeOrientationTiles(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 = gt6DTileDmvol(result.expected_dmvol, tiles);
    fprintf('\b\b: Done in %g s.\n- Merging Reconstructed Deviations..', toc(c))
    c = tic();
    post_result.domains_recon = gt6DTileDmvol(post_result.avg_orientations, tiles);

    fprintf('\b\b: Done in %g s.\n- Finding component-wise error Distance (quick)..', toc(c))
    c = tic();
    [gvdm_1, dmvol_size] = gtDefDmvol2Gvdm(result.expected_dmvol);
    [gvdm_2, ~] = gtDefDmvol2Gvdm(post_result.avg_orientations);
    if ~isempty(remove_axes)
        gvdm_2(remove_axes, :) = gvdm_1(remove_axes, :);
    end
    diffs = gtMathsRodSum(gvdm_1, -gvdm_2);
    diffs = gtDefGvdm2Dmvol(diffs, dmvol_size);
    post_result.distance_comp_deg = 2 * atand(abs(diffs));
    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(diffs .^ 2, 4)));
    post_result.distance_deg_alt = post_result.distance_deg_alt .* single(result.expected_seg);
    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 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