6D/Visualization: added a bunch of figure-making functions for generating...

6D/Visualization: added a bunch of figure-making functions for generating visual output of the reconstruction

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

4_grains/gtCalculateGrain.m

@@ -376,6 +376,6 @@ function sfDisplayFigure(parameters, app, grain, det_ind)
 
     hold(ax, 'Off');
 
-    ax_title = sprintf('Fsim for grain %d in phase %d', grain.id, grain.phaseid);
+    ax_title = sprintf('Fsim for grain %d in phase %d', grain.id(1), grain.phaseid);
     title(ax, ax_title, 'FontWeight', 'bold');
 end

zUtil_Deformation/gt6DCreateProjDataFromGrainCluster.m

@@ -168,7 +168,7 @@ function [refor, estim_space_bbox_pix, estim_orient_bbox] = gt6DCreateProjDataFr
     end
 
     refor = struct( ...
-        'id', refgr.id, 'phaseid', refgr.phaseid, ...
+        'id', refgr.id, 'phaseid', refgr.phaseid, 'gids', grs_list, ...
         'center', estim_space_bbox_mm(1:3) + bbox_size_mm / 2, ...
         'R_vector', estim_orient_bbox(1:3) + bbox_size_rod / 2 );
     refor = gtCalculateGrain(refor, p);

zUtil_Deformation/gt6DPlotReconstructionODF.m

-function gt6DPlotReconstructionODF(phase_id, grain_id, percent)
+function gt6DPlotReconstructionODF(phase_id, grain_id, percent, lims_mult)
     if (numel(grain_id) > 1)
         str_ids = sprintf('_%04d', grain_id);
         file_cluster_rec_path = fullfile('4_grains', ...
             sprintf('phase_%02d', phase_id), ...
             sprintf('grain_cluster_details%s.mat', str_ids));
         grain_rec = load(file_cluster_rec_path);
+
+        resolution_multiplier = 1;
     else
         grain_rec = gtLoadGrainRec(phase_id, grain_id);
+
+        resolution_multiplier = 0.2;
     end
+    grain = gtLoadGrain(phase_id, grain_id(1), 'fields', {'R_vector'});
 
     if (~isfield(grain_rec, 'SEG'))
         grain_rec.SEG = gtLoadGrain(phase_id, grain_id, ...
@@ -28,7 +33,7 @@ function gt6DPlotReconstructionODF(phase_id, grain_id, percent)
     gv_int = int_vol(:);
 
     p = gtLoadParameters();
-    resolution = tand(p.labgeo.omstep/2) * 0.2;
+    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), ...
@@ -42,48 +47,34 @@ function gt6DPlotReconstructionODF(phase_id, grain_id, percent)
     ax2 = subplot(3, 3, 4, 'parent', f);
     ax3 = subplot(3, 3, 7, 'parent', f);
 
