0001 function [cmdl, c2f]= mk_grid_model(fmdl, xvec, yvec, zvec);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
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
0034 cmdl = set_pixel_pos(cmdl,xvec,yvec);
0035
0036
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
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;
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
0082 nef= size( fmdl.elems,1);
0083
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
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
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;
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
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);
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;