gt6DPlotReconstructionODF.m

function varargout = gt6DPlotReconstructionODF(phase_id, grain_id, percent, lims_mult)
    if (numel(grain_id) > 1)
        grain_rec = gtLoadClusterRec(phase_id, grain_id);

        resolution_multiplier = 1;

        is_extended = false;
    else
        sample = GtSample.loadFromFile();
        is_extended = sample.phases{phase_id}.getUseExtended(grain_id);

        grain_rec = gtLoadGrainRec(phase_id, grain_id, 'is_extended', is_extended);

%         resolution_multiplier = 1;
        resolution_multiplier = 0.2;
    end
    grain = gtLoadGrain(phase_id, grain_id(1), 'fields', {'R_vector'}, 'is_extended', is_extended);

    if (~isfield(grain_rec, 'SEG'))
        grain_rec.SEG = gtLoadGrain(phase_id, grain_id, ...
            'fields', {'seg', 'segbb'});
    end

    internal_shift = grain_rec.SEG.segbb(1:3) - grain_rec.ODF6D.shift + 1;
    seg_vol_bb = [internal_shift, grain_rec.SEG.segbb(4:6)];

    dm_vol(:, :, :, 3) = gtCrop(grain_rec.ODF6D.voxels_avg_R_vectors(:, :, :, 3), seg_vol_bb);
    dm_vol(:, :, :, 2) = gtCrop(grain_rec.ODF6D.voxels_avg_R_vectors(:, :, :, 2), seg_vol_bb);
    dm_vol(:, :, :, 1) = gtCrop(grain_rec.ODF6D.voxels_avg_R_vectors(:, :, :, 1), seg_vol_bb);

    int_vol = gtCrop(grain_rec.ODF6D.intensity, seg_vol_bb);

    gv_dm = gtDefDmvol2Gvdm(dm_vol);
    gv_mask = logical(grain_rec.SEG.seg(:));
    gv_int = int_vol(:);

    p = gtLoadParameters();
    resolution = tand(p.labgeo.omstep/2) * resolution_multiplier;
    grid_points = [ ...
        min(grain_rec.ODF6D.R_vectors{1}, [], 1), ...
        max(grain_rec.ODF6D.R_vectors{1}, [], 1), ...
        resolution([1 1 1])];
    grid_points(1:3) = grid_points(1:3) - mod(grid_points(1:3), resolution);
    grid_points(4:6) = grid_points(4:6) + mod(resolution - mod(grid_points(4:6), resolution), resolution);
    odf = gtGetODFFromGvdm(gv_dm(:, gv_mask), grid_points, gv_int(gv_mask));

    f = figure();
    ax1 = subplot(3, 3, 1, 'parent', f);
    ax2 = subplot(3, 3, 4, 'parent', f);
    ax3 = subplot(3, 3, 7, 'parent', f);

    x_pos = 2 * atand(linspace(grid_points(1), grid_points(4), size(odf, 1)) - grain.R_vector(1));
    y_pos = 2 * atand(linspace(grid_points(2), grid_points(5), size(odf, 2)) - grain.R_vector(2));
    z_pos = 2 * atand(linspace(grid_points(3), grid_points(6), size(odf, 3)) - grain.R_vector(3));

    if (~exist('lims_mult', 'var') || isempty(lims_mult))
        lims_mult = 1;
    end

    if (true)
        imagesc(x_pos, y_pos, sum(odf, 3)', 'parent', ax1);
        imagesc(y_pos, z_pos, squeeze(sum(odf, 1))', 'parent', ax2);
        imagesc(x_pos, z_pos, squeeze(sum(odf, 2))', 'parent', ax3);

        set(ax1, 'CLim', lims_mult * get(ax1, 'CLim'));
        set(ax2, 'CLim', lims_mult * get(ax2, 'CLim'));
        set(ax3, 'CLim', lims_mult * get(ax3, 'CLim'));
    else
        contourf(x_pos, y_pos, sum(odf, 3)', 5, 'parent', ax1)
        contourf(y_pos, z_pos, squeeze(sum(odf, 1))', 5, 'parent', ax2)
        contourf(x_pos, z_pos, squeeze(sum(odf, 2))', 5, 'parent', ax3)
    end

    xlabel(ax1, 'X');
    ylabel(ax1, 'Y');
    xlabel(ax2, 'Y');
    ylabel(ax2, 'Z');
    xlabel(ax3, 'X');
    ylabel(ax3, 'Z');

    colormap(ax1, gray);
    colormap(ax2, gray);
    colormap(ax3, gray);

    ax4 = subplot(3, 3, [2 3 5 6 8 9], 'parent', f);
    if (~exist('percent', 'var') || isempty(percent))
        percent = 0.2;
    end
    thr_value = (max(odf(:)) - min(odf(:))) * percent * lims_mult + min(odf(:));
    p = patch(isosurface(x_pos, y_pos, z_pos, permute(odf, [2 1 3]), thr_value), 'parent', ax4);
    set(p, 'FaceColor', 'cyan', 'EdgeColor', 'none')
    axis(ax4, 'tight')
    camlight('right')
    camlight('left')
    lighting('gouraud')
    grid(ax4, 'on')

    xlabel(ax4, 'X');
    ylabel(ax4, 'Y');
    zlabel(ax4, 'Z');

    set(f, 'Units', 'centimeters')
    pos = get(f, 'Position');
    pos(3:4) = pos(3:4) * 1.5;
    set(f, 'position', pos);

    if (nargout >= 1)
        varargout{1} = odf;
    end
end