From 99432467dfd650eeaf6b78bd659b7d3057bfad78 Mon Sep 17 00:00:00 2001 From: Nicola Vigano <nicola.vigano@esrf.fr> Date: Tue, 24 Feb 2015 19:44:14 +0100 Subject: [PATCH] 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 | 2 +- .../gt6DCreateProjDataFromGrainCluster.m | 2 +- zUtil_Deformation/gt6DPlotReconstructionODF.m | 83 +++----- .../gt6DMakeClusterIPFFigureForPaper.m | 139 +++++++++++++ .../gt6DMakeClusterKAMFigureForPaper.m | 182 ++++++++++++++++++ .../gt6DMakeClusterODFFigureForPaper.m | 143 ++++++++++++++ 6 files changed, 494 insertions(+), 57 deletions(-) create mode 100644 zUtil_Deformation/plotting/gt6DMakeClusterIPFFigureForPaper.m create mode 100644 zUtil_Deformation/plotting/gt6DMakeClusterKAMFigureForPaper.m create mode 100644 zUtil_Deformation/plotting/gt6DMakeClusterODFFigureForPaper.m diff --git a/4_grains/gtCalculateGrain.m b/4_grains/gtCalculateGrain.m index f15d02ee..53e5c1eb 100644 --- a/4_grains/gtCalculateGrain.m +++ b/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 diff --git a/zUtil_Deformation/gt6DCreateProjDataFromGrainCluster.m b/zUtil_Deformation/gt6DCreateProjDataFromGrainCluster.m index 2ad8a71a..2b4696df 100644 --- a/zUtil_Deformation/gt6DCreateProjDataFromGrainCluster.m +++ b/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); diff --git a/zUtil_Deformation/gt6DPlotReconstructionODF.m b/zUtil_Deformation/gt6DPlotReconstructionODF.m index 76dbf16b..928bf36f 100644 --- a/zUtil_Deformation/gt6DPlotReconstructionODF.m +++ b/zUtil_Deformation/gt6DPlotReconstructionODF.m @@ -1,13 +1,18 @@ -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; diff --git a/zUtil_Deformation/plotting/gt6DMakeClusterIPFFigureForPaper.m b/zUtil_Deformation/plotting/gt6DMakeClusterIPFFigureForPaper.m new file mode 100644 index 00000000..ac5168c7 --- /dev/null +++ b/zUtil_Deformation/plotting/gt6DMakeClusterIPFFigureForPaper.m @@ -0,0 +1,139 @@ +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 diff --git a/zUtil_Deformation/plotting/gt6DMakeClusterKAMFigureForPaper.m b/zUtil_Deformation/plotting/gt6DMakeClusterKAMFigureForPaper.m new file mode 100644 index 00000000..bb684147 --- /dev/null +++ b/zUtil_Deformation/plotting/gt6DMakeClusterKAMFigureForPaper.m @@ -0,0 +1,182 @@ +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 diff --git a/zUtil_Deformation/plotting/gt6DMakeClusterODFFigureForPaper.m b/zUtil_Deformation/plotting/gt6DMakeClusterODFFigureForPaper.m new file mode 100644 index 00000000..c5962d3e --- /dev/null +++ b/zUtil_Deformation/plotting/gt6DMakeClusterODFFigureForPaper.m @@ -0,0 +1,143 @@ +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 -- GitLab