mk_grid_model

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

 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)
  removefn(xyzc) - function which elements to remove, given xyzctr
  removefn = 'outside' => remove elements outside the fmdl
    e.g. removefn = @(xyz) vecnorm(xyz(:,1:2),2,2)>1.0;

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

 Example: for circular model
  grid= linspace(-2,2,20);     % x,y grid
  [gmdl]= mk_grid_model([],grid,grid,[], @(xyz) vecnorm(xyz,2,2)>1);

 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{:});

 See also MK_COARSE_FINE_MAPPING, MK_PIXEL_SLICE

0001 function [cmdl, c2f]= mk_grid_model(fmdl, xvec, yvec, zvec, removefn);
0002 % MK_GRID_MODEL: Create reconstruction model on pixelated grid
0003 %  [cmdl,coarse2fine]= mk_grid_model(fmdl, xvec, yvec, zvec, removefn);
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 %  removefn(xyzc) - function which elements to remove, given xyzctr
0015 %  removefn = 'outside' => remove elements outside the fmdl
0016 %    e.g. removefn = @(xyz) vecnorm(xyz(:,1:2),2,2)>1.0;
0017 %
0018 % if fmdl == [], then just create the grid model without c2f
0019 %
0020 % Example: for circular model
0021 %  grid= linspace(-2,2,20);     % x,y grid
0022 %  [gmdl]= mk_grid_model([],grid,grid,[], @(xyz) vecnorm(xyz,2,2)>1);
0023 %
0024 % Example: for constructing an inverse model
0025 %  grid{1}= linspace(-2,2,20);     % x grid
0026 %  grid{2}= linspace(-0.5,+0.5,5); % y grid
0027 %  grid{3}= linspace(-2, 0,20);    % z grid
0028 %  imdl = select_imdl( fmdl, {'Basic GN dif'});
0029 %  [imdl.rec_model,imdl.fwd_model.coarse2fine]= mk_grid_model(fmdl,grid{:});
0030 %
0031 % See also MK_COARSE_FINE_MAPPING, MK_PIXEL_SLICE
0032 
0033 % ISSUES:
0034 %  Ensure that grids are defined from smallest to largest
0035 
0036 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0037 % $Id: mk_grid_model.m 7044 2024-11-30 15:53:22Z aadler $
0038 
0039 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0040 
0041 
0042 if nargin<5 % removefn not provided
0043    removefn = [];
0044 end
0045 if nargin == 3 || isempty(zvec)
0046    do3d= false;
0047 elseif nargin > 3
0048    do3d=true;
0049 else
0050    error('check nargin');
0051 end
0052 
0053 if ~do3d
0054    [cmdl,rmelems] = mk_2d_grid(xvec,yvec,removefn);
0055 else
0056    [cmdl,rmelems] = mk_3d_grid(xvec,yvec,zvec,removefn);
0057 end
0058 
0059 % this had too many side effects
0060 cmdl = set_pixel_pos(cmdl,xvec,yvec);% same for 2d and 3d
0061 
0062 % put in the centre (or near it)
0063 ctr = ones(num_nodes(cmdl),1)*mean(cmdl.nodes);
0064 dctr= sum( (cmdl.nodes - ctr).^2, 2);
0065 [jnk, c_idx] = min(dctr);
0066 cmdl.gnd_node = c_idx(1);
0067 
0068 if ~isempty( fmdl)
0069     if size(fmdl.nodes,2) == 2
0070         assert(~do3d);
0071         c2f= calc_c2f_2d( fmdl, xvec, yvec);
0072         
0073     else
0074         if ~do3d
0075             % here we could incorporate z_depth
0076             zvec = [ min(fmdl.nodes(:,3)) - 1; max(fmdl.nodes(:,3))+1 ];
0077             tmp = mk_3d_grid(xvec,yvec,zvec,removefn);
0078         else
0079             tmp = cmdl;
0080         end
0081         c2f = mk_grid_c2f(fmdl,tmp);
0082     end
0083 else
0084     c2f = [];
0085 end
0086 
0087 % now remove elements outsize shape
0088 if ischar(rmelems{2}) && strcmp(rmelems{2},'OUTSIDE')
0089    rmelems{2} = (sum(c2f,1) < eps);
0090    rmelems{1} = logical(kron(rmelems{2}, true(1,rmelems{1})));
0091 end
0092 if ~isempty(rmelems{1})
0093    cmdl.elems( rmelems{1}, : ) = [];
0094    cmdl.coarse2fine( rmelems{1},     :      ) = [];
0095    cmdl.coarse2fine(    :      , rmelems{2} ) = [];
0096    if ~isempty(c2f)
0097       c2f(:,rmelems{2}) = [];
0098    end
0099 %  cmdl.remover = rmelems; % debugging
0100 end
0101 
0102 cmdl.boundary = find_boundary( cmdl.elems);
0103 
0104 cmdl = mdl_normalize(cmdl, 'default');
0105 cmdl.solve =      @eidors_default;
0106 cmdl.system_mat = @eidors_default;
0107 cmdl.jacobian   = @eidors_default;
0108 
0109 % standard field order
0110 cmdl = eidors_obj('set', cmdl);
0111 
0112 function c2f= calc_c2f_2d( fmdl, xvec, yvec);
0113    nef= size( fmdl.elems,1);
0114    c2f= sparse(nef,0);
0115    mdl_pts = interp_mesh( fmdl, 3);
0116    x_pts = squeeze(mdl_pts(:,1,:));
0117    y_pts = squeeze(mdl_pts(:,2,:));
0118    for yi= 1:length(yvec)-1
0119          in_y_pts = y_pts >= yvec(yi) & y_pts < yvec(yi+1);
0120       for xi= 1:length(xvec)-1
0121           in_x_pts =  x_pts >= xvec(xi) & x_pts < xvec(xi+1);
0122           in_pts = mean( in_y_pts & in_x_pts , 2);
0123           c2f = [c2f,sparse(in_pts)];
0124       end
0125    end
0126 
0127 % Old implementation, replaced by mk_grid_c2f
0128 function c2f= calc_c2f_3d( fmdl, xvec, yvec, zvec);
0129 %  c2f= mk_coarse_fine_mapping( fmdl, cmdl);
0130    nef= size( fmdl.elems,1);
0131 %  c2f= sparse(nef,0);
0132    c2fiidx= [];
0133    c2fjidx= [];
0134    c2fdata= [];
0135    jidx= 0;
0136    mdl_pts = interp_mesh( fmdl, 3);
0137    x_pts = squeeze(mdl_pts(:,1,:));
0138    y_pts = squeeze(mdl_pts(:,2,:));
0139    z_pts = squeeze(mdl_pts(:,3,:));
0140    
0141    in_x_pts = calc_in_d_pts( x_pts, xvec);
0142    in_y_pts = calc_in_d_pts( y_pts, yvec);
0143    in_z_pts = calc_in_d_pts( z_pts, zvec);
0144 
0145    for zi= 1:length(zvec)-1
0146       for yi= 1:length(yvec)-1
0147              in_yz_pts = in_y_pts{yi} & in_z_pts{zi};
0148          for xi= 1:length(xvec)-1
0149              in_pts = mean( in_x_pts{xi} & in_yz_pts, 2);
0150              % c2f = [c2f,sparse(in_pts)];
0151              [ii,jj,vv] = find(in_pts);
0152              c2fiidx= [c2fiidx;ii];
0153              c2fjidx= [c2fjidx;jj+jidx]; jidx=jidx+1;
0154              c2fdata= [c2fdata;vv];
0155          end
0156       end
0157    end
0158    c2f= sparse(c2fiidx,c2fjidx,c2fdata, length(in_pts), jidx);
0159 
0160 function [cmdl,rmelems]= mk_2d_grid(xvec, yvec, removefn);
0161    xlen = length(xvec);
0162    ylen = length(yvec);
0163    cmdl= eidors_obj('fwd_model', ...
0164             sprintf('Grid model %d x %d', xlen, ylen) );
0165 
0166    [x,y]= ndgrid( xvec, yvec);
0167    cmdl.nodes= [x(:),y(:)];
0168    k= 1:xlen-1;
0169    elem_frac = [ k;k+1;k+xlen; ...
0170                  k+1;k+xlen;k+xlen+1];
0171    elem_frac= reshape(elem_frac, 3,[])';
0172    cmdl.elems=  [];
0173    for j=0:ylen-2
0174       cmdl.elems=  [cmdl.elems; elem_frac + xlen*j];
0175    end
0176 
0177    if isempty(removefn)
0178       rmelems{1} = [];
0179       rmelems{2} = [];
0180    elseif ischar(removefn) && strcmp(upper(removefn),'OUTSIDE') 
0181       rmelems{2} = 'OUTSIDE'; % solve later
0182       rmelems{1} = 2; % divide factor
0183    else % it's a function
0184       % Ensure we get centre of each grid square from triangles
0185       xyzc = interp_mesh(cmdl);
0186       xyzc = kron(xyzc(1:2:end,:) + xyzc(2:2:end,:),[1;1]/2);
0187       rmelems{1} = removefn(xyzc);
0188       rmelems{2} = rmelems{1}(1:2:end);
0189    end
0190 
0191 % assign one single parameter to each square element
0192    e= size(cmdl.elems,1);
0193    params= ceil(( 1:e )/2);
0194    cmdl.coarse2fine = sparse(1:e,params,1,e,max(params));
0195 
0196    cmdl.boundary = find_boundary( cmdl.elems);
0197 
0198 function [cmdl,rmelems]= mk_3d_grid(xvec, yvec, zvec,removefn)
0199    xlen = length(xvec);
0200    ylen = length(yvec);
0201    zlen = length(zvec);
0202    if zlen<2 || ylen<2 || xlen<2
0203       error('Need at least 2 components for each gridpoint')
0204    end
0205    cmdl= eidors_obj('fwd_model', ...
0206             sprintf('Grid model %d x %d x %d', xlen, ylen, zlen) );
0207 
0208    [x,y,z]= ndgrid( xvec, yvec, zvec);
0209    cmdl.nodes= [x(:),y(:),z(:)];
0210    k= 1:xlen-1;
0211    ac = xlen; up = xlen*ylen; % accross vs up
0212    elem_frac = [ k;     k+1 ;  k+ac;   k+up;  ...
0213                  k+1;   k+ac;  k+up;   k+up+1; ...
0214                  k+ac;  k+up;  k+up+1; k+up+ac; ...
0215                  k+1;   k+ac;  k+ac+1; k+up+1; ...
0216                  k+ac;  k+ac+1;k+up+1; k+up+ac; ...
0217                  k+ac+1;k+up+1;k+up+ac;k+up+ac+1];
0218    elem_frac= reshape(elem_frac, 4,[])';
0219    sz_elem_frac = size(elem_frac);
0220    row_frac =  zeros(sz_elem_frac .* [ylen-1,1]);
0221    for j=0:ylen-2
0222       idx = (1:sz_elem_frac(1)) + j*sz_elem_frac(1);
0223       row_frac(idx,:)=  elem_frac + ac*j;
0224    end
0225    
0226    sz_row_frac = size(row_frac);
0227    cmdl.elems=  zeros(sz_row_frac .* [zlen-1,1]);
0228    for k=0:zlen-2
0229       idx = (1:sz_row_frac(1)) + k*sz_row_frac(1);
0230       cmdl.elems(idx,:) =  row_frac + up*k;
0231    end
0232 
0233    if isempty(removefn)
0234       rmelems{1} = [];
0235       rmelems{2} = [];
0236    elseif ischar(removefn) && strcmp(upper(removefn),'OUTSIDE') 
0237       rmelems{2} = 'OUTSIDE'; % solve later
0238       rmelems{1} = 6; % divide factor
0239    else % it's a function
0240       % Ensure we get centre of each grid square from triangles
0241       xyzc = interp_mesh(cmdl);
0242       xyzc = kron(xyzc(1:6:end,:) + ...
0243                   xyzc(2:6:end,:) + ...
0244                   xyzc(3:6:end,:) + ...
0245                   xyzc(4:6:end,:) + ...
0246                   xyzc(5:6:end,:) + ...
0247                   xyzc(6:6:end,:),[1;1;1;1;1;1]/6);
0248       rmelems{1} = removefn(xyzc);
0249       rmelems{2} = rmelems{1}(1:6:end);
0250    end
0251 
0252 % assign one single parameter to each square element
0253    e= size(cmdl.elems,1);
0254    params= ceil(( 1:e )/6);
0255    cmdl.coarse2fine = sparse(1:e,params,1,e,max(params));
0256 
0257 function mdl = set_pixel_pos(mdl, xvec, yvec)
0258    x = xvec(1:end-1) + 0.5*diff(xvec);
0259    y = yvec(1:end-1) + 0.5*diff(yvec);
0260    y = y(end:-1:1); %get the medical orientation right
0261    mdl.mdl_slice_mapper.x_pts = x;
0262    mdl.mdl_slice_mapper.y_pts = y;
0263    
0264    
0265 function in_d_pts = calc_in_d_pts( d_pts, dvec);
0266    l1dvec= length(dvec)-1;
0267    in_d_pts = cell(l1dvec,1);
0268    for i= 1:l1dvec
0269       in_d_pts{i} = d_pts >= dvec(i) & d_pts < dvec(i+1);
0270    end
0271 
0272 function do_unit_test
0273 imdl = mk_common_model('b2c2',16); imdl.hyperparameter.value = 1e-3;
0274 img = mk_image(imdl,1);     vh= fwd_solve(img);
0275 img.elem_data([51,23])=1.1; vi= fwd_solve(img);
0276 subplot(3,4,1); show_fem(img);
0277 subplot(3,4,2); show_fem(inv_solve(imdl, vh, vi));
0278 
0279 grid = linspace(-1,1,33);
0280 [imdl.rec_model, imdl.fwd_model.coarse2fine] = ...
0281      mk_grid_model( imdl.fwd_model, grid, grid );
0282 subplot(3,4,3); show_fem(inv_solve(imdl, vh, vi));
0283 hold on; hh=show_fem(img); set(hh,'FaceAlpha',0,'EdgeColor',[0,0,1]); hold off;
0284 
0285 if 0
0286   outside = find(sum(imdl.fwd_model.coarse2fine,1) < eps);
0287   imdl.fwd_model.coarse2fine(:,outside) = [];
0288   imdl.rec_model.coarse2fine(:,outside) = [];
0289   rec_out = [2*outside-1,2*outside];
0290   imdl.rec_model.coarse2fine(rec_out,:) = [];
0291   imdl.rec_model.elems(rec_out,:) = [];
0292 else
0293   rmfn = 'outside';
0294   [imdl.rec_model, imdl.fwd_model.coarse2fine] = ...
0295        mk_grid_model( imdl.fwd_model, grid, grid, [],rmfn );
0296 end
0297 subplot(3,4,4); show_fem(inv_solve(imdl, vh, vi));
0298 hold on; hh=show_fem(img); set(hh,'FaceAlpha',0,'EdgeColor',[0,0,1]); hold off;
0299 
0300 rmfn = @(xyz) vecnorm(xyz,2,2)>1;
0301 [imdl.rec_model, imdl.fwd_model.coarse2fine] = ...
0302      mk_grid_model( imdl.fwd_model, grid, grid, [],rmfn );
0303 subplot(3,4,5); show_fem(inv_solve(imdl, vh, vi));
0304 hold on; hh=show_fem(img); set(hh,'FaceAlpha',0,'EdgeColor',[0,0,1]); hold off;
0305 
0306 subplot(3,4,6)
0307 imdl = mk_common_model('n3r2',[16,2]);
0308 A= [390 391 393 396 402 478 479 480 484 486 664 665 666 667 668 670 671 672 676 677 678 755 760 761];
0309 B= [318 319 321 324 330 439 440 441 445 447 592 593 594 595 596 598 599 600 604 605 606 716 721 722];
0310 img = mk_image(imdl.fwd_model,1); vh = fwd_solve(img);
0311 img.elem_data(A)=1.5; img.elem_data(B)=0.5; vi=fwd_solve(img);
0312 show_fem(img);
0313 
0314 subplot(3,4,7)
0315 [cmdl, c2f]= mk_grid_model(img.fwd_model, -.8:.1:.8, -.8:.1:.8, 0:.5:3);
0316 v = c2f'*(img.elem_data-1) + 1;
0317 img2 = mk_image(cmdl, v);
0318 show_fem(img2)
0319 
0320 subplot(3,4,8)
0321 rmfn = @(xyz) vecnorm(xyz(:,1:2),2,2)>1.0;
0322 [imdl.rec_model, imdl.fwd_model.coarse2fine] = ...
0323     mk_grid_model(img.fwd_model, -1:.1:1, -1:.1:1, 0:.5:3,rmfn );
0324 %v = imdl.fwd_model.coarse2fine'*(img.elem_data-1)+1;
0325 %img2 = mk_image(imdl.rec_model, v);
0326 img2= inv_solve(imdl,vh,vi);
0327 show_fem(img2)
0328 
0329 subplot(3,4,9)
0330 [imdl.rec_model, imdl.fwd_model.coarse2fine] = ...
0331     mk_grid_model(img.fwd_model, -1:.1:1, -1:.1:1, 0:.5:3,'outside');
0332 %v = imdl.fwd_model.coarse2fine'*(img.elem_data-1)+1;
0333 %img2 = mk_image(imdl.rec_model, v);
0334 img2= inv_solve(imdl,vh,vi);
0335 show_fem(img2)
0336 
0337 subplot(3,4,10)
0338    l48 = linspace(0.5,48.5,24+1); 
0339    l32 =-linspace(0.5,32.5,16+1) + 20; 
0340    rmfn = @(xyz) vecnorm(xyz,2,2)>40;
0341    fmdl = mk_grid_model([],l48,l32,l48,rmfn);
0342 show_fem(fmdl)