function varargout = gt6DPlotReconstructionODF(phase_id, grain_id, varargin)

    sample = GtSample.loadFromFile();

    if (numel(grain_id) > 1)
        grain = gtLoadCluster(phase_id, grain_id);
        grain_rec = gtLoadClusterRec(phase_id, grain_id);

        cl_info = sample.phases{phase_id}.getClusterInfo(grain_id);
        if (cl_info.cluster_type == 1)
            grain = grain.samp_ors;
        end
    else
        is_extended = sample.phases{phase_id}.getUseExtended(grain_id);

        grain = gtLoadGrain(phase_id, grain_id, 'fields', {'R_vector'}, 'is_extended', is_extended);
        grain_rec = gtLoadGrainRec(phase_id, grain_id, 'is_extended', is_extended);
    end

    conf = struct( ...
        'percent', 0.2, ...
        'lims_mult', 1, ...
        'resolution_multiplier', 0.1 );
    conf = parse_pv_pairs(conf, varargin);

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

    p = gtLoadParameters();
    om_step = gtAcqGetOmegaStep(p);

    num_orients = numel(grain);

    odf = cell(num_orients, 1);

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

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

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

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

        resolution = tand(om_step) * conf.resolution_multiplier;
        grid_points = [ ...
            min(grain_rec.ODF6D(ii_o).R_vectors, [], 1), ...
            max(grain_rec.ODF6D(ii_o).R_vectors, [], 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{ii_o} = gtGetODFFromGvdm(gv_dm(:, gv_mask), grid_points, gv_int(gv_mask));

        plot_odf(grain(ii_o), odf{ii_o}, grid_points, conf);
    end

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

function plot_odf(grain, odf, grid_points, conf)
    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 (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', conf.lims_mult * get(ax1, 'CLim'));
        set(ax2, 'CLim', conf.lims_mult * get(ax2, 'CLim'));
        set(ax3, 'CLim', conf.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);

    thr_value = (max(odf(:)) - min(odf(:))) * conf.percent * conf.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(1:2) = pos(1:2) * 0.5;
    pos(3:4) = pos(3:4) * 1.5;
    set(f, 'position', round(pos));
end