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 7070 2024-12-09 17:21:10Z 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
0040    opt = struct; 
0041 else
0042    opt.level = level;   
0043 end
0044 
0045 opt = parse_opts(fmdl,opt);
0046 
0047 [rmdl fmdl] = eidors_cache(@do_pixel_slice,{fmdl, opt},'mk_pixel_slice');
0048 
0049 switch imdl.type
0050    case 'inv_model'
0051       imdl.rec_model = rmdl;
0052       imdl.fwd_model = fmdl;
0053    case 'fwd_model'
0054       imdl = rmdl;
0055 end
0056 
0057 function [rmdl fmdl] = do_pixel_slice(fmdl, opt);
0058 [tmp, ~, ~, ifun] =level_model_slice(fmdl,opt.level);
0059 slc = mdl_slice_mesher(tmp,[inf inf 0]);
0060 slc = slc.fwd_model;
0061 mingrid = min(slc.nodes);
0062 maxgrid = max(slc.nodes);
0063 bnd = find_boundary(slc);
0064 % contour_boundary = order_loop(slc.nodes(unique(bnd),:));
0065 
0066 if opt.square_pixels ==1
0067     mdl_sz = maxgrid - mingrid;
0068     mdl_AR = mdl_sz(1)/mdl_sz(2);
0069     img_AR = opt.imgsz(1)/opt.imgsz(2);
0070     if mdl_AR < img_AR
0071         delta = (mdl_sz(2) * img_AR - mdl_sz(1)) /2;
0072         mingrid(1) = mingrid(1) - delta;
0073         maxgrid(1) = maxgrid(1) + delta;
0074     elseif mdl_AR > img_AR
0075         delta = (mdl_sz(1)/img_AR - mdl_sz(2)) / 2;
0076         mingrid(2) = mingrid(2) - delta;
0077         maxgrid(2) = maxgrid(2) + delta;
0078     end
0079 end
0080 
0081 xgrid = linspace(mingrid(1),maxgrid(1),opt.imgsz(1)+1);
0082 ygrid = linspace(mingrid(2),maxgrid(2),opt.imgsz(2)+1);
0083 rmdl = mk_grid_model([],xgrid,ygrid);
0084 x_pts = xgrid(1:end-1) + 0.5*diff(xgrid);
0085 y_pts = ygrid(1:end-1) + 0.5*diff(ygrid);
0086 % y_pts = fliplr(y_pts); %medical
0087 
0088 % NOTE: This controls the image resolution. If you want higher res, you
0089 % need to either specify it in opt.imgsz or manually overwrite (or remove)
0090 % the imdl.rec_model.mdl_slice_mapper.
0091 rmdl.mdl_slice_mapper.x_pts = x_pts;
0092 rmdl.mdl_slice_mapper.y_pts = y_pts;
0093 rmdl.mdl_slice_mapper.level = opt.level;
0094 rmdl.mdl_slice_mapper.model_2d = 1;
0095 x_avg = conv2(xgrid, [1,1]/2,'valid');
0096 y_avg = conv2(ygrid, [1,1]/2,'valid');
0097 [x,y] = ndgrid( x_avg, y_avg);
0098 
0099 % 20141119: The inpolygon approach fails on non-simply-connected domains
0100 % inside = inpolygon(x(:),y(:),contour_boundary(:,1),contour_boundary(:,2) );
0101 P = [x(:) y(:)]; % P(end,3) = 0;
0102 inside = any(point_in_triangle(P, slc.elems, slc.nodes(:,1:2)),2);
0103 % inside = full(inside);
0104 
0105 ff = find(~inside);
0106 
0107 if opt.do_coarse2fine
0108 %     rmdl.mk_coarse_fine_mapping.z_depth = opt.z_depth;
0109 %     fmdl.coarse2fine = mk_coarse_fine_mapping(tmp,rmdl);
0110 %     fmdl.coarse2fine(:,ff) = [];
0111     if isinf(opt.z_depth)
0112        zgrid = [min(tmp.nodes(:,3))-1 max(tmp.nodes(:,3))+1];
0113     else
0114        zgrid = [-opt.z_depth opt.z_depth];
0115     end
0116     [~, c2f] = mk_grid_model(tmp,xgrid,ygrid,zgrid);
0117     c2f(:,ff) = [];
0118     fmdl = eidors_obj('set',fmdl,'coarse2fine', c2f);
0119 end
0120 
0121 rmdl.elems([2*ff, 2*ff-1],:)= [];
0122 rmdl.coarse2fine([2*ff, 2*ff-1],:)= [];
0123 rmdl.coarse2fine(:,ff)= [];
0124 % rmdl.boundary = find_boundary(rmdl);
0125 % show individual elements (more like how the 2d grid models display)
0126 rmdl.boundary = rmdl.elems;
0127 rmdl.inside   = inside; % the inside array is useful in other functions
0128 
0129 
0130 if isfield(fmdl,'mat_idx')
0131    rmdl.mat_idx = calc_mat_idx(rmdl,fmdl,ff,opt);
0132 end
0133 
0134 % map electrodes
0135 rmdl.nodes(:,3) = 0;
0136 rmdl.nodes =  ifun(rmdl.nodes);
0137 slc.nodes =  ifun(slc.nodes);
0138 
0139 % no longer needed, level_model_slice does this
0140 % if ~isstruct(opt.level)
0141 %    isf = ~isinf(opt.level);
0142 %    if nnz(isf) == 1
0143 %       rmdl.nodes(:,isf) = opt.level(:,isf);
0144 %    end
0145 % end
0146 
0147 if isfield(slc, 'electrode')
0148    for i = flipud(1:numel(slc.electrode))
0149         tmp = rmfield(slc.electrode(i), 'nodes');
0150         x_elec = slc.nodes( [slc.electrode(i).nodes], 1);
0151         y_elec = slc.nodes( [slc.electrode(i).nodes], 2);
0152         z_elec = slc.nodes( [slc.electrode(i).nodes], 3);
0153         tmp.pos       = [x_elec, y_elec, z_elec];
0154         elec(i) = tmp;
0155     end
0156     rmdl.electrode = elec;
0157 end
0158    
0159 % Not needed any more (uses fwd_model.mdl_slice_mapper)
0160 % rmdl.show_slices.levels = opt.level;
0161 
0162 % not quite correct, but 'rec_model' has many side effects
0163 rmdl = eidors_obj('fwd_model', rmdl); 
0164       
0165 
0166 function mat_idx = calc_mat_idx(rmdl,fmdl,ff,opt)
0167    % calculate mat_idx for the rec_model
0168    fmdl.mdl_slice_mapper = rmfield(rmdl.mdl_slice_mapper,'model_2d');
0169    fmdl.mdl_slice_mapper.level = opt.level;
0170    img = mk_image(fmdl,0);
0171    for i = 1:length(fmdl.mat_idx);
0172       img.elem_data(fmdl.mat_idx{i}) = i;
0173    end
0174    slice = calc_slices(img); % level specified in fmdl.mdl_slice_mapper
0175    slice = slice';
0176    mat = reshape([slice(:)'; slice(:)'],1,[]);
0177    mat([2*ff, 2*ff-1])= [];
0178    mat_idx = cell(max(mat),1);
0179    for i = 1:max(mat)
0180       mat_idx(i) = {find(mat==i)'};
0181    end
0182 
0183 
0184  function opt = parse_opts(fmdl, opt)
0185     if ~isfield(opt, 'imgsz')     
0186         opt.imgsz = [32 32]; 
0187     end
0188     if ~isfield(opt, 'square_pixels')
0189         opt.square_pixels = 0;
0190     end
0191     if ~isfield(opt, 'level') || isempty(opt.level)
0192         opt.level = get_elec_level(fmdl);
0193     else
0194         if ~isstruct(opt.level) && numel(opt.level) ==1
0195             opt.level = [inf inf opt.level];
0196         end
0197     end
0198     if ~isfield(opt, 'do_coarse2fine')
0199         opt.do_coarse2fine = 1;
0200     end
0201     if ~isfield(opt, 'z_depth')
0202         opt.z_depth = inf;
0203     end
0204     
0205 function elec_lev = get_elec_level(fmdl)
0206     z_elec= fmdl.nodes( [fmdl.electrode(:).nodes], 3);
0207     min_e = min(z_elec); max_e = max(z_elec);
0208     elec_lev = [inf,inf,mean([min_e,max_e])];
0209 
0210     
0211 function do_unit_test
0212     imdl = mk_common_model('n3r2',[16,2]);
0213     opt.square_pixels = 1;
0214     opt.imgsz = [16 16];
0215     mdl = mk_pixel_slice(imdl.fwd_model,[inf 2 2.5], opt);
0216     img = mk_image(mdl,1);
0217     
0218     subplot(231)
0219     show_fem(imdl.fwd_model);
0220     view([-50 10])
0221 
0222     subplot(232)
0223     show_fem(img);
0224     zlim([0 3]);
0225     ylim([-1 1])
0226     xlim([-1 1]);
0227     view([-50 10])
0228     
0229     subplot(233)
0230     show_slices(img);
0231     
0232     subplot(234)
0233     imdl = mk_pixel_slice(imdl);
0234     img = mk_image(imdl.rec_model,1);
0235     show_fem(img);
0236     zlim([0 3]);
0237     ylim([-1 1])
0238     xlim([-1 1]);
0239     view([-50 10])
0240     
0241     subplot(235)
0242     mdl = mk_library_model('pig_23kg_16el_lungs');
0243     img = mk_image(mdl,1);
0244     for i = 1:length(mdl.mat_idx)
0245        img.elem_data(mdl.mat_idx{i}) = i;
0246     end
0247     show_fem(img)
0248     view(2)
0249     
0250     subplot(236)
0251     clear opt
0252     opt.imgsz = [64 64];
0253     opt.square_pixels = 1;
0254     opt.do_coarse2fine = 0;
0255     mdl = mk_library_model('pig_23kg_16el_lungs');
0256 %     rmdl = mk_pixel_slice(mdl,opt); % deprecated
0257     rmdl = mk_pixel_slice(mdl,[],opt);
0258     img = mk_image(rmdl,1);
0259     for i = 1:length(rmdl.mat_idx)
0260        img.elem_data(rmdl.mat_idx{i}) = i;
0261     end
0262     show_slices(img);