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 = either a vector of x-,y-, and z-intercepts of the cut plane or
          a single value interpreted as the height of a single cut in the
          z plane (by default, a horizontal slice at the average electrode
          height will be created)
  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

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 = either a vector of x-,y-, and z-intercepts of the cut plane or
0009 %          a single value interpreted as the height of a single cut in the
0010 %          z plane (by default, a horizontal slice at the average electrode
0011 %          height will be created)
0012 %  OPT   = an option structure with the following fields and defaults:
0013 %     opt.imgsz = [32 32];    % dimensions of the pixel grid
0014 %     opt.square_pixels = 0;  % adjust imgsz to get square pixels
0015 %     opt.do_coarse2fine = 1; % calcuate c2f on the forward model
0016 %     opt.z_depth = inf;      % z_depth to use with mk_coarse_fine_mapping
0017 %
0018 % Output depends on the type of model suplied. If MDL is a fwd_model
0019 % structure then OUT is a rec_model. If MDL is an inv_model, then OUT is a
0020 % modified version of it, with the pixel slice in inv_model.rec_model and
0021 % updated inv_model.fwd_model.coarse2fine
0022 %
0023 % [OUT FMDL] = MK_PIXEL_SLICE(MDL, ...) also returns the forward model
0024 % structure with the coarse2fine field.
0025 %
0026 % See also MK_COARSE_FINE_MAPPING, MK_GRID_MODEL
0027 
0028 % (C) 2013 Bartlomiej Grychtol. License: GPL version 2 or 3
0029 % $Id: mk_pixel_slice.m 5112 2015-06-14 13:00:41Z aadler $
0030 
0031 if ischar(imdl) && strcmp(imdl,'UNIT_TEST'),do_unit_test;return;end;
0032 
0033 switch(imdl.type)
0034     case 'inv_model'
0035         fmdl = imdl.fwd_model;
0036     case 'fwd_model'
0037         fmdl = imdl;
0038     otherwise
0039         error('An EIDORS inverse or forward model struct required');
0040 end
0041 
0042 if nargin < 2, opt = struct; end
0043 if nargin > 1 
0044    if ~isstruct(level)
0045       opt.level = level;
0046    else
0047       opt = level;
0048    end
0049 end
0050 opt = parse_opts(fmdl,opt);
0051 
0052 [rmdl fmdl] = eidors_cache(@do_pixel_slice,{fmdl, opt},'mk_pixel_slice');
0053 
0054 switch imdl.type
0055    case 'inv_model'
0056       imdl.rec_model = rmdl;
0057       imdl.fwd_model = fmdl;
0058    case 'fwd_model'
0059       imdl = rmdl;
0060 end
0061 
0062 function [rmdl fmdl] = do_pixel_slice(fmdl, opt);
0063 tmp = fmdl;
0064 [NODES R T]=level_model(fmdl,opt.level);
0065 tmp.nodes = NODES';
0066 slc = mdl_slice_mesher(tmp,[inf inf 0]);
0067 slc = slc.fwd_model;
0068 mingrid = min(slc.nodes);
0069 maxgrid = max(slc.nodes);
0070 bnd = find_boundary(slc);
0071 % contour_boundary = order_loop(slc.nodes(unique(bnd),:));
0072 
0073 if opt.square_pixels ==1
0074     mdl_sz = maxgrid - mingrid;
0075     mdl_AR = mdl_sz(1)/mdl_sz(2);
0076     img_AR = opt.imgsz(1)/opt.imgsz(2);
0077     if mdl_AR < img_AR
0078         delta = (mdl_sz(2) * img_AR - mdl_sz(1)) /2;
0079         mingrid(1) = mingrid(1) - delta;
0080         maxgrid(1) = maxgrid(1) + delta;
0081     elseif mdl_AR > img_AR
0082         delta = (mdl_sz(1)/img_AR - mdl_sz(2)) / 2;
0083         mingrid(2) = mingrid(2) - delta;
0084         maxgrid(2) = maxgrid(2) + delta;
0085     end
0086 end
0087 
0088 xgrid = linspace(mingrid(1),maxgrid(1),opt.imgsz(1)+1);
0089 ygrid = linspace(mingrid(2),maxgrid(2),opt.imgsz(2)+1);
0090 rmdl = mk_grid_model([],xgrid,ygrid);
0091 x_pts = xgrid(1:end-1) + 0.5*diff(xgrid);
0092 y_pts = ygrid(1:end-1) + 0.5*diff(ygrid);
0093 % y_pts = fliplr(y_pts); %medical
0094 
0095 % NOTE: This controls the image resolution. If you want higher res, you
0096 % need to either specify it in opt.imgsz or manually overwrite (or remove)
0097 % the imdl.rec_model.mdl_slice_mapper.
0098 rmdl.mdl_slice_mapper.x_pts = x_pts;
0099 rmdl.mdl_slice_mapper.y_pts = y_pts;
0100 rmdl.mdl_slice_mapper.level = opt.level;
0101 rmdl.mdl_slice_mapper.model_2d = 1;
0102 x_avg = conv2(xgrid, [1,1]/2,'valid');
0103 y_avg = conv2(ygrid, [1,1]/2,'valid');
0104 [x,y] = ndgrid( x_avg, y_avg);
0105 
0106 % 20141119: The inpolygon approach fails on non-simply-connected domains
0107 % inside = inpolygon(x(:),y(:),contour_boundary(:,1),contour_boundary(:,2) );
0108 P = [x(:) y(:)]; % P(end,3) = 0;
0109 inside = any(point_in_triangle(P, slc.elems, slc.nodes(:,1:2)),2);
0110 % inside = full(inside);
0111 
0112 ff = find(~inside);
0113 
0114 if opt.do_coarse2fine
0115 %     rmdl.mk_coarse_fine_mapping.z_depth = opt.z_depth;
0116 %     fmdl.coarse2fine = mk_coarse_fine_mapping(tmp,rmdl);
0117 %     fmdl.coarse2fine(:,ff) = [];
0118     if isinf(opt.z_depth)
0119        zgrid = [min(tmp.nodes(:,3))-1 max(tmp.nodes(:,3))+1];
0120     else
0121        zgrid = [-opt.z_depth opt.z_depth];
0122     end
0123     [jnk, fmdl.coarse2fine] = mk_grid_model(tmp,xgrid,ygrid,zgrid);
0124     fmdl.coarse2fine(:,ff);
0125 end
0126 
0127 rmdl.elems([2*ff, 2*ff-1],:)= [];
0128 rmdl.coarse2fine([2*ff, 2*ff-1],:)= [];
0129 rmdl.coarse2fine(:,ff)= [];
0130 % rmdl.boundary = find_boundary(rmdl);
0131 % show individual elements (more like how the 2d grid models display)
0132 rmdl.boundary = rmdl.elems;
0133 rmdl.inside   = inside; % the inside array is useful in other functions
0134 
0135 
0136 if isfield(fmdl,'mat_idx')
0137    rmdl.mat_idx = calc_mat_idx(rmdl,fmdl,ff,opt);
0138 end
0139 
0140 % map electrodes
0141 rmdl.nodes(:,3) = 0;
0142 rmdl.nodes =  (R\rmdl.nodes' + T'*ones(1,length(rmdl.nodes)))';
0143 slc.nodes =  (R\slc.nodes' + T'*ones(1,length(slc.nodes)))';
0144 
0145 isf = ~isinf(opt.level);
0146 if nnz(isf) == 1
0147    rmdl.nodes(:,isf) = opt.level(:,isf);
0148 end
0149 
0150 if isfield(slc, 'electrode')
0151    for i = flipud(1:numel(slc.electrode))
0152         tmp = rmfield(slc.electrode(i), 'nodes');
0153         x_elec = slc.nodes( [slc.electrode(i).nodes], 1);
0154         y_elec = slc.nodes( [slc.electrode(i).nodes], 2);
0155         z_elec = slc.nodes( [slc.electrode(i).nodes], 3);
0156         tmp.pos       = [x_elec, y_elec, z_elec];
0157         elec(i) = tmp;
0158     end
0159     rmdl.electrode = elec;
0160 end
0161    
0162 
0163 rmdl.show_slices.levels = opt.level;
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,opt.level);
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')
0192         opt.level = get_elec_level(fmdl);
0193     else
0194         if 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 [NODE R T] = level_model( fwd_model, level )
0206 
0207    vtx= fwd_model.nodes;
0208    [nn, dims] = size(vtx);
0209    if dims ==2 % 2D case
0210        NODE= vtx';
0211        return;
0212    end
0213 
0214    % Infinities tend to cause issues -> replace with realmax
0215    % Don't need to worry about the sign of the inf
0216    level( isinf(level) | isnan(level) ) = realmax;
0217    level( level==0 ) =     1e-10; %eps;
0218 
0219    % Step 1: Choose a centre point in the plane
0220    %  Weight the point by it's inv axis coords
0221    invlev= 1./level;
0222    ctr= invlev / sum( invlev.^2 );
0223 
0224    % Step 2: Choose basis vectors in the plane
0225    %  First is the axis furthest from ctr
0226    [jnk, s_ax]= sort( - abs(level - ctr) );
0227    v1= [0,0,0]; v1(s_ax(1))= level(s_ax(1));
0228    v1= v1 - ctr;
0229    v1= v1 / norm(v1);
0230 
0231    % Step 3: Get off-plane vector, by cross product
0232    v2= [0,0,0]; v2(s_ax(2))= level(s_ax(2));
0233    v2= v2 - ctr;
0234    v2= v2 / norm(v2);
0235    v3= cross(v1,v2);
0236 
0237    % Step 4: Get orthonormal basis. Replace v2
0238    v2= cross(v1,v3);
0239 
0240    % Step 5: Get bases to point in 'positive directions'
0241    v1= v1 * (1-2*(sum(v1)<0));
0242    v2= v2 * (1-2*(sum(v2)<0));
0243    v3= v3 * (1-2*(sum(v3)<0));
0244    
0245    R = [v1;v2;v3];
0246    T = ctr;
0247 
0248    NODE= R * (vtx' - T'*ones(1,nn) );
0249 
0250 function elec_lev = get_elec_level(fmdl)
0251     z_elec= fmdl.nodes( [fmdl.electrode(:).nodes], 3);
0252     min_e = min(z_elec); max_e = max(z_elec);
0253     elec_lev = [inf,inf,mean([min_e,max_e])];
0254 
0255     
0256 function do_unit_test
0257     imdl = mk_common_model('n3r2',[16,2]);
0258     opt.square_pixels = 1;
0259     opt.imgsz = [16 16];
0260     mdl = mk_pixel_slice(imdl.fwd_model,[inf 2 2.5], opt);
0261     img = mk_image(mdl,1);
0262     
0263     subplot(231)
0264     show_fem(imdl.fwd_model);
0265     view([-50 10])
0266 
0267     subplot(232)
0268     show_fem(img);
0269     zlim([0 3]);
0270     ylim([-1 1])
0271     xlim([-1 1]);
0272     view([-50 10])
0273     
0274     subplot(233)
0275     show_slices(img);
0276     
0277     subplot(234)
0278     imdl = mk_pixel_slice(imdl);
0279     img = mk_image(imdl.rec_model,1);
0280     show_fem(img);
0281     zlim([0 3]);
0282     ylim([-1 1])
0283     xlim([-1 1]);
0284     view([-50 10])
0285     
0286     subplot(235)
0287     mdl = mk_library_model('pig_23kg_16el_lungs');
0288     img = mk_image(mdl,1);
0289     for i = 1:length(mdl.mat_idx)
0290        img.elem_data(mdl.mat_idx{i}) = i;
0291     end
0292     show_fem(img)
0293     view(2)
0294     
0295     subplot(236)
0296     clear opt
0297     opt.imgsz = [64 64];
0298     opt.square_pixels = 1;
0299     opt.do_coarse2fine = 0;
0300     mdl = mk_library_model('pig_23kg_16el_lungs');
0301     rmdl = mk_pixel_slice(mdl,opt);
0302     img = mk_image(rmdl,1);
0303     for i = 1:length(rmdl.mat_idx)
0304        img.elem_data(rmdl.mat_idx{i}) = i;
0305     end
0306     show_slices(img);