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 6392 2022-08-30 13:54:01Z bgrychtol $
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 % Old implementation, replaced by mk_grid_c2f
0091 function c2f= calc_c2f_3d( fmdl, xvec, yvec, zvec);
0092 %  c2f= mk_coarse_fine_mapping( fmdl, cmdl);
0093    nef= size( fmdl.elems,1);
0094 %  c2f= sparse(nef,0);
0095    c2fiidx= [];
0096    c2fjidx= [];
0097    c2fdata= [];
0098    jidx= 0;
0099    mdl_pts = interp_mesh( fmdl, 3);
0100    x_pts = squeeze(mdl_pts(:,1,:));
0101    y_pts = squeeze(mdl_pts(:,2,:));
0102    z_pts = squeeze(mdl_pts(:,3,:));
0103    
0104    in_x_pts = calc_in_d_pts( x_pts, xvec);
0105    in_y_pts = calc_in_d_pts( y_pts, yvec);
0106    in_z_pts = calc_in_d_pts( z_pts, zvec);
0107 
0108    for zi= 1:length(zvec)-1
0109       for yi= 1:length(yvec)-1
0110              in_yz_pts = in_y_pts{yi} & in_z_pts{zi};
0111          for xi= 1:length(xvec)-1
0112              in_pts = mean( in_x_pts{xi} & in_yz_pts, 2);
0113              % c2f = [c2f,sparse(in_pts)];
0114              [ii,jj,vv] = find(in_pts);
0115              c2fiidx= [c2fiidx;ii];
0116              c2fjidx= [c2fjidx;jj+jidx]; jidx=jidx+1;
0117              c2fdata= [c2fdata;vv];
0118          end
0119       end
0120    end
0121    c2f= sparse(c2fiidx,c2fjidx,c2fdata, length(in_pts), jidx);
0122 
0123 function cmdl= mk_2d_grid(xvec, yvec);
0124    xlen = length(xvec);
0125    ylen = length(yvec);
0126    cmdl= eidors_obj('fwd_model', ...
0127             sprintf('Grid model %d x %d', xlen, ylen) );
0128 
0129    [x,y]= ndgrid( xvec, yvec);
0130    cmdl.nodes= [x(:),y(:)];
0131    k= 1:xlen-1;
0132    elem_frac = [ k;k+1;k+xlen; ...
0133                  k+1;k+xlen;k+xlen+1];
0134    elem_frac= reshape(elem_frac, 3,[])';
0135    cmdl.elems=  [];
0136    for j=0:ylen-2
0137       cmdl.elems=  [cmdl.elems; elem_frac + xlen*j];
0138    end
0139 
0140    cmdl.boundary = find_boundary( cmdl.elems);
0141 
0142 % assign one single parameter to each square element
0143    e= size(cmdl.elems,1);
0144    params= ceil(( 1:e )/2);
0145    cmdl.coarse2fine = sparse(1:e,params,1,e,max(params));
0146 
0147 
0148 function cmdl= mk_3d_grid(xvec, yvec, zvec)
0149    xlen = length(xvec);
0150    ylen = length(yvec);
0151    zlen = length(zvec);
0152    if zlen<2 || ylen<2 || xlen<2
0153       error('Need at least 2 components for each gridpoint')
0154    end
0155    cmdl= eidors_obj('fwd_model', ...
0156             sprintf('Grid model %d x %d x %d', xlen, ylen, zlen) );
0157 
0158    [x,y,z]= ndgrid( xvec, yvec, zvec);
0159    cmdl.nodes= [x(:),y(:),z(:)];
0160    k= 1:xlen-1;
0161    ac = xlen; up = xlen*ylen; % accross vs up
0162    elem_frac = [ k;     k+1 ;  k+ac;   k+up;  ...
0163                  k+1;   k+ac;  k+up;   k+up+1; ...
0164                  k+ac;  k+up;  k+up+1; k+up+ac; ...
0165                  k+1;   k+ac;  k+ac+1; k+up+1; ...
0166                  k+ac;  k+ac+1;k+up+1; k+up+ac; ...
0167                  k+ac+1;k+up+1;k+up+ac;k+up+ac+1];
0168    elem_frac= reshape(elem_frac, 4,[])';
0169    sz_elem_frac = size(elem_frac);
0170    row_frac =  zeros(sz_elem_frac .* [ylen-1,1]);
0171    for j=0:ylen-2
0172       idx = (1:sz_elem_frac(1)) + j*sz_elem_frac(1);
0173       row_frac(idx,:)=  elem_frac + ac*j;
0174    end
0175    
0176    sz_row_frac = size(row_frac);
0177    cmdl.elems=  zeros(sz_row_frac .* [zlen-1,1]);
0178    for k=0:zlen-2
0179       idx = (1:sz_row_frac(1)) + k*sz_row_frac(1);
0180       cmdl.elems(idx,:) =  row_frac + up*k;
0181    end
0182 
0183    cmdl.boundary = find_boundary( cmdl.elems);
0184 
0185 % assign one single parameter to each square element
0186    e= size(cmdl.elems,1);
0187    params= ceil(( 1:e )/6);
0188    cmdl.coarse2fine = sparse(1:e,params,1,e,max(params));
0189 
0190 function mdl = set_pixel_pos(mdl, xvec, yvec)
0191    x = xvec(1:end-1) + 0.5*diff(xvec);
0192    y = yvec(1:end-1) + 0.5*diff(yvec);
0193    y = y(end:-1:1); %get the medical orientation right
0194    mdl.mdl_slice_mapper.x_pts = x;
0195    mdl.mdl_slice_mapper.y_pts = y;
0196    
0197    
0198 function in_d_pts = calc_in_d_pts( d_pts, dvec);
0199    l1dvec= length(dvec)-1;
0200    in_d_pts = cell(l1dvec,1);
0201    for i= 1:l1dvec
0202       in_d_pts{i} = d_pts >= dvec(i) & d_pts < dvec(i+1);
0203    end
0204 
0205 function do_unit_test
0206 imdl = mk_common_model('b2c2',16); imdl.hyperparameter.value = 1e-3;
0207 img = mk_image(imdl,1);     vh= fwd_solve(img);
0208 img.elem_data([51,23])=1.1; vi= fwd_solve(img);
0209 subplot(321); show_fem(img);
0210 subplot(322); show_fem(inv_solve(imdl, vh, vi));
0211 
0212 grid = linspace(-1,1,33);
0213 [imdl.rec_model, imdl.fwd_model.coarse2fine] = ...
0214      mk_grid_model( imdl.fwd_model, grid, grid );
0215 subplot(323); show_fem(inv_solve(imdl, vh, vi));
0216 hold on; hh=show_fem(img); set(hh,'FaceAlpha',0,'EdgeColor',[0,0,1]); hold off;
0217 
0218 outside = find(sum(imdl.fwd_model.coarse2fine,1) < eps);
0219 imdl.fwd_model.coarse2fine(:,outside) = [];
0220 imdl.rec_model.coarse2fine(:,outside) = [];
0221 rec_out = [2*outside-1,2*outside];
0222 imdl.rec_model.coarse2fine(rec_out,:) = [];
0223 imdl.rec_model.elems(rec_out,:) = [];
0224 subplot(324); show_fem(inv_solve(imdl, vh, vi));
0225 hold on; hh=show_fem(img); set(hh,'FaceAlpha',0,'EdgeColor',[0,0,1]); hold off;
0226 
0227 subplot(325)
0228 imdl = mk_common_model('n3r2',[16,2]);
0229 img = mk_image(imdl.fwd_model,1:size(imdl.fwd_model.elems,1));
0230 show_fem(img);
0231 
0232 subplot(326)
0233 [cmdl, c2f]= mk_grid_model(img.fwd_model, -.5:.1:.5, -.5:.1:.5, 1:.1:2);
0234 v = c2f'*img.elem_data;
0235 img2 = mk_image(cmdl, v);
0236 show_fem(img2)
0237