0001 function c2f = mk_tri_c2f(fmdl,rmdl,opt)
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 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'), do_unit_test; return; end
0046 if nargin < 3
0047 opt = struct();
0048 end
0049
0050 f_elems = size(fmdl.elems,1);
0051 r_elems = size(rmdl.elems,1);
0052
0053 c2f = sparse(f_elems,r_elems);
0054 [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl);
0055
0056 if ~(any(fmdl_idx) && any(rmdl_idx))
0057 eidors_msg('@@: models do not overlap, returning all-zeros');
0058 return
0059 end
0060
0061 [fmdl,rmdl] = center_scale_models(fmdl,rmdl, opt);
0062
0063 opt = parse_opts(fmdl,rmdl, opt);
0064
0065
0066 copt.fstr = 'mk_tri_c2f';
0067
0068 c2f(fmdl_idx,rmdl_idx) = eidors_cache(@do_mk_tri_c2f,{fmdl,rmdl,opt},copt);
0069 end
0070
0071 function c2f = do_mk_tri_c2f(fmdl,rmdl,opt)
0072 DEBUG = eidors_debug('query','mk_tri_c2f');
0073
0074 c2f = sparse(0,0);
0075 progress_msg('Prepare fine model...');
0076 fmdl = prepare_tri_mdl(fmdl);
0077 progress_msg(Inf);
0078
0079 progress_msg('Prepare course model...');
0080 rmdl = prepare_tri_mdl(rmdl);
0081 progress_msg(Inf);
0082
0083 progress_msg('Find edge2edge intersections...');
0084 [intpts,fedge2redge,fedge2intpts,redge2intpts] = ...
0085 find_edge2edge_intersections( fmdl.edges,fmdl.nodes,...
0086 rmdl.edges,rmdl.nodes);
0087 progress_msg(sprintf('Found %d',size(intpts,1)),Inf);
0088
0089
0090 progress_msg('Find c_nodes in f_elems...');
0091 rnode2ftri = point_in_triangle(rmdl.nodes, fmdl.elems, fmdl.nodes, opt.tol_node2tri);
0092 progress_msg(sprintf('Found %d', nnz(rnode2ftri)), Inf);
0093
0094 progress_msg('Find c_elems in f_elems...')
0095 rtri_in_ftri = (double(rmdl.node2elem') * rnode2ftri) == 3;
0096 progress_msg(sprintf('Found %d',nnz(rtri_in_ftri)), Inf);
0097
0098 progress_msg('Find f_nodes in c_elems...');
0099 fnode2rtri = point_in_triangle(fmdl.nodes, rmdl.elems, rmdl.nodes, opt.tol_node2tri);
0100 progress_msg(sprintf('Found %d', nnz(fnode2rtri)), Inf);
0101
0102 progress_msg('Find f_elems in c_elems...')
0103 ftri_in_rtri = (double(fmdl.node2elem') * fnode2rtri) == 3;
0104 progress_msg(sprintf('Found %d',nnz(ftri_in_rtri)), Inf);
0105
0106 progress_msg('Find total intersections...');
0107 rtri2ftri = double(rmdl.edge2elem') * fedge2redge' * fmdl.edge2elem;
0108
0109 rtri2ftri = rtri2ftri & ~rtri_in_ftri & ~ftri_in_rtri';
0110 progress_msg(sprintf('Found %d',nnz(rtri2ftri)), Inf);
0111
0112 progress_msg('Calculate intersection volumes...');
0113
0114 rtri2intpt = logical(rmdl.edge2elem'*redge2intpts)';
0115 ftri2intpt = logical(fmdl.edge2elem'*fedge2intpts)';
0116
0117 rtri_todo = find(sum(rtri2ftri,2)>0);
0118 C = []; F = []; V = [];
0119
0120 id = 0; N = length(rtri_todo);
0121 mint = ceil(N/100);
0122 for v = rtri_todo'
0123 id = id+1;
0124 if mod(id,mint)==0, progress_msg(id/N); end
0125 tri_todo = find(rtri2ftri(v,:));
0126
0127 common_intpts = bsxfun(@and,rtri2intpt(:,v), ftri2intpt(:,tri_todo));
0128
0129 f_nodes = bsxfun(@and,fnode2rtri(:,v), fmdl.node2elem(:,tri_todo));
0130 r_nodes = bsxfun(@and,rnode2ftri(:,tri_todo), rmdl.node2elem(:,v));
0131
0132 C = [C; v*ones(numel(tri_todo),1)];
0133 F = [F; tri_todo'];
0134 last_v = numel(V);
0135 V = [V; zeros(numel(tri_todo),1)];
0136
0137 for t = 1:numel(tri_todo)
0138 pts = [ intpts(common_intpts(:,t),:);
0139 fmdl.nodes(f_nodes(:,t),:);
0140 rmdl.nodes(r_nodes(:,t),:)];
0141 last_v = last_v + 1;
0142 if size(pts,1) < 3, continue, end
0143 try
0144
0145
0146 ctr = mean(pts);
0147 pts = bsxfun(@minus,pts,ctr);
0148 scale = max(abs(pts(:)));
0149 if scale == 0
0150 continue
0151 end
0152
0153 pts = pts ./ scale;
0154
0155
0156 [K, V(last_v)] = convhulln(pts,{'Qt Pp Qs'});
0157 V(last_v) = max(V(last_v),0);
0158 V(last_v) = V(last_v) * scale^2;
0159 catch err
0160 ok = false;
0161 switch err.identifier
0162 case {'MATLAB:qhullmx:DegenerateData', 'MATLAB:qhullmx:UndefinedError'}
0163
0164
0165
0166 u = uniquetol(pts*scale,1e-14,'ByRows',true,'DataScale', 1);
0167 ok = ok | (size(u,1) < 3);
0168 if ~ok
0169
0170 cp = bsxfun(@minus, u(2:end,:), u(1,:));
0171 l = sqrt(sum(cp.^2,2));
0172 cp = bsxfun(@rdivide, cp, l);
0173 u = uniquetol(cp,1e-14,'ByRows',true,'DataScale',1);
0174 ok = ok | size(u,1) == 1;
0175 end
0176 end
0177 if ~ok
0178 if DEBUG || eidors_debug('query','mk_tet_c2f:convhulln');
0179 tri.nodes = fmdl.nodes;
0180 vox.nodes = rmdl.nodes;
0181 tri.type = 'fwd_model';
0182 vox.type = 'fwd_model';
0183 vox.elems = rmdl.elems(v,:);
0184 vox.boundary = vox.elems;
0185 tri.elems = fmdl.elems(tri_todo(t),:);
0186 clf
0187 show_fem(vox)
0188 hold on
0189 h = show_fem(tri);
0190 set(h,'EdgeColor','b')
0191 pts = bsxfun(@plus,pts*scale,ctr);
0192 plot(pts(:,1),pts(:,2),'o');
0193 hold off
0194 axis auto
0195 keyboard
0196 else
0197 fprintf('\n');
0198 eidors_msg(['convhulln has thrown an error. ' ...
0199 'Enable eidors_debug on mk_tri_c2f and re-run to see a debug plot'],0);
0200 rethrow(err);
0201 end
0202 end
0203 end
0204
0205 end
0206 end
0207 progress_msg(Inf);
0208
0209 c2f = sparse(F,C,V,size(fmdl.elems,1),size(rmdl.elems,1));
0210
0211
0212 try rmdl = rmfield(rmdl,'coarse2fine'); end
0213 c2f = c2f + bsxfun(@times, sparse(rtri_in_ftri), get_elem_volume(rmdl))';
0214
0215
0216 vol = get_elem_volume(fmdl);
0217 c2f = bsxfun(@rdivide,c2f,vol);
0218
0219
0220 ftri_in_rtri(rtri_in_ftri') = 0;
0221
0222
0223 c2f = c2f + ftri_in_rtri;
0224
0225 end
0226
0227
0228 function [pts,FE2CE,FE2pts,CE2pts] = find_edge2edge_intersections(FE,FN,CE,CN)
0229 P1 = FN(FE(:,1),:);
0230 P2 = FN(FE(:,2),:);
0231 P3 = CN(CE(:,1),:);
0232 P4 = CN(CE(:,2),:);
0233
0234 P21 = P2 - P1;
0235 P43 = P4 - P3;
0236
0237 invden = ( bsxfun(@times, P21(:,1), P43(:,2)') - ...
0238 bsxfun(@times, P21(:,2), P43(:,1)') ).^-1;
0239 P13_x = bsxfun(@minus,P1(:,1),P3(:,1)');
0240 P13_y = bsxfun(@minus,P1(:,2),P3(:,2)');
0241 s = ( bsxfun(@times,-P21(:,2), P13_x) + ...
0242 bsxfun(@times, P21(:,1), P13_y)) .* invden;
0243 t = ( bsxfun(@times,-P43(:,2)',P13_x) + ...
0244 bsxfun(@times, P43(:,1)',P13_y)) .* invden;
0245
0246 FE2CE= s >= 0 & s <= 1 & t >= 0 & t <= 1;
0247
0248 [fe, ce] = find(FE2CE);
0249 N_pts = size(fe,1);
0250 pts = zeros(N_pts,2);
0251 for i = 1:N_pts
0252 pts(i,:) = P1(fe(i),:) + t(fe(i),ce(i)) * P21(fe(i),:);
0253 end
0254 FE2CE = sparse(FE2CE);
0255 N_ce = size(CE,1);
0256 N_fe = size(FE,1);
0257 FE2pts = sparse(fe, 1:N_pts, ones(N_pts,1), N_fe, N_pts);
0258 CE2pts = sparse(ce, 1:N_pts, ones(N_pts,1), N_ce, N_pts);
0259
0260
0261 end
0262
0263
0264
0265 function fmdl = prepare_tri_mdl(fmdl)
0266 fmopt.elem2edge = true;
0267 fmopt.edge2elem = true;
0268
0269 fmopt.node2elem = true;
0270
0271 fmopt.linear_reorder = false;
0272 ll = eidors_msg('log_level',1);
0273 fmdl = fix_model(fmdl,fmopt);
0274 eidors_msg('log_level',ll);
0275 fmdl.node2elem = logical(fmdl.node2elem);
0276
0277 end
0278
0279
0280
0281 function [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl)
0282
0283
0284
0285
0286 f_min = min(fmdl.nodes);
0287 f_max = max(fmdl.nodes);
0288 r_min = min(rmdl.nodes);
0289 r_max = max(rmdl.nodes);
0290
0291
0292 f_gt = bsxfun(@gt, fmdl.nodes, r_max);
0293 f_lt = bsxfun(@lt, fmdl.nodes, r_min);
0294 r_gt = bsxfun(@gt, rmdl.nodes, f_max);
0295 r_lt = bsxfun(@lt, rmdl.nodes, f_min);
0296
0297
0298 re_gt = any(reshape(all(reshape(r_gt(rmdl.elems',:),3,[])),[],2),2);
0299 re_lt = any(reshape(all(reshape(r_lt(rmdl.elems',:),3,[])),[],2),2);
0300 fe_gt = any(reshape(all(reshape(f_gt(fmdl.elems',:),3,[])),[],2),2);
0301 fe_lt = any(reshape(all(reshape(f_lt(fmdl.elems',:),3,[])),[],2),2);
0302
0303
0304 rmdl_idx = ~(re_gt | re_lt);
0305 fmdl_idx = ~(fe_gt | fe_lt);
0306
0307
0308 rmdl.elems = rmdl.elems(rmdl_idx,:);
0309 fmdl.elems = fmdl.elems(fmdl_idx,:);
0310
0311
0312 [r_used_nodes,jnk,r_n] = unique(rmdl.elems(:));
0313 [f_used_nodes,jnk,f_n] = unique(fmdl.elems(:));
0314
0315 r_idx = 1:numel(r_used_nodes);
0316 f_idx = 1:numel(f_used_nodes);
0317
0318 rmdl.elems = reshape(r_idx(r_n),[],3);
0319 fmdl.elems = reshape(f_idx(f_n),[],3);
0320
0321 rmdl.nodes = rmdl.nodes(r_used_nodes,:);
0322 fmdl.nodes = fmdl.nodes(f_used_nodes,:);
0323
0324
0325 try, rmdl = rmfield(rmdl,'boundary'); end
0326 try, fmdl = rmfield(fmdl,'boundary'); end
0327
0328 end
0329
0330
0331
0332 function[fmdl,rmdl] = center_scale_models(fmdl,rmdl, opt)
0333 ctr = mean([min(rmdl.nodes);max(rmdl.nodes)]);
0334 rmdl.nodes = bsxfun(@minus,rmdl.nodes,ctr);
0335 fmdl.nodes = bsxfun(@minus,fmdl.nodes,ctr);
0336 if isfield(opt,'do_not_scale') && opt.do_not_scale
0337 return
0338 end
0339 maxnode = min( max(abs(rmdl.nodes(:))), max(abs(fmdl.nodes(:))));
0340 scale = 1/maxnode;
0341 rmdl.nodes = scale*rmdl.nodes;
0342 fmdl.nodes = scale*fmdl.nodes;
0343 eidors_msg('@@ models scaled by %g', scale,2);
0344 end
0345
0346
0347
0348 function opt = parse_opts(fmdl,rmdl, opt)
0349
0350 if ~isfield(opt, 'tol_node2tri');
0351 opt.tol_node2tri = eps;
0352 end
0353
0354 eidors_msg('@@ node2tri tolerance = %g', opt.tol_node2tri,2);
0355
0356 end
0357
0358 function do_unit_test
0359 do_small_test
0360 end
0361
0362 function do_small_test
0363 imdl = mk_common_model('a2c',16);
0364 rmdl = imdl.fwd_model;
0365
0366
0367
0368
0369 imdl = mk_common_model('d2c',16);
0370 fmdl = imdl.fwd_model;
0371
0372
0373
0374
0375
0376 c2f = mk_tri_c2f(fmdl,rmdl);
0377
0378
0379 clf
0380 img1 = mk_image(fmdl,0);
0381 img1.elem_data = sum(c2f,2);
0382 img1.calc_colous.ref_level = .5;
0383 img1.calc_colours.clim = .5;
0384
0385 subplot(121)
0386 show_fem(img1);
0387 img2 = mk_image(rmdl,0);
0388 img2.elem_data = (c2f' * get_elem_volume(fmdl)) ./ get_elem_volume(rmdl);
0389 img2.calc_colous.ref_level = .5;
0390 img2.calc_colours.clim = .55;
0391 subplot(122)
0392 show_fem(img2);
0393
0394 unit_test_cmp('Check C2F size', size(c2f),[length(fmdl.elems), length(rmdl.elems)]);
0395 unit_test_cmp('Check C2F max==1', max(c2f(:)), 1);
0396 unit_test_cmp('Check C2F min==0', min(c2f(:)), 0);
0397
0398 f2c = bsxfun(@rdivide, bsxfun(@times, c2f, get_elem_volume(fmdl))', get_elem_volume(rmdl));
0399 unit_test_cmp('Check F2C max==1', max(sum(f2c,2)), 1, 1e-15);
0400 unit_test_cmp('Check F2C min==0', min(f2c(:)), 0);
0401 end