Home > eidors > models > mk_grid_model.m

mk_grid_model

PURPOSE ^

MK_GRID_MODEL: Create reconstruction model on pixelated grid

SYNOPSIS ^

function [cmdl, c2f]= mk_grid_model(fmdl, xvec, yvec, zvec);

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Fri 01-Jun-2018 15:59:55 by m2html © 2005