CALC_VOXELS Calculate volumetric data from an image
V = CALC_VOXELS(IMG) calculates volumetric representation of the image
IMG using a 3D grid with a maximum dimension of 32. V is a 3-D or 4-D
matrix with dimensions Y x X x Z x N, where N is the number of frames in
the image data (time), and X, Y, Z refer to the axes of the image node
space.
V = CALC_VOXELS(IMG, OPT) allows specifying an options struct:
OPT.method - specifies the method of determining value in each voxel:
'centre' - the value at its centre (default)
'average' - the average value (expensive)
The average method can only be used for images with
values specified on elements (elem_data). It produces
smoother images that simulate the partial volume effect
and look better at lower resolution than the 'centre'
method.
OPT.imgsz - specifies the target grid size (default [32 32 32]).
See DEFINE_VOXELS for details, alternatives and
additional options.
voxels supported by DEFINE_VOXELS can also be used.
OPT.save_memory - only useful if method is 'average', reduces memory
consumption, at the expense of time. See MK_GRID_C2F.
All other options of MK_GRID_C2F may also be specified.
NOTE: The OPT struct is passed in its entirety to both DEFINE_VOXELS and
MK_GRID_C2F.
V = CALC_VOXELS(IMG, OPT, LEVEL) specifies the orientaion of the XY plane
of the voxel grid. Any definition of a single slice accepted by
LEVEL_MODEL_SLICE can be used. OPT can be an empty struct.
[V, RIMG] = CALC_VOXELS(_) also returns a grid-based image as produced
by MK_GRID_MODEL, suitable for use with SHOW_FEM.
CALC_VOXELS supports parametrization through DATA_MAPPER.
Example:
fmdl = ng_mk_ellip_models([1,1.5,0.8,0.1],[],[]);
fmdl.nodes(:,3) = fmdl.nodes(:,3) - .25;
fmdl = fix_model(fmdl, struct('elem_centre',true));
img = mk_image(fmdl, prod(fmdl.elem_centre,2));
[V, rimg] = calc_voxels(img);
subplot(221), show_fem(img)
subplot(222), show_fem(rimg);
subplot(223), show_slices(rimg,3);
subplot(224), show_slices(V);
See also DEFINE_VOXELS, CALC_GRID_POINTS, MK_GRID_C2F, MK_GRID_MODEL,
CALC_SLICES, DATA_MAPPER0001 function [V, rimg] = calc_voxels(img, opt, level) 0002 %CALC_VOXELS Calculate volumetric data from an image 0003 % V = CALC_VOXELS(IMG) calculates volumetric representation of the image 0004 % IMG using a 3D grid with a maximum dimension of 32. V is a 3-D or 4-D 0005 % matrix with dimensions Y x X x Z x N, where N is the number of frames in 0006 % the image data (time), and X, Y, Z refer to the axes of the image node 0007 % space. 0008 % 0009 % V = CALC_VOXELS(IMG, OPT) allows specifying an options struct: 0010 % OPT.method - specifies the method of determining value in each voxel: 0011 % 'centre' - the value at its centre (default) 0012 % 'average' - the average value (expensive) 0013 % The average method can only be used for images with 0014 % values specified on elements (elem_data). It produces 0015 % smoother images that simulate the partial volume effect 0016 % and look better at lower resolution than the 'centre' 0017 % method. 0018 % OPT.imgsz - specifies the target grid size (default [32 32 32]). 0019 % See DEFINE_VOXELS for details, alternatives and 0020 % additional options. 0021 % voxels supported by DEFINE_VOXELS can also be used. 0022 % OPT.save_memory - only useful if method is 'average', reduces memory 0023 % consumption, at the expense of time. See MK_GRID_C2F. 0024 % All other options of MK_GRID_C2F may also be specified. 0025 % NOTE: The OPT struct is passed in its entirety to both DEFINE_VOXELS and 0026 % MK_GRID_C2F. 0027 % 0028 % V = CALC_VOXELS(IMG, OPT, LEVEL) specifies the orientaion of the XY plane 0029 % of the voxel grid. Any definition of a single slice accepted by 0030 % LEVEL_MODEL_SLICE can be used. OPT can be an empty struct. 0031 % 0032 % [V, RIMG] = CALC_VOXELS(_) also returns a grid-based image as produced 0033 % by MK_GRID_MODEL, suitable for use with SHOW_FEM. 0034 % 0035 % CALC_VOXELS supports parametrization through DATA_MAPPER. 0036 % 0037 % Example: 0038 % fmdl = ng_mk_ellip_models([1,1.5,0.8,0.1],[],[]); 0039 % fmdl.nodes(:,3) = fmdl.nodes(:,3) - .25; 0040 % fmdl = fix_model(fmdl, struct('elem_centre',true)); 0041 % img = mk_image(fmdl, prod(fmdl.elem_centre,2)); 0042 % [V, rimg] = calc_voxels(img); 0043 % subplot(221), show_fem(img) 0044 % subplot(222), show_fem(rimg); 0045 % subplot(223), show_slices(rimg,3); 0046 % subplot(224), show_slices(V); 0047 % 0048 % See also DEFINE_VOXELS, CALC_GRID_POINTS, MK_GRID_C2F, MK_GRID_MODEL, 0049 % CALC_SLICES, DATA_MAPPER 0050 0051 % (C) 2024 Bartlomiej Grychtol 0052 % License: GPL version 2 or 3 0053 % $Id: calc_voxels.m 7094 2024-12-20 23:05:22Z bgrychtol $ 0054 0055 if ischar(img) && strcmp(img, 'UNIT_TEST'), do_unit_test; return, end 0056 0057 if nargin < 2 0058 opt = struct(); 0059 end 0060 0061 if ~isfield(opt, 'method') 0062 opt.method = 'centre'; 0063 end 0064 0065 img = data_mapper(img); % support parametrization 0066 0067 do_rimg = false; 0068 0069 switch opt.method 0070 case 'centre' 0071 do_c2f = false; 0072 case 'average' 0073 do_rimg = true; 0074 do_c2f = true; 0075 otherwise 0076 error('EIDORS:WrongInput', 'Method must be centre or average.'); 0077 end 0078 0079 if isfield(img, 'node_data') && strcmp(opt.method, 'average') 0080 error('Only the centre method can be used with node data.'); 0081 end 0082 0083 if nargout > 1 0084 do_rimg = true; 0085 end 0086 0087 if nargin > 2 0088 img.fwd_model = level_model_slice(img.fwd_model, level); 0089 end 0090 0091 [~, ~, ~, opt] = define_voxels(img.fwd_model, opt); 0092 0093 if do_c2f % method 'average' 0094 rmdl = mk_grid_model([],opt.xvec,opt.yvec,opt.zvec); 0095 c2f = mk_grid_c2f(img.fwd_model, rmdl, opt); 0096 fvol = get_elem_volume(img.fwd_model,'no_c2f'); 0097 % need bsxfun because of octave, and it's super slow 0098 f2c = bsxfun(@times, c2f, fvol)'; 0099 %f2c = f2c ./ sum(f2c,2); % MATLAB:array:SizeLimitExceeded 0100 f2c = spdiag(1./sum(f2c,2)) * f2c; 0101 rimg = mk_image(rmdl,f2c*img.elem_data); 0102 V = calc_grid_points(rimg, opt.x_pts, opt.y_pts, opt.z_pts); 0103 else % method 'centre' 0104 V = calc_grid_points(img, opt.x_pts, opt.y_pts, opt.z_pts); 0105 if do_rimg 0106 rmdl = mk_grid_model([],opt.xvec,opt.yvec,opt.zvec); 0107 num_elems = size(rmdl.coarse2fine,2); 0108 rimg = mk_image(rmdl, reshape(V,num_elems,[])); 0109 end 0110 end 0111 0112 V = permute(V, [2 1 3 4]); 0113 0114 function do_unit_test 0115 do_cases 0116 do_example 0117 0118 function do_cases 0119 0120 imdl= mk_common_model('n3r2',[16,2]); 0121 fmdl = imdl.fwd_model; 0122 0123 obj01 = [390,391,393,396,402,478,479,480,484,486, ... 0124 664,665,666,667,668,670,671,672,676,677, ... 0125 678,755,760,761]; 0126 obj02 = [318,319,321,324,330,439,440,441,445,447, ... 0127 592,593,594,595,596,598,599,600,604,605, ... 0128 606,716,721,722]; 0129 0130 test = 1; 0131 while 1 0132 img = mk_image(fmdl, 1); 0133 opt = struct; 0134 error_expected = false; 0135 switch test 0136 case 1 0137 opt.imgsz = [32 32 32]; 0138 opt.cube_voxels = 1; 0139 img.elem_data(obj01) = 1.15; 0140 img.elem_data(obj02) = 0.8; 0141 expect_out = [22 22 32]; 0142 case 2 0143 opt.imgsz = [32 32 32]; 0144 opt.cube_voxels = 1; 0145 img.elem_data(:,2) = img.elem_data; 0146 img.elem_data(obj01,1) = 1.15; 0147 img.elem_data(obj02,2) = 0.8; 0148 expect_out = [22 22 32 2]; 0149 case 3 0150 opt.imgsz = [32 32 32]; 0151 opt.cube_voxels = 1; 0152 img.elem_data(:,2) = img.elem_data; 0153 img.elem_data(obj01,1) = 1.15; 0154 img.elem_data(obj02,2) = 0.8; 0155 opt.method = 'average'; 0156 expect_out = [22 22 32 2]; 0157 case 4 0158 opt.imgsz = [32 32 32]; 0159 opt.cube_voxels = 1; 0160 img = mk_image(fmdl, 1:num_nodes(fmdl)); 0161 expect_out = [22 22 32]; 0162 case 5 0163 opt.imgsz = [32 32 32]; 0164 opt.cube_voxels = 1; 0165 img = mk_image(fmdl, 1:num_nodes(fmdl)); 0166 opt.method = 'average'; 0167 error_expected = true; 0168 case 6 % empty options 0169 img.elem_data(obj01) = 1.15; 0170 img.elem_data(obj02) = 0.8; 0171 expect_out = [22 22 32]; 0172 otherwise 0173 break 0174 end 0175 test_txt = sprintf('calc_voxels T#%2d:', test); 0176 try 0177 [V, rimg] = calc_voxels(img, opt); 0178 unit_test_cmp(test_txt, size(V), expect_out); 0179 % fprintf('Test #%2d: %s\n', test, mat2str(size(V))); 0180 catch e 0181 if error_expected 0182 unit_test_cmp([test_txt, ' expected error'], error_expected,true); 0183 else 0184 fprintf('Test #%2d: error\n', test); 0185 end 0186 end 0187 0188 show_test(img, V, rimg) 0189 test = test + 1; 0190 end 0191 0192 0193 function show_test(img, V, rimg) 0194 s = warning('query','EIDORS:FirstImageOnly'); 0195 warning('off','EIDORS:FirstImageOnly'); 0196 try 0197 clf 0198 subplot(131) 0199 show_fem(img) 0200 subplot(132) 0201 show_fem(rimg) 0202 subplot(133) 0203 show_slices(V) 0204 catch e 0205 warning(s, 'EIDORS:FirstImageOnly'); 0206 rethrow(e) 0207 end 0208 0209 function do_example 0210 fmdl = ng_mk_ellip_models([1,1.5,0.8,0.1],[],[]); 0211 fmdl.nodes(:,3) = fmdl.nodes(:,3) - .25; 0212 fmdl = fix_model(fmdl, struct('elem_centre',true)); 0213 img = mk_image(fmdl, prod(fmdl.elem_centre,2)); 0214 [V, rimg] = calc_voxels(img); 0215 subplot(221), show_fem(img) 0216 subplot(222), show_fem(rimg); 0217 subplot(223), show_slices(rimg,3); 0218 subplot(224), show_slices(V);