Description of mk_voxel

0001 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 5749 2018-05-16 18:39:36Z alistair_boyle $
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     x_pts = opt.xvec(1:end-1) + 0.5*diff(opt.xvec);
0109     y_pts = opt.yvec(1:end-1) + 0.5*diff(opt.yvec);
0110     z_pts = opt.zvec(1:end-1) + 0.5*diff(opt.zvec);
0111 
0112     rmdl.mdl_slice_mapper.x_pts = x_pts;
0113     rmdl.mdl_slice_mapper.y_pts = y_pts;
0114     rmdl.mdl_slice_mapper.z_pts = z_pts;
0115     
0116     
0117 %-------------------------------------------------------------------------%
0118 % Assemble a reference object for caching
0119 function cache_obj = get_cache_obj(fmdl, opt)
0120     tmp = struct;
0121     flds = {'nodes','elems'};
0122     for f = flds;
0123         tmp.(f{1}) = fmdl.(f{1});
0124     end
0125     cache_obj = {tmp, opt};
0126     
0127 
0128 %-------------------------------------------------------------------------%
0129 % Parse option struct
0130  function opt = parse_opts(fmdl, opt)
0131     if ~isfield(opt, 'imgsz'),     
0132         opt.imgsz = [32 32 4]; 
0133     end
0134     if ~isfield(opt, 'square_pixels')
0135         opt.square_pixels = 0;
0136     end
0137     if ~isfield(opt, 'cube_voxels')
0138         opt.cube_voxels = 0;
0139     end
0140     if ~isfield(opt, 'xvec')
0141         opt.xvec = [];
0142     end
0143     if ~isfield(opt, 'yvec')
0144         opt.yvec = [];
0145     end
0146     if ~isfield(opt, 'zvec')
0147         opt.zvec = [];
0148     end
0149     if ~isfield(opt, 'prune_model')
0150         opt.prune_model = true;
0151     end
0152     if isempty(opt.xvec) && isempty(opt.imgsz)
0153         error('EIDORS:WrongInput', 'opt.imgsz must not be empty if opt.xvec is empty or absent');
0154     end
0155     if isempty(opt.yvec) && numel(opt.imgsz) < 2
0156         error('EIDORS:WrongInput', 'opt.imgsz must have at least 2 elements if opt.yvec is empty or absent');
0157     end
0158     if isempty(opt.zvec) && numel(opt.imgsz) < 3
0159         error('EIDORS:WrongInput', 'opt.imgsz must have 3 elements if opt.zvec is empty or absent');
0160     end
0161     
0162     mingrid = min(fmdl.nodes);
0163     maxgrid = max(fmdl.nodes);
0164     mdl_sz = maxgrid - mingrid;
0165     
0166     allempty = isempty(opt.xvec) & isempty(opt.yvec) & isempty(opt.zvec);
0167     if opt.cube_voxels && ~allempty
0168         warning('EIDORS:IncompatibleOptions','Option cube_voxels is ignored when xvec, yvec or zvec is specifed');
0169     end
0170     if opt.cube_voxels && allempty
0171         side_sz = max(mdl_sz ./ opt.imgsz(1:3));
0172         n_steps = ceil(mdl_sz / side_sz);
0173         mdl_ctr = mingrid + mdl_sz/2;
0174         mingrid = mdl_ctr - n_steps/2 * side_sz;
0175         maxgrid = mdl_ctr + n_steps/2 * side_sz;
0176         opt.imgsz = n_steps;
0177  
0178     elseif opt.square_pixels
0179         if ~isempty(opt.xvec) || ~isempty(opt.yvec)
0180             warning('EIDORS:IncompatibleOptions','Option square_pixels is ignored when xvec or yvec is specifed');
0181         else
0182             mdl_AR = mdl_sz(1)/mdl_sz(2);
0183             img_AR = opt.imgsz(1)/opt.imgsz(2);
0184             if mdl_AR < img_AR
0185                 delta = (mdl_sz(2) * img_AR - mdl_sz(1)) /2;
0186                 mingrid(1) = mingrid(1) - delta;
0187                 maxgrid(1) = maxgrid(1) + delta;
0188             elseif mdl_AR > img_AR
0189                 delta = (mdl_sz(1)/img_AR - mdl_sz(2)) / 2;
0190                 mingrid(2) = mingrid(2) - delta;
0191                 maxgrid(2) = maxgrid(2) + delta;
0192             end
0193         end
0194     end
0195     if isempty(opt.xvec)
0196         opt.xvec = linspace(mingrid(1), maxgrid(1), opt.imgsz(1)+1);
0197     end      
0198     if isempty(opt.yvec)
0199         opt.yvec = linspace(mingrid(2), maxgrid(2), opt.imgsz(2)+1);
0200     end
0201     if isempty(opt.zvec)
0202         opt.zvec = linspace(mingrid(3), maxgrid(3), opt.imgsz(3)+1);
0203     end
0204     
0205     if ~isfield(opt, 'do_coarse2fine')
0206         opt.do_coarse2fine = 1;
0207     end
0208     if ~isfield(opt, 'z_depth')
0209         opt.z_depth = inf;
0210     end
0211 
0212 %-------------------------------------------------------------------------%
0213 % Perfom unit tests
0214 function do_unit_test
0215 
0216 fmdl = mk_library_model('adult_male_16el');
0217 % fmdl= ng_mk_cyl_models([2,2,.4],[16,1],[.1,0,.025]);
0218 opt.square_pixels = 1;
0219 [rmdl, fmdl] = mk_voxel_volume(fmdl, opt);
0220 
0221 
0222 subplot(121)
0223 rimg = mk_image(rmdl,0);
0224 rimg.elem_data = zeros(size(rmdl.coarse2fine,2),1);
0225 idx = round(rand(5,1)* length(rimg.elem_data));
0226 rimg.elem_data(idx) = 1;
0227 show_fem(rimg);
0228 
0229 subplot(122)
0230 img = mk_image(fmdl,0);
0231 img.elem_data = fmdl.coarse2fine * rimg.elem_data;
0232 img.calc_colours.ref_level = 0;
0233 img.calc_colours.transparency_thresh = 1e-2;
0234 
0235 show_fem(img);
0236 
0237 
0238 
0239 
0240
mk_voxel_volume

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
