0001 function [mapping, outside] = mk_approx_c2f( f_mdl, c_mdl );
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050 if ischar(f_mdl) && strcmp(f_mdl, 'UNIT_TEST'); do_unit_test; return; end
0051
0052 if ischar(f_mdl) && strcmp(f_mdl, 'LOAD'); load; return; end
0053
0054 [c_mdl, f_mdl] = assign_defaults( c_mdl, f_mdl );
0055
0056 copt.cache_obj = cache_obj(c_mdl, f_mdl);
0057 copt.fstr = 'mk_coarse_fine_mapping';
0058
0059
0060 [mapping, outside] = eidors_cache(@mapping_calc,{f_mdl,c_mdl},copt);
0061
0062
0063
0064 function [mapping, outside] = mapping_calc(f_mdl, c_mdl)
0065 f_mdl= offset_and_project( f_mdl, c_mdl);
0066 z_depth = c_mdl.mk_coarse_fine_mapping.z_depth;
0067
0068 f_elems = all_contained_elems( f_mdl, c_mdl, z_depth);
0069 mapping = contained_elems_i( f_mdl, c_mdl, f_elems, z_depth);
0070
0071 if isfield(c_mdl,'coarse2fine')
0072 mapping = mapping*c_mdl.coarse2fine;
0073 end
0074
0075
0076 if nargout>1;
0077 outside = 1 - sum(mapping,2);
0078 end
0079
0080
0081 function c_obj = cache_obj(c_mdl, f_mdl)
0082 c_obj = {c_mdl.nodes, c_mdl.elems, c_mdl.mk_coarse_fine_mapping, ...
0083 f_mdl.nodes, f_mdl.elems, f_mdl.interp_mesh};
0084
0085
0086
0087 function c_elems = all_contained_elems( fm, cm, z_depth)
0088 [nf,ef]= size(fm.elems);
0089 [nc,ec]= size(cm.elems);
0090 fm_pts = zeros(nf*ef,3);
0091
0092
0093
0094
0095 s_fac= .9;
0096 for dim= 1:ef-1
0097
0098 fm_pt= reshape(fm.nodes(fm.elems,dim),nf,ef);
0099 fm_ctr= mean(fm_pt,2)*ones(1,ef);
0100 fm_pt = s_fac*(fm_pt-fm_ctr) + fm_ctr;
0101 fm_pts(:,dim) = fm_pt(:);
0102 end
0103
0104 tsn= search_fm_pts_in_cm(cm, fm_pts, z_depth);
0105 tsn= reshape( tsn, [], ef);
0106
0107
0108
0109 c_elems= all(diff(tsn,1,2)==0,2) .* tsn(:,1);
0110 c_elems(all(tsn==-1,2))= -1;
0111
0112
0113
0114
0115
0116 function tsn= search_fm_pts_in_cm(cm, fm_pts, z_depth);
0117 dc= size(cm.elems,2)-1;
0118 df= size(fm_pts,2);
0119
0120 tsn= -ones(size(fm_pts,1),1);
0121 not_oor= (tsn==-1);
0122
0123 if dc==2
0124
0125 if df==3
0126
0127 not_oor= not_oor & any( abs(fm_pts(:,3) ) <= z_depth , 2);
0128 end
0129 dims=1:2;
0130
0131 elseif dc==3
0132
0133 dims=1:3;
0134
0135 else
0136 error('coarse model must be 2 or 3D');
0137 end
0138
0139 tsn(not_oor)= tsearchn(cm.nodes(:,dims), cm.elems, fm_pts(not_oor,dims));
0140 tsn(isnan(tsn))= -1;
0141
0142
0143
0144
0145
0146
0147 function c_elems = contained_elems_i( fm, cm, idx, z_depth)
0148 [nc,dc]= size(cm.elems);
0149 [nf,df]= size(fm.elems);
0150
0151 fidx= find(idx==0);
0152 l_fidx= length(fidx);
0153
0154 c_e_i= []; c_e_j=[]; c_e_v=[];
0155
0156
0157
0158 n_interp = 7-df;
0159 interp_mdl.type= 'fwd_model';
0160 interp_mdl.nodes= fm.nodes;
0161 interp_mdl.elems= fm.elems(fidx,:);
0162 interp_mdl.interp_mesh.n_interp = fm.interp_mesh.n_interp;
0163
0164
0165 fm_pts = interp_mesh( interp_mdl, n_interp);
0166 l_interp = size(fm_pts,3);
0167
0168 fm_pts = permute( fm_pts, [3,1,2]);
0169 fm_pts = reshape(fm_pts, [], df-1);
0170
0171 tsn= search_fm_pts_in_cm(cm, fm_pts, z_depth);
0172
0173 tsn_idx= ones(l_interp,1)*fidx(:)';
0174 tsn_idx= tsn_idx(:);
0175
0176 outside_idx= tsn==-1;
0177 tsn(outside_idx) = [];
0178 tsn_idx(outside_idx) = [];
0179
0180 in_idx= find(idx<=0);
0181 ridx= 1:nf; ridx(in_idx)=[];
0182 idx(in_idx)=[];
0183
0184
0185
0186 c_elems = sparse(ridx,idx,1,nf,nc) + ...
0187 sparse(tsn_idx,tsn,1,nf,nc)/l_interp;
0188
0189
0190
0191 function xyz = interpxyz( xyzmin, xyzmax, n_interp)
0192 xyzdelta= xyzmax - xyzmin;
0193 xyz_interp = 1 + floor(n_interp * xyzdelta / max(xyzdelta) );
0194 xspace = linspace(xyzmin(1), xyzmax(1), xyz_interp(1) );
0195 yspace = linspace(xyzmin(2), xyzmax(2), xyz_interp(2) );
0196 zspace = linspace(xyzmin(3), xyzmax(3), xyz_interp(3) );
0197 [xx3,yy3,zz3] = ndgrid( xspace, yspace, zspace );
0198 xyz= [xx3(:), yy3(:), zz3(:)];
0199
0200
0201 function f_mdl= offset_and_project( f_mdl, c_mdl)
0202 [fn,fd]= size(f_mdl.nodes);
0203 T= c_mdl.mk_coarse_fine_mapping.f2c_offset;
0204 M= c_mdl.mk_coarse_fine_mapping.f2c_project;
0205
0206 f_mdl.nodes= (M*( f_mdl.nodes - ones(fn,1)*T )')';
0207
0208 function [c_mdl f_mdl] = assign_defaults( c_mdl, f_mdl )
0209 [fn,fd]= size(f_mdl.nodes);
0210 try
0211 c_mdl.mk_coarse_fine_mapping.f2c_offset;
0212 catch
0213 c_mdl.mk_coarse_fine_mapping.f2c_offset= zeros(1,fd);
0214 end
0215 try
0216 c_mdl.mk_coarse_fine_mapping.f2c_project;
0217 catch
0218 c_mdl.mk_coarse_fine_mapping.f2c_project= speye(fd);
0219 end
0220 try
0221 c_mdl.mk_coarse_fine_mapping.z_depth;
0222 catch
0223 c_mdl.mk_coarse_fine_mapping.z_depth= inf;
0224 end
0225 try
0226 f_mdl.interp_mesh.n_interp;
0227
0228
0229 catch
0230 f_mdl.interp_mesh.n_interp = 7 - size(f_mdl.elems,2);
0231 end
0232
0233
0234 function do_unit_test
0235 fmdl = mk_circ_tank(2,[],2); fmdl.nodes = fmdl.nodes*2;
0236 cmdl = mk_circ_tank(2,[],2); cmdl.nodes = cmdl.nodes*2;
0237 c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0238 unit_test_cmp('t1',c2f,eye(16),1e-15)
0239
0240 fmdl = mk_circ_tank(3,[],2);
0241 fmdl.nodes = fmdl.nodes*3;
0242 c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0243 unit_test_cmp('t2 analytic',c2f,[eye(16);zeros(20,16)],1e-15)
0244 c2f = mk_approx_c2f( fmdl, cmdl);
0245 unit_test_cmp('t2 approx',c2f,[eye(16);zeros(20,16)],1e-15)
0246
0247 fmdl = mk_circ_tank(2,[],2); fmdl.nodes = fmdl.nodes*2;
0248 cmdl = mk_circ_tank(1,[],2); cmdl.nodes = cmdl.nodes*1;
0249 c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0250 unit_test_cmp('t3 analytic',c2f,[eye(4);zeros(12,4)],1e-15)
0251 c2f = mk_approx_c2f( fmdl, cmdl);
0252 unit_test_cmp('t3 approx',c2f,[eye(4);zeros(12,4)],1e-15)
0253
0254 fac = 0.8;
0255 cmdl = mk_circ_tank(1,[],2); cmdl.nodes = cmdl.nodes*fac;
0256 c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0257 unit_test_cmp('t4 analytic',c2f,[eye(4)*fac^2;zeros(12,4)],1e-15)
0258 c2f = mk_approx_c2f( fmdl, cmdl);
0259 unit_test_cmp('t4 approx',c2f,[eye(4)*2/3;zeros(12,4)],1e-15)
0260
0261 fac=1.2;
0262 cmdl = mk_circ_tank(1,[],2); cmdl.nodes = cmdl.nodes*fac;
0263 c2f = mk_approx_c2f( fmdl, cmdl);
0264 unit_test_cmp('t5',c2f,[eye(4);eye(4)/3;kron(eye(4),[1;1])/15],1e-15);
0265
0266 fmdl = mk_circ_tank(10,[],2);
0267 cmdl = mk_circ_tank(8,[],2);
0268 c2f = mk_approx_c2f( fmdl, cmdl);
0269 unit_test_cmp('t6',sum(c2f'),ones(1,size(c2f,1)),1e-14);
0270
0271 cmdl.nodes = cmdl.nodes*0.95;
0272
0273 c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0274
0275 function load
0276
0277
0278 electrodes_per_plane = 16;
0279 number_of_planes = 2;
0280 tank_radius = 0.2;
0281 tank_height = 0.5;
0282 fine_mdl = ng_mk_cyl_models([tank_height,tank_radius],...
0283 [electrodes_per_plane,0.15,0.35],[0.01]);
0284
0285
0286 coarse_mdl_maxh = 0.07;
0287 coarse_mdl = ng_mk_cyl_models([tank_height,tank_radius,coarse_mdl_maxh],[0],[]);
0288
0289 disp('Calculating coarse2fine mapping ...');
0290 inv3d.fwd_model.coarse2fine = ...
0291 mk_coarse_fine_mapping( fine_mdl, coarse_mdl);
0292 disp(' ... done');