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

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 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0018 % $Id: mk_grid_model.html 2819 2011-09-07 16:43:11Z aadler $
0019 
0020 if nargin == 3
0021    cmdl = mk_2d_grid(xvec,yvec);
0022 elseif nargin ==4
0023    cmdl = mk_3d_grid(xvec,yvec,zvec);
0024 else
0025    error('check nargin');
0026 end
0027 
0028 if ~isempty( fmdl)
0029    if nargin ==3
0030       c2f= calc_c2f_2d( fmdl, xvec, yvec);
0031    elseif nargin ==4
0032       c2f= calc_c2f_3d( fmdl, xvec, yvec, zvec);
0033    end
0034 end
0035 
0036 function c2f= calc_c2f_2d( fmdl, xvec, yvec);
0037    nef= size( fmdl.elems,1);
0038    c2f= sparse(nef,0);
0039    mdl_pts = interp_mesh( fmdl, 3);
0040    x_pts = squeeze(mdl_pts(:,1,:));
0041    y_pts = squeeze(mdl_pts(:,2,:));
0042    for yi= 1:length(yvec)-1
0043          in_y_pts = y_pts >= yvec(yi) & y_pts < yvec(yi+1);
0044       for xi= 1:length(xvec)-1
0045           in_x_pts =  x_pts >= xvec(xi) & x_pts < xvec(xi+1);
0046           in_pts = mean( in_y_pts & in_x_pts , 2);
0047           c2f = [c2f,sparse(in_pts)];
0048       end
0049    end
0050 
0051 function c2f= calc_c2f_3d( fmdl, xvec, yvec, zvec);
0052 %  c2f= mk_coarse_fine_mapping( fmdl, cmdl);
0053    nef= size( fmdl.elems,1);
0054 %  c2f= sparse(nef,0);
0055    c2fiidx= [];
0056    c2fjidx= [];
0057    c2fdata= [];
0058    jidx= 0;
0059    mdl_pts = interp_mesh( fmdl, 3);
0060    x_pts = squeeze(mdl_pts(:,1,:));
0061    y_pts = squeeze(mdl_pts(:,2,:));
0062    z_pts = squeeze(mdl_pts(:,3,:));
0063    
0064    in_x_pts = calc_in_d_pts( x_pts, xvec);
0065    in_y_pts = calc_in_d_pts( y_pts, yvec);
0066    in_z_pts = calc_in_d_pts( z_pts, zvec);
0067 
0068    for zi= 1:length(zvec)-1
0069       for yi= 1:length(yvec)-1
0070              in_yz_pts = in_y_pts{yi} & in_z_pts{zi};
0071          for xi= 1:length(xvec)-1
0072              in_pts = mean( in_x_pts{xi} & in_yz_pts, 2);
0073              % c2f = [c2f,sparse(in_pts)];
0074              [ii,jj,vv] = find(in_pts);
0075              c2fiidx= [c2fiidx;ii];
0076              c2fjidx= [c2fjidx;jj+jidx]; jidx=jidx+1;
0077              c2fdata= [c2fdata;vv];
0078          end
0079       end
0080    end
0081    c2f= sparse(c2fiidx,c2fjidx,c2fdata, length(in_pts), jidx);
0082 
0083 function cmdl= mk_2d_grid(xvec, yvec);
0084    xlen = length(xvec);
0085    ylen = length(yvec);
0086    cmdl= eidors_obj('fwd_model', ...
0087             sprintf('Grid model %d x %d', xlen, ylen) );
0088 
0089    [x,y]= ndgrid( xvec, yvec);
0090    cmdl.nodes= [x(:),y(:)];
0091    k= 1:xlen-1;
0092    elem_frac = [ k;k+1;k+xlen; ...
0093                  k+1;k+xlen;k+xlen+1];
0094    elem_frac= reshape(elem_frac, 3,[])';
0095    cmdl.elems=  [];
0096    for j=0:ylen-2
0097       cmdl.elems=  [cmdl.elems; elem_frac + xlen*j];
0098    end
0099 
0100    cmdl.boundary = find_boundary( cmdl.elems);
0101 
0102 % assign one single parameter to each square element
0103    e= size(cmdl.elems,1);
0104    params= ceil(( 1:e )/2);
0105    cmdl.coarse2fine = sparse(1:e,params,1,e,max(params));
0106 
0107 function cmdl= mk_3d_grid(xvec, yvec, zvec);
0108    xlen = length(xvec);
0109    ylen = length(yvec);
0110    zlen = length(zvec);
0111    cmdl= eidors_obj('fwd_model', ...
0112             sprintf('Grid model %d x %d x %d', xlen, ylen, zlen) );
0113 
0114    [x,y,z]= ndgrid( xvec, yvec, zvec);
0115    cmdl.nodes= [x(:),y(:),z(:)];
0116    k= 1:xlen-1;
0117    ac = xlen; up = xlen*ylen; % accross vs up
0118    elem_frac = [ k;     k+1 ;  k+ac;   k+up;  ...
0119                  k+1;   k+ac;  k+up;   k+up+1; ...
0120                  k+ac;  k+up;  k+up+1; k+up+ac; ...
0121                  k+1;   k+ac;  k+ac+1; k+up+1; ...
0122                  k+ac;  k+ac+1;k+up+1; k+up+ac; ...
0123                  k+ac+1;k+up+1;k+up+ac;k+up+ac+1];
0124    elem_frac= reshape(elem_frac, 4,[])';
0125 
0126    row_frac =  [];
0127    for j=0:ylen-2
0128       row_frac=  [row_frac; elem_frac + ac*j];
0129    end
0130 
0131    cmdl.elems=  [];
0132    for k=0:zlen-2
0133       cmdl.elems=  [cmdl.elems; row_frac + up*k];
0134    end
0135 
0136    cmdl.boundary = find_boundary( cmdl.elems);
0137 
0138 % assign one single parameter to each square element
0139    e= size(cmdl.elems,1);
0140    params= ceil(( 1:e )/6);
0141    cmdl.coarse2fine = sparse(1:e,params,1,e,max(params));
0142 
0143 
0144 function in_d_pts = calc_in_d_pts( d_pts, dvec);
0145    l1dvec= length(dvec)-1;
0146    in_d_pts = cell(l1dvec,1);
0147    for i= 1:l1dvec
0148       in_d_pts{i} = d_pts >= dvec(i) & d_pts < dvec(i+1);
0149    end