mk_voxel_volume

PURPOSE ^

MK_VOXEL_VOLUME create a voxel model to reconstruct on

SYNOPSIS ^

function [imdl, fmdl] = mk_voxel_volume(varargin)

DESCRIPTION ^

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 suplied. 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_C2F

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005