DEFINE VOXELS define edges of a voxel grid covering a model [XVEC, YVEC, ZVEC, OPT] = DEFINE_VOXELS(FMDL) [XVEC, YVEC, ZVEC, OPT] = DEFINE_VOXELS(FMDL, OPT) Inputs: FMDL = an EIDORS forward model structure OPT = an option structure with the following fields and defaults: opt.imgsz = [32 32 32]; % X, Y and Z dimensions of the voxel grid. % If opt.zvec is provided, Z need not be % specified. If imgsz is a scalar, X = Y = Z % and cube_voxels defaults to 1; 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.xlim = [] % 2-element vector specifying the minimum and % maximum extent in the X direction. % Default: extent of FMDL.nodes opt.ylim = [] % 2-element vector specifying the minimum and % maximum extent in the X direction. % Default: extent of FMDL.nodes opt.zlim = [] % 2-element vector specifying the minimum and % maximum extent in the X direction. % Default: extent of FMDL.nodes opt.square_pixels = 1; % adjust imgsz to get square pixels (in XY) % Default: 1 if imgsz is scalar, 0 otherwise opt.cube_voxels = 1; % adjust imgsz to get cube voxels (in XYZ) % Default: 1 if imgsz is scalar, 0 otherwise Outputs: XVEC, YVEC, ZVEC = edges of voxels suitable for use with NDGRID OPT = an option structure as above with all fields present and consistent, and the additional fields x_pts, y_pts, and z_pts defining cordinates of voxel centers (to be passed to mdl_slice_mapper) DEFINE_VOXELS(FMDL, 'Key', Value, ...) is an alternative way to specify options. Example: [xvec, yvec, zvec] = define_voxels(img.fwd_model); [~,~,~,opt] = define_voxels(fmdl,'imgsz',32,'zvec',[0,1,2,3]); [~,~,~,opt] = define_voxels(fmdl,'imgsz',32, ... 'cube_voxels', 0, 'square_pixels', 0); See also MK_VOXEL_VOLUME, CALC_VOXELS
0001 function [xvec, yvec, zvec, opt] = define_voxels(fmdl, varargin) 0002 %DEFINE VOXELS define edges of a voxel grid covering a model 0003 % [XVEC, YVEC, ZVEC, OPT] = DEFINE_VOXELS(FMDL) 0004 % [XVEC, YVEC, ZVEC, OPT] = DEFINE_VOXELS(FMDL, OPT) 0005 % 0006 % Inputs: 0007 % FMDL = an EIDORS forward model structure 0008 % OPT = an option structure with the following fields and defaults: 0009 % opt.imgsz = [32 32 32]; % X, Y and Z dimensions of the voxel grid. 0010 % % If opt.zvec is provided, Z need not be 0011 % % specified. If imgsz is a scalar, X = Y = Z 0012 % % and cube_voxels defaults to 1; 0013 % opt.xvec = [] % Specific X cut-planes between voxels 0014 % % A scalar means the number of planes 0015 % % Takes precedence over other options 0016 % opt.yvec = [] % Specific Y cut-planes between voxels 0017 % % A scalar means the number of planes 0018 % % Takes precedence over other options 0019 % opt.zvec = [] % Specific Z cut-planes between voxels 0020 % % A scalar means the number of planes 0021 % % Takes precedence over other options 0022 % opt.xlim = [] % 2-element vector specifying the minimum and 0023 % % maximum extent in the X direction. 0024 % % Default: extent of FMDL.nodes 0025 % opt.ylim = [] % 2-element vector specifying the minimum and 0026 % % maximum extent in the X direction. 0027 % % Default: extent of FMDL.nodes 0028 % opt.zlim = [] % 2-element vector specifying the minimum and 0029 % % maximum extent in the X direction. 0030 % % Default: extent of FMDL.nodes 0031 % opt.square_pixels = 1; % adjust imgsz to get square pixels (in XY) 0032 % % Default: 1 if imgsz is scalar, 0 otherwise 0033 % opt.cube_voxels = 1; % adjust imgsz to get cube voxels (in XYZ) 0034 % % Default: 1 if imgsz is scalar, 0 otherwise 0035 % 0036 % Outputs: 0037 % XVEC, YVEC, ZVEC = edges of voxels suitable for use with NDGRID 0038 % OPT = an option structure as above with all fields present 0039 % and consistent, and the additional fields 0040 % x_pts, y_pts, and z_pts defining cordinates of voxel 0041 % centers (to be passed to mdl_slice_mapper) 0042 % 0043 % DEFINE_VOXELS(FMDL, 'Key', Value, ...) is an alternative way to specify 0044 % options. 0045 % 0046 % Example: 0047 % [xvec, yvec, zvec] = define_voxels(img.fwd_model); 0048 % [~,~,~,opt] = define_voxels(fmdl,'imgsz',32,'zvec',[0,1,2,3]); 0049 % [~,~,~,opt] = define_voxels(fmdl,'imgsz',32, ... 0050 % 'cube_voxels', 0, 'square_pixels', 0); 0051 % See also MK_VOXEL_VOLUME, CALC_VOXELS 0052 0053 % (C) 2015-2024 Bartlomiej Grychtol - all rights reserved by Swisstom AG 0054 % License: GPL version 2 or 3 0055 % $Id: define_voxels.m 6877 2024-05-10 13:26:26Z bgrychtol $ 0056 0057 % >> SWISSTOM CONTRIBUTION << 0058 0059 switch fmdl.type 0060 case 'fwd_model' 0061 %nothing 0062 case 'image' 0063 fmdl = fmdl.fwd_model; 0064 otherwise 0065 error('EIDORS:WrongInput', 'Expected EIDORS image or fwd_model.') 0066 end 0067 0068 if nargin < 2 0069 opt = struct; 0070 elseif ischar(varargin{1}) 0071 opt = struct(varargin{:}); 0072 else 0073 opt = varargin{1}; 0074 end 0075 0076 if ~isfield(opt, 'imgsz') 0077 opt.imgsz = 32; 0078 end 0079 fields = {'xvec', 'yvec', 'zvec', 'xlim', 'ylim', 'zlim'}; 0080 for f = 1: numel(fields) 0081 if ~isfield(opt, fields{f}) 0082 opt.(fields{f}) = []; 0083 end 0084 end 0085 if isempty(opt.xvec) && isempty(opt.imgsz) 0086 error('EIDORS:WrongInput', 'opt.imgsz must not be empty if opt.xvec is empty or absent'); 0087 end 0088 if isscalar(opt.imgsz) 0089 opt.imgsz(2:3) = opt.imgsz; 0090 if ~isfield(opt, 'square_pixels') 0091 opt.square_pixels = 1; 0092 end 0093 if ~isfield(opt, 'cube_voxels') 0094 opt.cube_voxels = 1; 0095 end 0096 end 0097 if isempty(opt.zvec) && numel(opt.imgsz) == 2 0098 error('EIDORS:WrongInput', 'opt.imgsz must have 3 elements if opt.zvec is empty or absent'); 0099 end 0100 0101 if ~isfield(opt, 'square_pixels') 0102 opt.square_pixels = 0; 0103 end 0104 if ~isfield(opt, 'cube_voxels') 0105 opt.cube_voxels = 0; 0106 end 0107 0108 mingrid = min(fmdl.nodes); 0109 maxgrid = max(fmdl.nodes); 0110 fields = {'xlim', 'ylim', 'zlim'}; 0111 for f = 1:3 0112 if ~isempty(opt.(fields{f})) 0113 mingrid(f) = opt.(fields{f})(1); 0114 maxgrid(f) = opt.(fields{f})(2); 0115 else 0116 opt.(fields{f})(1) = mingrid(f); 0117 opt.(fields{f})(2) = maxgrid(f); 0118 end 0119 end 0120 0121 mdl_sz = maxgrid - mingrid; 0122 0123 allempty = isempty(opt.xvec) & isempty(opt.yvec) & isempty(opt.zvec); 0124 if opt.cube_voxels && ~allempty 0125 warning('EIDORS:IncompatibleOptions','Option cube_voxels is set to 0 when xvec, yvec or zvec is specifed'); 0126 opt.cube_voxels = 0; 0127 end 0128 if opt.cube_voxels && allempty 0129 side_sz = max(mdl_sz ./ opt.imgsz(1:3)); 0130 n_steps = ceil(mdl_sz / side_sz); 0131 mdl_ctr = mingrid + mdl_sz/2; 0132 mingrid = mdl_ctr - n_steps/2 * side_sz; 0133 maxgrid = mdl_ctr + n_steps/2 * side_sz; 0134 opt.imgsz = n_steps; 0135 0136 elseif opt.square_pixels 0137 if ~isempty(opt.xvec) || ~isempty(opt.yvec) 0138 warning('EIDORS:IncompatibleOptions','Option square_pixels is set to 0 when xvec or yvec is specifed'); 0139 opt.square_pixels = 0; 0140 else 0141 mdl_AR = mdl_sz(1)/mdl_sz(2); 0142 img_AR = opt.imgsz(1)/opt.imgsz(2); 0143 if mdl_AR < img_AR 0144 delta = (mdl_sz(2) * img_AR - mdl_sz(1)) /2; 0145 mingrid(1) = mingrid(1) - delta; 0146 maxgrid(1) = maxgrid(1) + delta; 0147 elseif mdl_AR > img_AR 0148 delta = (mdl_sz(1)/img_AR - mdl_sz(2)) / 2; 0149 mingrid(2) = mingrid(2) - delta; 0150 maxgrid(2) = maxgrid(2) + delta; 0151 end 0152 end 0153 end 0154 0155 fields = {'xvec', 'yvec', 'zvec'}; 0156 pts_fields = {'x_pts', 'y_pts', 'z_pts'}; 0157 for i = 1:3 0158 if isscalar(opt.(fields{i})) 0159 opt.imgsz(i) = opt.(fields{i}); % make imgsz consistent 0160 opt.(fields{i}) = []; 0161 end 0162 if isempty(opt.(fields{i})) 0163 opt.(fields{i}) = linspace(mingrid(i), maxgrid(i), opt.imgsz(i)+1); 0164 else 0165 opt.imgsz(i) = numel(opt.(fields{i})) - 1; % make imgsz consistent 0166 end 0167 opt.(pts_fields{i}) = opt.(fields{i})(1:end-1) + diff(opt.(fields{i}))/2; 0168 end 0169 0170 0171 xvec = opt.xvec; 0172 yvec = opt.yvec; 0173 zvec = opt.zvec; 0174