mk_pixel_slice

MK_PIXEL_SLICE create a pixel model to reconstruct on
 OUT = MK_PIXEL_SLICE(MDL, LEVEL, OPT) creates a slice of pixels as a
 model to reconstruct on. 

 Inputs:
  MDL   = an EIDORS forward or inverse model structure
  LEVEL = any definition accepted by LEVEL_MODEL_SLICE for a single slice
  OPT   = an option structure with the following fields and defaults:
     opt.imgsz = [32 32];    % dimensions of the pixel grid
     opt.square_pixels = 0;  % adjust imgsz to get square pixels
     opt.do_coarse2fine = 1; % calcuate c2f on the forward model
     opt.z_depth = inf;      % z_depth to use with mk_coarse_fine_mapping

 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 pixel slice in inv_model.rec_model and
 updated inv_model.fwd_model.coarse2fine
 
 [OUT FMDL] = MK_PIXEL_SLICE(MDL, ...) also returns the forward model
 structure with the coarse2fine field.

 See also MK_COARSE_FINE_MAPPING, MK_GRID_MODEL, LEVEL_MODEL_SLICE

0001 function [imdl fmdl] = mk_pixel_slice(imdl,level,opt)
0002 %MK_PIXEL_SLICE create a pixel model to reconstruct on
0003 % OUT = MK_PIXEL_SLICE(MDL, LEVEL, OPT) creates a slice of pixels as a
0004 % model to reconstruct on.
0005 %
0006 % Inputs:
0007 %  MDL   = an EIDORS forward or inverse model structure
0008 %  LEVEL = any definition accepted by LEVEL_MODEL_SLICE for a single slice
0009 %  OPT   = an option structure with the following fields and defaults:
0010 %     opt.imgsz = [32 32];    % dimensions of the pixel grid
0011 %     opt.square_pixels = 0;  % adjust imgsz to get square pixels
0012 %     opt.do_coarse2fine = 1; % calcuate c2f on the forward model
0013 %     opt.z_depth = inf;      % z_depth to use with mk_coarse_fine_mapping
0014 %
0015 % Output depends on the type of model suplied. If MDL is a fwd_model
0016 % structure then OUT is a rec_model. If MDL is an inv_model, then OUT is a
0017 % modified version of it, with the pixel slice in inv_model.rec_model and
0018 % updated inv_model.fwd_model.coarse2fine
0019 %
0020 % [OUT FMDL] = MK_PIXEL_SLICE(MDL, ...) also returns the forward model
0021 % structure with the coarse2fine field.
0022 %
0023 % See also MK_COARSE_FINE_MAPPING, MK_GRID_MODEL, LEVEL_MODEL_SLICE
0024 
0025 % (C) 2013 Bartlomiej Grychtol. License: GPL version 2 or 3
0026 % $Id: mk_pixel_slice.m 6373 2022-05-13 16:46:17Z bgrychtol $
0027 
0028 if ischar(imdl) && strcmp(imdl,'UNIT_TEST'),do_unit_test;return;end;
0029 
0030 switch(imdl.type)
0031     case 'inv_model'
0032         fmdl = imdl.fwd_model;
0033     case 'fwd_model'
0034         fmdl = imdl;
0035     otherwise
0036         error('An EIDORS inverse or forward model struct required');
0037 end
0038 
0039 if nargin < 2, opt = struct; end
0040 if nargin > 1 
0041    if ~isstruct(level)
0042       opt.level = level;
0043    else
0044       opt = level;
0045    end
0046 end
0047 opt = parse_opts(fmdl,opt);
0048 
0049 [rmdl fmdl] = eidors_cache(@do_pixel_slice,{fmdl, opt},'mk_pixel_slice');
0050 
0051 switch imdl.type
0052    case 'inv_model'
0053       imdl.rec_model = rmdl;
0054       imdl.fwd_model = fmdl;
0055    case 'fwd_model'
0056       imdl = rmdl;
0057 end
0058 
0059 function [rmdl fmdl] = do_pixel_slice(fmdl, opt);
0060 [tmp, T, R] =level_model_slice(fmdl,opt.level);
0061 slc = mdl_slice_mesher(tmp,[inf inf 0]);
0062 slc = slc.fwd_model;
0063 mingrid = min(slc.nodes);
0064 maxgrid = max(slc.nodes);
0065 bnd = find_boundary(slc);
0066 % contour_boundary = order_loop(slc.nodes(unique(bnd),:));
0067 
0068 if opt.square_pixels ==1
0069     mdl_sz = maxgrid - mingrid;
0070     mdl_AR = mdl_sz(1)/mdl_sz(2);
0071     img_AR = opt.imgsz(1)/opt.imgsz(2);
0072     if mdl_AR < img_AR
0073         delta = (mdl_sz(2) * img_AR - mdl_sz(1)) /2;
0074         mingrid(1) = mingrid(1) - delta;
0075         maxgrid(1) = maxgrid(1) + delta;
0076     elseif mdl_AR > img_AR
0077         delta = (mdl_sz(1)/img_AR - mdl_sz(2)) / 2;
0078         mingrid(2) = mingrid(2) - delta;
0079         maxgrid(2) = maxgrid(2) + delta;
0080     end
0081 end
0082 
0083 xgrid = linspace(mingrid(1),maxgrid(1),opt.imgsz(1)+1);
0084 ygrid = linspace(mingrid(2),maxgrid(2),opt.imgsz(2)+1);
0085 rmdl = mk_grid_model([],xgrid,ygrid);
0086 x_pts = xgrid(1:end-1) + 0.5*diff(xgrid);
0087 y_pts = ygrid(1:end-1) + 0.5*diff(ygrid);
0088 % y_pts = fliplr(y_pts); %medical
0089 
0090 % NOTE: This controls the image resolution. If you want higher res, you
0091 % need to either specify it in opt.imgsz or manually overwrite (or remove)
0092 % the imdl.rec_model.mdl_slice_mapper.
0093 rmdl.mdl_slice_mapper.x_pts = x_pts;
0094 rmdl.mdl_slice_mapper.y_pts = y_pts;
0095 rmdl.mdl_slice_mapper.level = opt.level;
0096 rmdl.mdl_slice_mapper.model_2d = 1;
0097 x_avg = conv2(xgrid, [1,1]/2,'valid');
0098 y_avg = conv2(ygrid, [1,1]/2,'valid');
0099 [x,y] = ndgrid( x_avg, y_avg);
0100 
0101 % 20141119: The inpolygon approach fails on non-simply-connected domains
0102 % inside = inpolygon(x(:),y(:),contour_boundary(:,1),contour_boundary(:,2) );
0103 P = [x(:) y(:)]; % P(end,3) = 0;
0104 inside = any(point_in_triangle(P, slc.elems, slc.nodes(:,1:2)),2);
0105 % inside = full(inside);
0106 
0107 ff = find(~inside);
0108 
0109 if opt.do_coarse2fine
0110 %     rmdl.mk_coarse_fine_mapping.z_depth = opt.z_depth;
0111 %     fmdl.coarse2fine = mk_coarse_fine_mapping(tmp,rmdl);
0112 %     fmdl.coarse2fine(:,ff) = [];
0113     if isinf(opt.z_depth)
0114        zgrid = [min(tmp.nodes(:,3))-1 max(tmp.nodes(:,3))+1];
0115     else
0116        zgrid = [-opt.z_depth opt.z_depth];
0117     end
0118     [jnk, fmdl.coarse2fine] = mk_grid_model(tmp,xgrid,ygrid,zgrid);
0119     fmdl.coarse2fine(:,ff);
0120 end
0121 
0122 rmdl.elems([2*ff, 2*ff-1],:)= [];
0123 rmdl.coarse2fine([2*ff, 2*ff-1],:)= [];
0124 rmdl.coarse2fine(:,ff)= [];
0125 % rmdl.boundary = find_boundary(rmdl);
0126 % show individual elements (more like how the 2d grid models display)
0127 rmdl.boundary = rmdl.elems;
0128 rmdl.inside   = inside; % the inside array is useful in other functions
0129 
0130 
0131 if isfield(fmdl,'mat_idx')
0132    rmdl.mat_idx = calc_mat_idx(rmdl,fmdl,ff,opt);
0133 end
0134 
0135 % map electrodes
0136 rmdl.nodes(:,3) = 0;
0137 rmdl.nodes =  (R\rmdl.nodes' + T'*ones(1,length(rmdl.nodes)))';
0138 slc.nodes =  (R\slc.nodes' + T'*ones(1,length(slc.nodes)))';
0139 
0140 isf = ~isinf(opt.level);
0141 if nnz(isf) == 1
0142    rmdl.nodes(:,isf) = opt.level(:,isf);
0143 end
0144 
0145 if isfield(slc, 'electrode')
0146    for i = flipud(1:numel(slc.electrode))
0147         tmp = rmfield(slc.electrode(i), 'nodes');
0148         x_elec = slc.nodes( [slc.electrode(i).nodes], 1);
0149         y_elec = slc.nodes( [slc.electrode(i).nodes], 2);
0150         z_elec = slc.nodes( [slc.electrode(i).nodes], 3);
0151         tmp.pos       = [x_elec, y_elec, z_elec];
0152         elec(i) = tmp;
0153     end
0154     rmdl.electrode = elec;
0155 end
0156    
0157 
0158 rmdl.show_slices.levels = opt.level;
0159       
0160 
0161 function mat_idx = calc_mat_idx(rmdl,fmdl,ff,opt)
0162    % calculate mat_idx for the rec_model
0163    fmdl.mdl_slice_mapper = rmfield(rmdl.mdl_slice_mapper,'model_2d');
0164    fmdl.mdl_slice_mapper.level = opt.level;
0165    img = mk_image(fmdl,0);
0166    for i = 1:length(fmdl.mat_idx);
0167       img.elem_data(fmdl.mat_idx{i}) = i;
0168    end
0169    slice = calc_slices(img,opt.level);
0170    slice = slice';
0171    mat = reshape([slice(:)'; slice(:)'],1,[]);
0172    mat([2*ff, 2*ff-1])= [];
0173    mat_idx = cell(max(mat),1);
0174    for i = 1:max(mat)
0175       mat_idx(i) = {find(mat==i)'};
0176    end
0177 
0178 
0179  function opt = parse_opts(fmdl, opt)
0180     if ~isfield(opt, 'imgsz'),     
0181         opt.imgsz = [32 32]; 
0182     end
0183     if ~isfield(opt, 'square_pixels')
0184         opt.square_pixels = 0;
0185     end
0186     if ~isfield(opt, 'level')
0187         opt.level = get_elec_level(fmdl);
0188     else
0189         if numel(opt.level) ==1
0190             opt.level = [inf inf opt.level];
0191         end
0192     end
0193     if ~isfield(opt, 'do_coarse2fine')
0194         opt.do_coarse2fine = 1;
0195     end
0196     if ~isfield(opt, 'z_depth')
0197         opt.z_depth = inf;
0198     end
0199     
0200 function elec_lev = get_elec_level(fmdl)
0201     z_elec= fmdl.nodes( [fmdl.electrode(:).nodes], 3);
0202     min_e = min(z_elec); max_e = max(z_elec);
0203     elec_lev = [inf,inf,mean([min_e,max_e])];
0204 
0205     
0206 function do_unit_test
0207     imdl = mk_common_model('n3r2',[16,2]);
0208     opt.square_pixels = 1;
0209     opt.imgsz = [16 16];
0210     mdl = mk_pixel_slice(imdl.fwd_model,[inf 2 2.5], opt);
0211     img = mk_image(mdl,1);
0212     
0213     subplot(231)
0214     show_fem(imdl.fwd_model);
0215     view([-50 10])
0216 
0217     subplot(232)
0218     show_fem(img);
0219     zlim([0 3]);
0220     ylim([-1 1])
0221     xlim([-1 1]);
0222     view([-50 10])
0223     
0224     subplot(233)
0225     show_slices(img);
0226     
0227     subplot(234)
0228     imdl = mk_pixel_slice(imdl);
0229     img = mk_image(imdl.rec_model,1);
0230     show_fem(img);
0231     zlim([0 3]);
0232     ylim([-1 1])
0233     xlim([-1 1]);
0234     view([-50 10])
0235     
0236     subplot(235)
0237     mdl = mk_library_model('pig_23kg_16el_lungs');
0238     img = mk_image(mdl,1);
0239     for i = 1:length(mdl.mat_idx)
0240        img.elem_data(mdl.mat_idx{i}) = i;
0241     end
0242     show_fem(img)
0243     view(2)
0244     
0245     subplot(236)
0246     clear opt
0247     opt.imgsz = [64 64];
0248     opt.square_pixels = 1;
0249     opt.do_coarse2fine = 0;
0250     mdl = mk_library_model('pig_23kg_16el_lungs');
0251     rmdl = mk_pixel_slice(mdl,opt);
0252     img = mk_image(rmdl,1);
0253     for i = 1:length(rmdl.mat_idx)
0254        img.elem_data(rmdl.mat_idx{i}) = i;
0255     end
0256     show_slices(img);