mk_grid_model

 MK_GRID_MODEL: Create reconstruction model on pixelated grid 
  [cmdl,coarse2fine]= mk_grid_model(fmdl, xvec, yvec, zvec);

 Outputs:
  cmdl - eidors reconstruction model (coarse model)
  coarse2fine - c2f mapping to put onto fmdl (specify [] to not use)

 Inputs:
  fmdl - fine model (forward model) to create coarse2fine mapping
  xvec - x edges
  yvec - y edges
  zvec - z edges (optional - to create 3D model)

 if fmdl == [], then just create the grid model without c2f

 See also MK_COARSE_FINE_MAPPING, MK_PIXEL_SLICE

 Example: for constructing an inverse model
  grid{1}= linspace(-2,2,20);     % x grid
  grid{2}= linspace(-0.5,+0.5,5); % y grid
  grid{3}= linspace(-2, 0,20);    % z grid
  imdl = select_imdl( fmdl, {'Basic GN dif'});
  [imdl.rec_model,imdl.fwd_model.coarse2fine]= mk_grid_model(fmdl,grid{:});

 ISSUES:
  Ensure that grids are defined from smallest to largest

0001 function [cmdl, c2f]= mk_grid_model(fmdl, xvec, yvec, zvec);
0002 % MK_GRID_MODEL: Create reconstruction model on pixelated grid
0003 %  [cmdl,coarse2fine]= mk_grid_model(fmdl, xvec, yvec, zvec);
0004 %
0005 % Outputs:
0006 %  cmdl - eidors reconstruction model (coarse model)
0007 %  coarse2fine - c2f mapping to put onto fmdl (specify [] to not use)
0008 %
0009 % Inputs:
0010 %  fmdl - fine model (forward model) to create coarse2fine mapping
0011 %  xvec - x edges
0012 %  yvec - y edges
0013 %  zvec - z edges (optional - to create 3D model)
0014 %
0015 % if fmdl == [], then just create the grid model without c2f
0016 %
0017 % See also MK_COARSE_FINE_MAPPING, MK_PIXEL_SLICE
0018 %
0019 % Example: for constructing an inverse model
0020 %  grid{1}= linspace(-2,2,20);     % x grid
0021 %  grid{2}= linspace(-0.5,+0.5,5); % y grid
0022 %  grid{3}= linspace(-2, 0,20);    % z grid
0023 %  imdl = select_imdl( fmdl, {'Basic GN dif'});
0024 %  [imdl.rec_model,imdl.fwd_model.coarse2fine]= mk_grid_model(fmdl,grid{:});
0025 %
0026 % ISSUES:
0027 %  Ensure that grids are defined from smallest to largest
0028 
0029 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0030 % $Id: mk_grid_model.m 5902 2019-02-04 10:18:28Z aadler $
0031 
0032 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0033 
0034 
0035 if nargin == 3
0036    cmdl = mk_2d_grid(xvec,yvec);
0037 elseif nargin ==4
0038    cmdl = mk_3d_grid(xvec,yvec,zvec);
0039 else
0040    error('check nargin');
0041 end
0042 
0043 % this had too many side effects
0044 cmdl = set_pixel_pos(cmdl,xvec,yvec);% same for 2d and 3d
0045 
0046 % put in the centre (or near it)
0047 ctr = ones(num_nodes(cmdl),1)*mean(cmdl.nodes);
0048 dctr= sum( (cmdl.nodes - ctr).^2, 2);
0049 [jnk, c_idx] = min(dctr);
0050 cmdl.gnd_node = c_idx(1);
0051 
0052 if ~isempty( fmdl)
0053     if size(fmdl.nodes,2) == 2
0054         assert(nargin==3);
0055         c2f= calc_c2f_2d( fmdl, xvec, yvec);
0056         
0057     else
0058         if nargin == 3
0059             % here we could incorporate z_depth
0060             zvec = [ min(fmdl.nodes(:,3)) - 1; max(fmdl.nodes(:,3))+1 ];
0061             tmp = mk_3d_grid(xvec,yvec,zvec);
0062         elseif nargin == 4
0063             tmp = cmdl;
0064         end
0065         c2f = mk_grid_c2f(fmdl,tmp);
0066     end
0067 end
0068 
0069 cmdl.normalize_measurements = 0;% @eidors_default;
0070 cmdl.solve =      @eidors_default;
0071 cmdl.system_mat = @eidors_default;
0072 cmdl.jacobian   = @eidors_default;
0073 
0074 
0075 function c2f= calc_c2f_2d( fmdl, xvec, yvec);
0076    nef= size( fmdl.elems,1);
0077    c2f= sparse(nef,0);
0078    mdl_pts = interp_mesh( fmdl, 3);
0079    x_pts = squeeze(mdl_pts(:,1,:));
0080    y_pts = squeeze(mdl_pts(:,2,:));
0081    for yi= 1:length(yvec)-1
0082          in_y_pts = y_pts >= yvec(yi) & y_pts < yvec(yi+1);
0083       for xi= 1:length(xvec)-1
0084           in_x_pts =  x_pts >= xvec(xi) & x_pts < xvec(xi+1);
0085           in_pts = mean( in_y_pts & in_x_pts , 2);
0086           c2f = [c2f,sparse(in_pts)];
0087       end
0088    end
0089 
0090 function c2f= calc_c2f_3d( fmdl, xvec, yvec, zvec);
0091 %  c2f= mk_coarse_fine_mapping( fmdl, cmdl);
0092    nef= size( fmdl.elems,1);
0093 %  c2f= sparse(nef,0);
0094    c2fiidx= [];
0095    c2fjidx= [];
0096    c2fdata= [];
0097    jidx= 0;
0098    mdl_pts = interp_mesh( fmdl, 3);
0099    x_pts = squeeze(mdl_pts(:,1,:));
0100    y_pts = squeeze(mdl_pts(:,2,:));
0101    z_pts = squeeze(mdl_pts(:,3,:));
0102    
0103    in_x_pts = calc_in_d_pts( x_pts, xvec);
0104    in_y_pts = calc_in_d_pts( y_pts, yvec);
0105    in_z_pts = calc_in_d_pts( z_pts, zvec);
0106 
0107    for zi= 1:length(zvec)-1
0108       for yi= 1:length(yvec)-1
0109              in_yz_pts = in_y_pts{yi} & in_z_pts{zi};
0110          for xi= 1:length(xvec)-1
0111              in_pts = mean( in_x_pts{xi} & in_yz_pts, 2);
0112              % c2f = [c2f,sparse(in_pts)];
0113              [ii,jj,vv] = find(in_pts);
0114              c2fiidx= [c2fiidx;ii];
0115              c2fjidx= [c2fjidx;jj+jidx]; jidx=jidx+1;
0116              c2fdata= [c2fdata;vv];
0117          end
0118       end
0119    end
0120    c2f= sparse(c2fiidx,c2fjidx,c2fdata, length(in_pts), jidx);
0121 
0122 function cmdl= mk_2d_grid(xvec, yvec);
0123    xlen = length(xvec);
0124    ylen = length(yvec);
0125    cmdl= eidors_obj('fwd_model', ...
0126             sprintf('Grid model %d x %d', xlen, ylen) );
0127 
0128    [x,y]= ndgrid( xvec, yvec);
0129    cmdl.nodes= [x(:),y(:)];
0130    k= 1:xlen-1;
0131    elem_frac = [ k;k+1;k+xlen; ...
0132                  k+1;k+xlen;k+xlen+1];
0133    elem_frac= reshape(elem_frac, 3,[])';
0134    cmdl.elems=  [];
0135    for j=0:ylen-2
0136       cmdl.elems=  [cmdl.elems; elem_frac + xlen*j];
0137    end
0138 
0139    cmdl.boundary = find_boundary( cmdl.elems);
0140 
0141 % assign one single parameter to each square element
0142    e= size(cmdl.elems,1);
0143    params= ceil(( 1:e )/2);
0144    cmdl.coarse2fine = sparse(1:e,params,1,e,max(params));
0145 
0146 
0147 function cmdl= mk_3d_grid(xvec, yvec, zvec);
0148    xlen = length(xvec);
0149    ylen = length(yvec);
0150    zlen = length(zvec);
0151    cmdl= eidors_obj('fwd_model', ...
0152             sprintf('Grid model %d x %d x %d', xlen, ylen, zlen) );
0153 
0154    [x,y,z]= ndgrid( xvec, yvec, zvec);
0155    cmdl.nodes= [x(:),y(:),z(:)];
0156    k= 1:xlen-1;
0157    ac = xlen; up = xlen*ylen; % accross vs up
0158    elem_frac = [ k;     k+1 ;  k+ac;   k+up;  ...
0159                  k+1;   k+ac;  k+up;   k+up+1; ...
0160                  k+ac;  k+up;  k+up+1; k+up+ac; ...
0161                  k+1;   k+ac;  k+ac+1; k+up+1; ...
0162                  k+ac;  k+ac+1;k+up+1; k+up+ac; ...
0163                  k+ac+1;k+up+1;k+up+ac;k+up+ac+1];
0164    elem_frac= reshape(elem_frac, 4,[])';
0165 
0166    row_frac =  [];
0167    for j=0:ylen-2
0168       row_frac=  [row_frac; elem_frac + ac*j];
0169    end
0170 
0171    cmdl.elems=  [];
0172    for k=0:zlen-2
0173       cmdl.elems=  [cmdl.elems; row_frac + up*k];
0174    end
0175 
0176    cmdl.boundary = find_boundary( cmdl.elems);
0177 
0178 % assign one single parameter to each square element
0179    e= size(cmdl.elems,1);
0180    params= ceil(( 1:e )/6);
0181    cmdl.coarse2fine = sparse(1:e,params,1,e,max(params));
0182 
0183 function mdl = set_pixel_pos(mdl, xvec, yvec)
0184    x = xvec(1:end-1) + 0.5*diff(xvec);
0185    y = yvec(1:end-1) + 0.5*diff(yvec);
0186    y = y(end:-1:1); %get the medical orientation right
0187    mdl.mdl_slice_mapper.x_pts = x;
0188    mdl.mdl_slice_mapper.y_pts = y;
0189    
0190    
0191 function in_d_pts = calc_in_d_pts( d_pts, dvec);
0192    l1dvec= length(dvec)-1;
0193    in_d_pts = cell(l1dvec,1);
0194    for i= 1:l1dvec
0195       in_d_pts{i} = d_pts >= dvec(i) & d_pts < dvec(i+1);
0196    end
0197 
0198 function do_unit_test
0199 imdl = mk_common_model('b2c2',16); imdl.hyperparameter.value = 1e-3;
0200 img = mk_image(imdl,1);     vh= fwd_solve(img);
0201 img.elem_data([51,23])=1.1; vi= fwd_solve(img);
0202 subplot(221); show_fem(img);
0203 subplot(222); show_fem(inv_solve(imdl, vh, vi));
0204 
0205 grid = linspace(-1,1,33);
0206 [imdl.rec_model, imdl.fwd_model.coarse2fine] = ...
0207      mk_grid_model( imdl.fwd_model, grid, grid );
0208 subplot(223); show_fem(inv_solve(imdl, vh, vi));
0209 hold on; hh=show_fem(img); set(hh,'FaceAlpha',0,'EdgeColor',[0,0,1]); hold off;
0210 
0211 outside = find(sum(imdl.fwd_model.coarse2fine,1) < eps);
0212 imdl.fwd_model.coarse2fine(:,outside) = [];
0213 imdl.rec_model.coarse2fine(:,outside) = [];
0214 rec_out = [2*outside-1,2*outside];
0215 imdl.rec_model.coarse2fine(rec_out,:) = [];
0216 imdl.rec_model.elems(rec_out,:) = [];
0217 subplot(224); show_fem(inv_solve(imdl, vh, vi));
0218 hold on; hh=show_fem(img); set(hh,'FaceAlpha',0,'EdgeColor',[0,0,1]); hold off;