-    if (false)
-        imagesc(sum(odf, 3)', 'parent', ax1);
-        imagesc(squeeze(sum(odf, 1))', 'parent', ax2);
-        imagesc(squeeze(sum(odf, 2))', 'parent', ax3);
-    else
-        contourf(sum(odf, 3)', 5, 'parent', ax1)
-        contourf(squeeze(sum(odf, 1))', 5, 'parent', ax2)
-        contourf(squeeze(sum(odf, 2))', 5, 'parent', ax3)
-    end
+    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));
 
-    x_ticks = linspace(1, size(odf, 1), 5);
-    x_tick_labels = 2 * atand(linspace(grid_points(1), grid_points(4), 5) - sum(grid_points([1 4]))/2);
-    x_tick_labels = arrayfun(@(x){sprintf('%5.2f', x)}, x_tick_labels);
+    if (~exist('lims_mult', 'var') || isempty(lims_mult))
+        lims_mult = 1;
+    end
 
-    y_ticks = linspace(1, size(odf, 2), 5);
-    y_tick_labels = 2 * atand(linspace(grid_points(2), grid_points(5), 5) - sum(grid_points([2 5]))/2);
-    y_tick_labels = arrayfun(@(x){sprintf('%5.2f', x)}, y_tick_labels);
+    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);
 
-    z_ticks = linspace(1, size(odf, 3), 5);
-    z_tick_labels = 2 * atand(linspace(grid_points(3), grid_points(6), 5) - sum(grid_points([3 6]))/2);
-    z_tick_labels = arrayfun(@(x){sprintf('%5.2f', x)}, z_tick_labels);
+        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');
-    set(ax1, 'XTick', x_ticks);
-    set(ax1, 'XTickLabel', x_tick_labels);
-    set(ax1, 'YTick', y_ticks);
-    set(ax1, 'YTickLabel', y_tick_labels);
-
     xlabel(ax2, 'Y');
     ylabel(ax2, 'Z');
-    set(ax2, 'XTick', y_ticks);
-    set(ax2, 'XTickLabel', y_tick_labels);
-    set(ax2, 'YTick', z_ticks);
-    set(ax2, 'YTickLabel', z_tick_labels);
-
     xlabel(ax3, 'X');
     ylabel(ax3, 'Z');
-    set(ax3, 'XTick', x_ticks);
-    set(ax3, 'XTickLabel', x_tick_labels);
-    set(ax3, 'YTick', z_ticks);
-    set(ax3, 'YTickLabel', z_tick_labels);
 
     colormap(ax1, gray);
     colormap(ax2, gray);
@@ -93,37 +84,19 @@ function gt6DPlotReconstructionODF(phase_id, grain_id, percent)
     if (~exist('percent', 'var') || isempty(percent))
         percent = 0.2;
     end
-    thr_value = (max(odf(:)) - min(odf(:))) * percent + min(odf(:));
-    p = patch(isosurface(permute(odf, [2 1 3]), thr_value), 'parent', ax4);
+    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')
-
-    x_ticks = linspace(1, size(odf, 1), 15);
-    x_tick_labels = 2 * atand(linspace(grid_points(1), grid_points(4), 15) - sum(grid_points([1 4]))/2);
-    x_tick_labels = arrayfun(@(x){sprintf('%5.2f', x)}, x_tick_labels);
-
-    y_ticks = linspace(1, size(odf, 2), 15);
-    y_tick_labels = 2 * atand(linspace(grid_points(2), grid_points(5), 15) - sum(grid_points([2 5]))/2);
-    y_tick_labels = arrayfun(@(x){sprintf('%5.2f', x)}, y_tick_labels);
-
-    z_ticks = linspace(1, size(odf, 3), 15);
-    z_tick_labels = 2 * atand(linspace(grid_points(3), grid_points(6), 15) - sum(grid_points([3 6]))/2);
-    z_tick_labels = arrayfun(@(x){sprintf('%5.2f', x)}, z_tick_labels);
+    grid(ax4, 'on')
 
     xlabel(ax4, 'X');
     ylabel(ax4, 'Y');
     zlabel(ax4, 'Z');
 
-    set(ax4, 'XTick', x_ticks);
-    set(ax4, 'XTickLabel', x_tick_labels);
-    set(ax4, 'YTick', y_ticks);
-    set(ax4, 'YTickLabel', y_tick_labels);
-    set(ax4, 'ZTick', z_ticks);
-    set(ax4, 'ZTickLabel', z_tick_labels);
-
     set(f, 'Units', 'centimeters')
     pos = get(f, 'Position');
     pos(3:4) = pos(3:4) * 1.5;

zUtil_Deformation/plotting/gt6DMakeClusterIPFFigureForPaper.m

+function gt6DMakeClusterIPFFigureForPaper(cluster, cluster_rec, varargin)
+% function gt6DMakeClusterIPFFigureForPaper(cluster, cluster_rec, varargin)
+
+    conf = struct( ...
+        'plane_normal', [0 0 1], ...
+        'kam_bound', 0.2, ...
+        'filename', '', ...
+        'plane', 'xy', ...
+        'flip_ud', false, ...
+        'flip_lr', false, ...
+        'slice', 1, ...
+        'clims', [0 1], ...
+        'show_labels', false, ...
+        'pixel_to_cm', 0.1, ...
+        'borders', 'no', ...
+        'colorbar', 'no' );
+    conf = parse_pv_pairs(conf, varargin);
+
+    if (~exist('cluster_rec', 'var') || isempty(cluster_rec))
+        cluster_rec = gtLoadClusterRec(cluster.phaseid, cluster.gids);
+    end
+
+    internal_shift = cluster_rec.SEG.segbb(1:3) - cluster_rec.ODF6D.shift + 1;
+    seg_vol_bb = [internal_shift, cluster_rec.SEG.segbb(4:6)];
+
+    seg_vol = logical(cluster_rec.SEG.seg);
+
+    dm_vol(:, :, :, 3) = gtCrop(cluster_rec.ODF6D.voxels_avg_R_vectors(:, :, :, 3), seg_vol_bb) .* seg_vol;
+    dm_vol(:, :, :, 2) = gtCrop(cluster_rec.ODF6D.voxels_avg_R_vectors(:, :, :, 2), seg_vol_bb) .* seg_vol;
+    dm_vol(:, :, :, 1) = gtCrop(cluster_rec.ODF6D.voxels_avg_R_vectors(:, :, :, 1), seg_vol_bb) .* seg_vol;
+
+    slice_rod = get_slice(dm_vol, conf.plane, conf.slice);
+    slice_seg = get_slice(seg_vol, conf.plane, conf.slice);
+    r_vecs = reshape(slice_rod, [], 3);
+    r_vecs = r_vecs(slice_seg(:), :);
+
+    [cmap, ~, ~] = gtIPFCmap(cluster.phaseid, conf.plane_normal, ...
+        'r_vectors', r_vecs, 'background', false);
+
+    slice_R = zeros(size(slice_seg));
+    slice_G = zeros(size(slice_seg));
+    slice_B = zeros(size(slice_seg));
+
+    slice_R(slice_seg(:)) = cmap(:, 1);
+    slice_G(slice_seg(:)) = cmap(:, 2);
+    slice_B(slice_seg(:)) = cmap(:, 3);
+
+    int_vol = gtCrop(cluster_rec.ODF6D.intensity, seg_vol_bb) .* seg_vol;
+    slice_int_vol = get_slice(int_vol, conf.plane, conf.slice, true);
+    slice_dm_vol = get_slice(dm_vol, conf.plane, conf.slice, true);
+
+    % needed for the high anglular boundaries
+    [slice_kam, slice_gam] = gtDefComputeKernelAverageMisorientation(slice_dm_vol, slice_int_vol);
+    slice_kam = squeeze(slice_kam);
+    high_ang = slice_kam > conf.kam_bound;
+
+    slice_R(high_ang) = 0;
+    slice_G(high_ang) = 0;
+    slice_B(high_ang) = 0;
+
+    %
+
+    % Note the inversion between lr/up because of matlab's transpose problem
+    if (conf.flip_ud)
+        slice_R = fliplr(slice_R);
+        slice_G = fliplr(slice_G);
+        slice_B = fliplr(slice_B);
+
+        slice_seg = fliplr(slice_seg);
+    end
+    if (conf.flip_lr)
+        slice_R = flipud(slice_R);
+        slice_G = flipud(slice_G);
+        slice_B = flipud(slice_B);
+
+        slice_seg = flipud(slice_seg);
+    end
+
+    f = figure();
+    ax = axes('parent', f);
+    imagesc(cat(3, slice_R', slice_G', slice_B'), 'parent', ax, ...
+        'AlphaDataMapping', 'scaled', 'AlphaData', slice_seg' );
+    colormap(ax, jet)
+
+    set_size(f, ax, conf, slice_kam);
+
+    drawnow();
+end
+
+function slice = get_slice(vol, plane, slice_idx, no_squeeze)
+    switch(plane)
+        case {'xy', 'XY'}
+            slice = vol(:, :, slice_idx, :);
+        case {'yx', 'YX'}
+            slice = permute(vol(:, :, slice_idx, :), [2 1 3 4]);
+        case {'yz', 'YZ'}
+            slice = vol(slice_idx, :, :, :);
+        case {'zy', 'ZY'}
+            slice = permute(vol(slice_idx, :, :, :), [1 3 2 4]);
+        case {'xz', 'XZ'}
+            slice = vol(:, slice_idx, :, :);
+        case {'zx', 'ZX'}
+            slice = permute(vol(:, slice_idx, :, :), [3 2 1 4]);
+    end
+    if (~exist('no_squeeze', 'var') || isempty(no_squeeze) || ~no_squeeze)
+        slice = squeeze(slice);
+    end
+end
+
+function set_size(f, ax, conf, slice_kam)
+    set(f, 'Units', 'centimeters')
+    img_size = size(slice_kam) * conf.pixel_to_cm;
+    if (strcmpi(conf.borders, 'no'))
+        position = [0, 0, img_size];
+        set(f, 'Position', position)
+        set(f, 'Paperposition', position)
+
+        set(ax, 'Units', 'normalized')
+        set(ax, 'Position', [0, 0, 1, 1])
+    else
+        set(ax, 'Units', 'centimeters')
+        if (strcmpi(conf.colorbar, 'no'))
+            position = [0, 0, img_size] + [0, 0, 2, 1.5];
+            set(f, 'Position', position)
+            set(f, 'Paperposition', position)
+        else
+            position = [0, 0, img_size] + [0, 0, 3, 1.5];
+            set(f, 'Position', position)
+            set(f, 'Paperposition', position)
+
+            cb = colorbar('peer', ax, 'location', 'EastOutside');
+            set(cb, 'Units', 'centimeters')
+            position = [img_size(1)+1.35, 1, 0.65, img_size(2)];
+            set(cb, 'Position', position)
+        end
+        position = [1, 1, img_size];
+        set(ax, 'Position', position)
+    end
+end

zUtil_Deformation/plotting/gt6DMakeClusterKAMFigureForPaper.m

+function gt6DMakeClusterKAMFigureForPaper(cluster, cluster_rec, grs, varargin)
+% function gt6DMakeClusterKAMFigureForPaper(cluster, cluster_rec, grs, varargin)
+
+    conf = struct( ...
+        'percent', 0.0025, ...
+        'filename', '', ...
+        'plane', 'xy', ...
+        'flip_ud', false, ...
+        'flip_lr', false, ...
+        'slice', 1, ...
+        'clims', [0 1], ...
+        'show_labels', false, ...
+        'pixel_to_cm', 0.1, ...
+        'borders', 'yes', ...
+        'colorbar', 'yes' );
+    conf = parse_pv_pairs(conf, varargin);
+
+    if (~exist('cluster_rec', 'var') || isempty(cluster_rec))
+        cluster_rec = gtLoadClusterRec(cluster.phaseid, cluster.gids);
+    end
+    if (~exist('grs', 'var') || isempty(grs))
+        grs = gtLoadGrain(cluster.phaseid, cluster.gids);
+    end
+
+    internal_shift = cluster_rec.SEG.segbb(1:3) - cluster_rec.ODF6D.shift + 1;
+    seg_vol_bb = [internal_shift, cluster_rec.SEG.segbb(4:6)];
+
+    seg_vol = logical(cluster_rec.SEG.seg);
+
+    dm_vol(:, :, :, 3) = gtCrop(cluster_rec.ODF6D.voxels_avg_R_vectors(:, :, :, 3), seg_vol_bb) .* seg_vol;
+    dm_vol(:, :, :, 2) = gtCrop(cluster_rec.ODF6D.voxels_avg_R_vectors(:, :, :, 2), seg_vol_bb) .* seg_vol;
+    dm_vol(:, :, :, 1) = gtCrop(cluster_rec.ODF6D.voxels_avg_R_vectors(:, :, :, 1), seg_vol_bb) .* seg_vol;
+
+    int_vol = gtCrop(cluster_rec.ODF6D.intensity, seg_vol_bb) .* seg_vol;
+
+    [kam, gam] = gtDefComputeKernelAverageMisorientation(dm_vol, int_vol);
+    gam
+    [igm, gos] = gtDefComputeIntraGranularMisorientation(dm_vol, int_vol, 'R_vector', grs(1).R_vector);
+    gos
+
+    slice_kam = get_slice(kam, conf.plane, conf.slice);
+    slice_igm = get_slice(igm, conf.plane, conf.slice);
+    slice_seg = get_slice(seg_vol, conf.plane, conf.slice);
+
+    %
+
+    % Note the inversion between lr/up because of matlab's transpose problem
+    if (conf.flip_ud)
+        slice_kam = fliplr(slice_kam);
+        slice_igm = fliplr(slice_igm);
+        slice_seg = fliplr(slice_seg);
+    end
+    if (conf.flip_lr)
+        slice_kam = flipud(slice_kam);
+        slice_igm = flipud(slice_igm);
+        slice_seg = flipud(slice_seg);
+    end
+
+    clims_kam = [conf.clims(1), min(conf.clims(2), max(slice_kam(:)))];
+
+    f = figure();
+    ax = axes('parent', f);
+    imagesc(slice_kam', 'parent', ax, ...
+        'AlphaDataMapping', 'scaled', 'AlphaData', slice_seg', ...
+        clims_kam );
+    colormap(ax, jet)
+
+    if (conf.show_labels)
+        t = uicontrol('Style', 'Text', 'parent', f, ...
+        'String', sprintf('Grain Average Misorientation: %.3g°', gam), ...
+        'FontWeight', 'bold');
+
+        set_size(f, ax, t, conf, slice_kam);
+    else
+        set_size(f, ax, [], conf, slice_kam);
+    end
+
+    drawnow();
+
+    clims_igm = [conf.clims(1), max(slice_igm(:))];
+
+    f = figure();
+    ax = axes('parent', f);
+    imagesc(slice_igm', 'parent', ax, ...
+        'AlphaDataMapping', 'scaled', 'AlphaData', slice_seg', ...
+        clims_igm );
+    colormap(ax, jet);
+
+    if (conf.show_labels)
+        t = uicontrol('Style', 'Text', 'parent', f, ...
+            'String', sprintf('Grain Orientation Spread: %.3g°', gos), ...
+            'FontWeight', 'bold');
+
+        set_size(f, ax, t, conf, slice_kam);
+    else
+        set_size(f, ax, [], conf, slice_kam);
+    end
+
+    drawnow();
+
+    %
+
+%     if (~isempty(conf.filename) && strcmpi(conf.filename, 'preview'))
+%         [filename, filedir] = uigetfile({'*.png'; '*.eps'});
+%         if (filename)
+%             conf.filename = fullfile(filedir, filename);
+%         else
+%             conf.filename = [];
+%         end
+%     end
+% 
+%     if (~isempty(conf.filename))
+%         [~, ~, ext] = fileparts(conf.filename);
+%         if (isempty(ext))
+%             saveas(f, [conf.filename '.eps']);
+%             saveas(f, [conf.filename '.png']);
+%         else
+%             saveas(f, conf.filename);
+%         end
+%     end
+end
+
+function slice = get_slice(vol, plane, slice_idx)
+    switch(plane)
+        case {'xy', 'XY'}
+            slice = vol(:, :, slice_idx, :);
+        case {'yx', 'YX'}
+            slice = permute(vol(:, :, slice_idx, :), [2 1 3]);
+        case {'yz', 'YZ'}
+            slice = squeeze(vol(slice_idx, :, :, :));
+        case {'zy', 'ZY'}
+            slice = permute(squeeze(vol(slice_idx, :, :, :)), [2 1 3]);
+        case {'xz', 'XZ'}
+            slice = squeeze(vol(:, slice_idx, :, :));
+        case {'zx', 'ZX'}
+            slice = permute(squeeze(vol(:, slice_idx, :, :)), [2 1 3]);
+    end
+end
+
+function set_size(f, ax, t, conf, slice_kam)
+    set(f, 'Units', 'centimeters')
+    img_size = size(slice_kam) * conf.pixel_to_cm;
+    if (strcmpi(conf.borders, 'no'))
+        position = [0, 0, img_size];
+        set(f, 'Position', position)
+        set(f, 'Paperposition', position)
+
+        set(ax, 'Units', 'normalized')
+        set(ax, 'Position', [0, 0, 1, 1])
+    else
+        if (conf.show_labels)
+            upper_padding = 2;
+            lower_offset = 1.75;
+        else
+            upper_padding = 1.5;
+            lower_offset = 1;
+        end
+        set(ax, 'Units', 'centimeters')
+        if (strcmpi(conf.colorbar, 'no'))
+            position = [0, 0, img_size] + [0, 0, 2, upper_padding];
+            set(f, 'Position', position)
+            set(f, 'Paperposition', position)
+        else
+            position = [0, 0, img_size] + [0, 0, 3, upper_padding];
+            set(f, 'Position', position)
+            set(f, 'Paperposition', position)
+
+            cb = colorbar('peer', ax, 'location', 'EastOutside');
+            set(cb, 'Units', 'centimeters')
+            position = [img_size(1)+1.35, lower_offset, 0.65, img_size(2)];
+            set(cb, 'Position', position)
+        end
+        position = [1, lower_offset, img_size];
+        set(ax, 'Position', position)
+
+        if (conf.show_labels)
+            set(t, 'Units', 'centimeters')
+            position = [1, 0.25, img_size(1), 0.75];
+            set(t, 'Position', position)
+        end
+    end
+end

zUtil_Deformation/plotting/gt6DMakeClusterODFFigureForPaper.m

+function gt6DMakeClusterODFFigureForPaper(cluster, cluster_rec, grs, varargin)
+% function gt6DMakeClusterODFFigureForPaper(cluster, cluster_rec, grs, varargin)
+
+    conf = struct( ...
+        'percent', 0.0025, ...
+        'filename', '', ...
+        'clims', [], ...
+        'pixel_to_cm', 0.5, ...
+        'colorbar', 'no' );
+    conf = parse_pv_pairs(conf, varargin);
+
+    if (~exist('cluster_rec', 'var') || isempty(cluster_rec))
+        cluster_rec = gtLoadClusterRec(cluster.phaseid, cluster.gids);
+    end
+    if (~exist('grs', 'var') || isempty(grs))
+        grs = gtLoadGrain(cluster.phaseid, cluster.gids);
+    end
+
+    internal_shift = cluster_rec.SEG.segbb(1:3) - cluster_rec.ODF6D.shift + 1;
+    seg_vol_bb = [internal_shift, cluster_rec.SEG.segbb(4:6)];
+
+    dm_vol(:, :, :, 3) = gtCrop(cluster_rec.ODF6D.voxels_avg_R_vectors(:, :, :, 3), seg_vol_bb);
+    dm_vol(:, :, :, 2) = gtCrop(cluster_rec.ODF6D.voxels_avg_R_vectors(:, :, :, 2), seg_vol_bb);
+    dm_vol(:, :, :, 1) = gtCrop(cluster_rec.ODF6D.voxels_avg_R_vectors(:, :, :, 1), seg_vol_bb);
+
+    int_vol = gtCrop(cluster_rec.ODF6D.intensity, seg_vol_bb);
+
+    gv_dm = gtDefDmvol2Gvdm(dm_vol);
+    gv_mask = logical(cluster_rec.SEG.seg(:));
+    gv_int = int_vol(:);
+
+    p = gtLoadParameters();
+    resolution = tand(p.labgeo.omstep/2);
+    grid_points = [ ...
+        min(cluster_rec.ODF6D.R_vectors{1}, [], 1), ...
+        max(cluster_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));
+
+%     x_pos = 2 * atand(linspace(grid_points(1), grid_points(4), size(odf, 1)) - grs(1).R_vector(1));
+%     y_pos = 2 * atand(linspace(grid_points(2), grid_points(5), size(odf, 2)) - grs(1).R_vector(2));
+%     z_pos = 2 * atand(linspace(grid_points(3), grid_points(6), size(odf, 3)) - grs(1).R_vector(3));
+
+    x_pos = linspace(grid_points(1), grid_points(4), size(odf, 1));
+    y_pos = linspace(grid_points(2), grid_points(5), size(odf, 2));
+    z_pos = linspace(grid_points(3), grid_points(6), size(odf, 3));
+
+    f = figure();
+    ax = axes('parent', f);
+    thr_value = (max(odf(:)) - min(odf(:))) * conf.percent + min(odf(:));
+    p = patch( ...
+        isosurface(x_pos, y_pos, z_pos, permute(odf, [2 1 3]), thr_value), ...
+        'FaceVertexAlphaData', 0.2, 'FaceAlpha', 'flat', 'parent', ax);
+    set(p, 'FaceColor', 'cyan', 'EdgeColor', 'none')
+    camlight('right')
+    camlight('left')
+    lighting('gouraud')
+    grid(ax, 'on')
+    view(ax, 3);
+
+    r_vecs = cat(1, grs(:).R_vector);
+    hold(ax, 'on')
+    scatter3(ax, r_vecs(:, 1), r_vecs(:, 2), r_vecs(:, 3), 60, 'red', 'filled')
+    hold(ax, 'off')
+
+    xlabel(ax, 'X');
+    ylabel(ax, 'Y');
+    zlabel(ax, 'Z');
+
+    for ii = 1:numel(grs)
+        text(r_vecs(ii, 1), r_vecs(ii, 2), r_vecs(ii, 3), ...
+            sprintf('  \\leftarrow %d', grs(ii).id), 'parent', ax)
+    end
+
+%     for ii = 1:numel(grs)
+%         if (grs(ii).id == 171)
+%             text(r_vecs(ii, 1), r_vecs(ii, 2), r_vecs(ii, 3), ...
+%                 sprintf('%d \\rightarrow  ', grs(ii).id), 'parent', ax, ...
+%                 'HorizontalAlignment', 'right')
+%         elseif (grs(ii).id == 3)
+%             text(r_vecs(ii, 1), r_vecs(ii, 2), r_vecs(ii, 3), ...
+%                 sprintf('%d \\rightarrow            ', grs(ii).id), 'parent', ax, ...
+%                 'HorizontalAlignment', 'right')
+%         elseif (grs(ii).id == 67)
+%             text(r_vecs(ii, 1), r_vecs(ii, 2), r_vecs(ii, 3), ...
+%                 sprintf('   \\leftarrow %d', grs(ii).id), 'parent', ax)
+%         elseif (grs(ii).id == 200)
+%             text(r_vecs(ii, 1), r_vecs(ii, 2), r_vecs(ii, 3), ...
+%                 sprintf('   \\leftarrow %d', grs(ii).id), 'parent', ax)
+%         else
+%             text(r_vecs(ii, 1), r_vecs(ii, 2), r_vecs(ii, 3), ...
+%                 sprintf('  \\leftarrow %d', grs(ii).id), 'parent', ax)
+%         end
+%     end
+
+    set_size(f, ax, conf, odf);
+
+    drawnow();
+
+    if (~isempty(conf.filename) && strcmpi(conf.filename, 'preview'))
+        [filename, filedir] = uigetfile({'*.png'; '*.eps'});
+        if (filename)
+            conf.filename = fullfile(filedir, filename);
+        else
+            conf.filename = [];
+        end
+    end
+
+    if (~isempty(conf.filename))
+        [~, ~, ext] = fileparts(conf.filename);
+        if (isempty(ext))
+            saveas(f, [conf.filename '.eps']);
+            saveas(f, [conf.filename '.png']);
+        else
+            saveas(f, conf.filename);
+        end
+    end
+end
+
+function set_size(f, ax, conf, odf)
+    set(f, 'Units', 'centimeters')
+    img_size = mean(size(odf)) * conf.pixel_to_cm;
+    img_size = img_size([1 1]);
+    set(ax, 'Units', 'centimeters')
+    if (strcmpi(conf.colorbar, 'no'))
+        position = [0, 0, img_size] + [0, 0, 4, 1.5];
+        set(f, 'Position', position)
+        set(f, 'Paperposition', position)
+    else
+        position = [0, 0, img_size] + [0, 0, 5, 2];
+        set(f, 'Position', position)
+        set(f, 'Paperposition', position)
+
+        cb = colorbar('peer', ax, 'location', 'EastOutside');
+        set(cb, 'Units', 'centimeters')
+        position = [img_size(1)+2.35, 1, 0.65, img_size(2)];
+        set(cb, 'Position', position)
+    end
+    position = [2, 1, img_size];
+    set(ax, 'Position', position)
+end