MK_VOXEL_VOLUME create a voxel model to reconstruct on
OUT = MK_VOXEL_VOLUME(MDL)
OUT = MK_VOXEL_VOLUME(MDL, OPT)
Inputs:
MDL = an EIDORS forward or inverse model structure
OPT = an option structure with the following fields and defaults:
opt.imgsz = [32 32 4]; % X, Y and Z dimensions of the voxel grid
opt.xvec = [] % Specific X cut-planes between voxels
% A scalar means the number of planes
% Takes precedence over other options
opt.yvec = [] % Specific Y cut-planes between voxels
% A scalar means the number of planes
% Takes precedence over other options
opt.zvec = [] % Specific Z cut-planes between voxels
% A scalar means the number of planes
% Takes precedence over other options
opt.square_pixels = 0; % adjust imgsz to get square pixels (in XY)
opt.cube_voxels = 0; % adjust imgsz to get cube voxels (in XYZ)
opt.prune_model = true % removes voxels outside the supplied MDL
% This runs mk_grid_c2f. For simple
% geometries, such a cylinder, it is much
% quicker to set to false and prune manually.
opt.save_memory % passed to mk_grid_c2f
Output depends on the type of model supplied. If MDL is a fwd_model
structure then OUT is a rec_model. If MDL is an inv_model, then OUT is a
modified version of it, with the voxel volume in inv_model.rec_model and
updated inv_model.fwd_model.coarse2fine
[OUT FMDL] = MK_VOXEL_VOLUME(MDL, ...) also returns the forward model
structure with the coarse2fine field.
See also MK_PIXEL_SLICE, MK_GRID_MODEL, MK_GRID_C2F0001 function [imdl, fmdl] = mk_voxel_volume(varargin) 0002 %MK_VOXEL_VOLUME create a voxel model to reconstruct on 0003 % OUT = MK_VOXEL_VOLUME(MDL) 0004 % OUT = MK_VOXEL_VOLUME(MDL, OPT) 0005 % 0006 % Inputs: 0007 % MDL = an EIDORS forward or inverse model structure 0008 % OPT = an option structure with the following fields and defaults: 0009 % opt.imgsz = [32 32 4]; % X, Y and Z dimensions of the voxel grid 0010 % opt.xvec = [] % Specific X cut-planes between voxels 0011 % % A scalar means the number of planes 0012 % % Takes precedence over other options 0013 % opt.yvec = [] % Specific Y cut-planes between voxels 0014 % % A scalar means the number of planes 0015 % % Takes precedence over other options 0016 % opt.zvec = [] % Specific Z cut-planes between voxels 0017 % % A scalar means the number of planes 0018 % % Takes precedence over other options 0019 % opt.square_pixels = 0; % adjust imgsz to get square pixels (in XY) 0020 % opt.cube_voxels = 0; % adjust imgsz to get cube voxels (in XYZ) 0021 % opt.prune_model = true % removes voxels outside the supplied MDL 0022 % % This runs mk_grid_c2f. For simple 0023 % % geometries, such a cylinder, it is much 0024 % % quicker to set to false and prune manually. 0025 % opt.save_memory % passed to mk_grid_c2f 0026 % 0027 % Output depends on the type of model supplied. If MDL is a fwd_model 0028 % structure then OUT is a rec_model. If MDL is an inv_model, then OUT is a 0029 % modified version of it, with the voxel volume in inv_model.rec_model and 0030 % updated inv_model.fwd_model.coarse2fine 0031 % 0032 % [OUT FMDL] = MK_VOXEL_VOLUME(MDL, ...) also returns the forward model 0033 % structure with the coarse2fine field. 0034 % 0035 % See also MK_PIXEL_SLICE, MK_GRID_MODEL, MK_GRID_C2F 0036 0037 % (C) 2015 Bartlomiej Grychtol - all rights reserved by Swisstom AG 0038 % License: GPL version 2 or 3 0039 % $Id: mk_voxel_volume.m 6809 2024-04-23 14:14:43Z aadler $ 0040 0041 % >> SWISSTOM CONTRIBUTION << 0042 0043 imdl = varargin{1}; 0044 % if input is 'UNIT_TEST', run tests 0045 if ischar(imdl) && strcmp(imdl,'UNIT_TEST') 0046 do_unit_test; clear imdl 0047 return; 0048 end 0049 0050 if nargin < 2 0051 opt = struct; 0052 else 0053 opt = varargin{2}; 0054 end 0055 0056 switch(imdl.type) 0057 case 'inv_model' 0058 fmdl = imdl.fwd_model; 0059 case 'fwd_model' 0060 fmdl = imdl; 0061 otherwise 0062 error('An EIDORS inverse or forward model struct required'); 0063 end 0064 0065 opt = parse_opts(fmdl,opt); 0066 0067 copt.fstr = 'mk_voxel_volume'; 0068 copt.cache_obj = get_cache_obj(fmdl, opt); 0069 [rmdl, c2f] = eidors_cache(@do_voxel_volume,{fmdl, opt},copt); 0070 0071 if ~isempty(c2f) 0072 fmdl.coarse2fine = c2f; 0073 end 0074 0075 switch imdl.type 0076 case 'inv_model' 0077 imdl.rec_model = rmdl; 0078 imdl.fwd_model = fmdl; 0079 case 'fwd_model' 0080 imdl = rmdl; 0081 end 0082 0083 0084 %-------------------------------------------------------------------------% 0085 % The main function 0086 function [rmdl, c2f] = do_voxel_volume(fmdl,opt) 0087 0088 rmdl = mk_grid_model([],opt.xvec,opt.yvec,opt.zvec); 0089 0090 0091 c2f = []; 0092 if ~opt.prune_model, return, end 0093 0094 % fmdl.elems = fmdl.elems( 210714,:); 0095 [c2f, m] = mk_grid_c2f(fmdl, rmdl, opt); 0096 0097 inside = any(c2f,1); 0098 c2f(:,~inside) = []; 0099 rm = ~logical(rmdl.coarse2fine*inside'); 0100 rmdl.elems(rm,:) = []; 0101 rmdl.coarse2fine(rm,:) = []; 0102 rmdl.coarse2fine(:,~inside) = []; 0103 0104 0105 bnd_fcs = ones(1,nnz(inside))*m.vox2face(inside,:) == 1; 0106 rmdl.boundary = m.faces(bnd_fcs,:); 0107 rmdl.inside = inside; % the inside array is useful in other functions 0108 0109 rmdl.mdl_slice_mapper.x_pts = opt.x_pts; 0110 rmdl.mdl_slice_mapper.y_pts = opt.y_pts; 0111 rmdl.mdl_slice_mapper.z_pts = opt.z_pts; 0112 0113 0114 %-------------------------------------------------------------------------% 0115 % Assemble a reference object for caching 0116 function cache_obj = get_cache_obj(fmdl, opt) 0117 tmp = struct; 0118 flds = {'nodes','elems'}; 0119 for f = flds; 0120 tmp.(f{1}) = fmdl.(f{1}); 0121 end 0122 cache_obj = {tmp, opt}; 0123 0124 0125 %-------------------------------------------------------------------------% 0126 % Parse option struct 0127 function opt = parse_opts(fmdl, opt) 0128 if ~isfield(opt, 'imgsz') 0129 opt.imgsz = [32 32 4]; 0130 end 0131 if ~isfield(opt, 'square_pixels') 0132 opt.square_pixels = 0; 0133 end 0134 if ~isfield(opt, 'cube_voxels') 0135 opt.cube_voxels = 0; 0136 end 0137 if ~isfield(opt, 'xvec') 0138 opt.xvec = []; 0139 end 0140 if ~isfield(opt, 'yvec') 0141 opt.yvec = []; 0142 end 0143 if ~isfield(opt, 'zvec') 0144 opt.zvec = []; 0145 end 0146 if ~isfield(opt, 'prune_model') 0147 opt.prune_model = true; 0148 end 0149 0150 [~,~,~,opt] = define_voxels(fmdl, opt); 0151 0152 if ~isfield(opt, 'do_coarse2fine') 0153 opt.do_coarse2fine = 1; 0154 end 0155 if ~isfield(opt, 'z_depth') 0156 opt.z_depth = inf; 0157 end 0158 0159 %-------------------------------------------------------------------------% 0160 % Perfom unit tests 0161 function do_unit_test 0162 % Check for convhull errors 0163 vopt.imgsz = [10,10,10]; 0164 imdl = mk_common_model('n3r2',[16,2]); 0165 [img.fwd_model,cmdl] = mk_voxel_volume(imdl.fwd_model,vopt); 0166 unit_test_cmp('Simple',cmdl.elems(3,:),[64,65,96,33]); 0167 0168 0169 fmdl = mk_library_model('adult_male_16el'); 0170 % fmdl= ng_mk_cyl_models([2,2,.4],[16,1],[.1,0,.025]); 0171 opt.square_pixels = 1; 0172 [rmdl, fmdl] = mk_voxel_volume(fmdl, opt); 0173 0174 0175 subplot(121) 0176 rimg = mk_image(rmdl,0); 0177 rimg.elem_data = zeros(size(rmdl.coarse2fine,2),1); 0178 idx = round(rand(5,1)* length(rimg.elem_data)); 0179 rimg.elem_data(idx) = 1; 0180 show_fem(rimg); 0181 0182 subplot(122) 0183 img = mk_image(fmdl,0); 0184 img.elem_data = fmdl.coarse2fine * rimg.elem_data; 0185 img.calc_colours.ref_level = 0; 0186 img.calc_colours.transparency_thresh = 1e-2; 0187 0188 show_fem(img); 0189 0190 0191 0192 0193