0001 function [nimg out] = mdl_slice_mesher(fmdl,level,varargin)
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 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return, end;
0041
0042 if isempty(varargin)
0043 try
0044 varargin{1}.calc_colours = fmdl.calc_colours;
0045 end
0046 end
0047
0048 switch fmdl.type
0049 case 'image'
0050 img = fmdl;
0051 fmdl = fmdl.fwd_model;
0052 case 'fwd_model'
0053 img = mk_image(fmdl,1);
0054 otherwise; error('Unknown object type');
0055 end
0056 opt.cache_obj = {fmdl.nodes, fmdl.elems, level};
0057 if isfield(fmdl,'electrode');
0058 opt.cache_obj(end+1) = {fmdl.electrode};
0059 end
0060 opt.fstr = 'mdl_slice_mesher';
0061 opt.log_level = 4;
0062
0063 [nmdl f2c p_struct] = eidors_cache(@do_mdl_slice_mesher,{fmdl, level},opt);
0064 nimg = build_image(nmdl, f2c, img);
0065
0066 switch nargout
0067 case 2
0068 out = draw_patch(p_struct, nimg.fwd_model, get_img_data(img), varargin{:});
0069 case 0
0070 out = draw_patch(p_struct, nimg.fwd_model, get_img_data(img), varargin{:});
0071 cmap_type = calc_colours('cmap_type');
0072 try
0073 calc_colours('cmap_type',varargin{1}.calc_colours.cmap_type);
0074 end
0075 colormap(calc_colours('colourmap'));
0076 patch(out);
0077 calc_colours('cmap_type',cmap_type);
0078 clear nimg;
0079 end
0080
0081
0082
0083
0084 function fmdl = extrude_3d_if_reqd( fmdl );
0085 if size(fmdl.nodes,2)==3; return; end
0086 nn = fmdl.nodes; N= size(nn,1);
0087 ee = fmdl.elems; E= size(ee,1);
0088 oN= ones(N,1);
0089 oE= ones(E,1) + N;
0090 fmdl.nodes = [nn,-10*oN;
0091 nn,+10*oN];
0092 fmdl.elems = [ee(:,[1,2,3,1]) + oE*[0,0,0,1];
0093 ee(:,[1,2,3,2]) + oE*[1,0,0,1];
0094 ee(:,[1,2,3,3]) + oE*[1,1,0,1]];
0095
0096 function [nmdl f2c out] = do_mdl_slice_mesher(fmdl,level)
0097
0098 mdl = fmdl;
0099 opt.edge2elem = true;
0100 opt.node2elem = true;
0101 mdl = fix_model(mdl,opt);
0102 edges = mdl.edges;
0103 edge2elem = mdl.edge2elem;
0104 tmp = mdl;
0105 tmp.nodes = level_model( tmp, level )';
0106 [nodeval nodedist] = nodes_above_or_below(tmp,0);
0107
0108 e_nodes = zeros(length(mdl.nodes),1);
0109 try
0110 for i = 1:length(mdl.electrode)
0111 e_nodes(mdl.electrode(i).nodes) = i;
0112 end
0113 end
0114 e_edges = (sum(e_nodes(edges),2)/ 2) .* (e_nodes(edges(:,1)) == e_nodes(edges(:,2)));
0115
0116
0117 idx = (sum(nodeval(edges),2) == 0) & (nodeval(edges(:,1)) ~= 0) ;
0118 dist = (nodedist(edges(idx,2)) - nodedist(edges(idx,1)));
0119 t = -nodedist(edges(idx,1))./dist;
0120
0121 nodes = mdl.nodes(edges(idx,1),:) + ...
0122 repmat(t,1,3).*(mdl.nodes(edges(idx,2),:) - mdl.nodes(edges(idx,1),:));
0123
0124 if any(idx)
0125 [nn els] = find(edge2elem(idx,:));
0126 else
0127 nn = []; els = [];
0128 end
0129 els_edge = els;
0130
0131 electrode_node = e_edges(idx);
0132
0133 idx = nodeval == 0;
0134 ln = length(nodes);
0135 nodes = [nodes; mdl.nodes(idx,:)];
0136 electrode_node = [electrode_node; e_nodes(idx)];
0137
0138 [nnn eee] = find(mdl.node2elem(idx,:));
0139 nnn = nnn + ln;
0140
0141 [ueee jnk n1] = unique(eee,'last');
0142 nodes_per_elem = jnk;
0143 nodes_per_elem(2:end) = diff(jnk);
0144
0145 add = find(nodes_per_elem == 3);
0146 if ~isempty(add)
0147
0148 addel = ueee(add*[1 1 1])';
0149 els = [els; addel(:)];
0150 for i = 1:length(add)
0151 addnd = nnn(n1 == add(i));
0152 nn = [nn; addnd(:)];
0153 end
0154 [els idx] = sort(els);
0155 nn = nn(idx);
0156 end
0157
0158 [uels jnk n] = unique(els_edge,'last');
0159
0160
0161 [idx ia ib] = intersect(ueee, uels);
0162 for i = 1:length(ia)
0163 newnodes = nnn(n1==ia(i));
0164 nn = [nn; newnodes];
0165 els = [els; repmat(uels(ib(i)),length(newnodes),1)];
0166 end
0167 [els idx] = sort(els);
0168 nn = nn(idx);
0169 [uels jnk n] = unique(els,'last');
0170 nodes_per_elem = jnk;
0171 nodes_per_elem(2:end) = diff(jnk);
0172
0173 n_tri = length(uels) + sum(nodes_per_elem==4);
0174
0175 nmdl.type = 'fwd_model';
0176 nmdl.nodes = nodes;
0177 nmdl.elems = zeros(n_tri,3);
0178
0179 if n_tri == 0
0180 error('EIDORS:NoIntersection',...
0181 'No intersection found between the cut plane [%.2f %.2f %.2f] and the model.', ...
0182 level(1),level(2),level(3));
0183 end
0184 n_el_data = size(fmdl.elems,1);
0185 f2c = sparse(n_el_data,length(uels));
0186 c = 1;
0187
0188 for i = 1:length(uels)
0189 switch nodes_per_elem(i)
0190 case 3
0191 nmdl.elems(c,:) = nn(n==i);
0192 f2c(uels(i),c) = 1;
0193 c = c + 1;
0194 case 4
0195 nds = nn(n==i);
0196 nmdl.elems(c,:) = nds(1:3);
0197 f2c(uels(i),c) = 1;
0198 nmdl.elems(c+1,:) = nds(2:4);
0199 f2c(uels(i),c+1) = 1;
0200 c = c + 2;
0201 end
0202 end
0203
0204 nmdl.elems = sort(nmdl.elems,2);
0205 [nmdl.elems idx] = sortrows(nmdl.elems);
0206 f2c = f2c(:,idx);
0207 [nmdl.elems n idx] = unique(nmdl.elems, 'rows');
0208 [x y] = find(f2c);
0209
0210
0211 f2c = sparse(x,idx,1);
0212
0213
0214 if size(f2c,1) < n_el_data
0215 f2c(n_el_data,end) = 0;
0216 end
0217 n_src_els = sum(f2c,1);
0218 f2c = f2c * spdiag(1./n_src_els);
0219
0220
0221
0222 try
0223 for i = 1:length(mdl.electrode)
0224 nmdl.electrode(i) = mdl.electrode(i);
0225 nmdl.electrode(i).nodes = find(electrode_node == i);
0226 end
0227 end
0228
0229 out.uels = uels;
0230 out.els = els;
0231 out.nn = nn;
0232
0233
0234 function nimg = build_image(nmdl, f2c, img)
0235 nimg = mk_image(nmdl,1);
0236 if isempty(nimg.elem_data)
0237 return
0238 end
0239 nimg.elem_data = (get_img_data(img)' * f2c)';
0240 try
0241 nimg.calc_colours = img.calc_colours;
0242 end
0243
0244
0245 function out = draw_patch(in, nmdl, elem_data, varargin)
0246
0247 uels = in.uels;
0248 els = in.els;
0249 nn = in.nn;
0250 nodes = nmdl.nodes;
0251
0252 img.elem_data = elem_data;
0253
0254 out.Vertices = nodes;
0255 out.Faces = NaN(length(uels),4);
0256 ed = zeros(length(uels),1);
0257
0258
0259 for i = 1:length(uels)
0260 idx = els == uels(i);
0261 id = find(idx);
0262 nnn = nodes(nn(idx),:);
0263 if nnz(idx)==4
0264 if intersection_test(nnn(1,:),nnn(4,:),nnn(2,:),nnn(3,:))
0265 id(3:4) = id([4 3]);
0266 elseif intersection_test(nnn(1,:),nnn(2,:),nnn(3,:),nnn(4,:))
0267 id(2:3) = id([3 2]);
0268 end
0269 end
0270
0271 out.Faces(i,1:length(id)) = nn(id);
0272 ed(i) = img.elem_data(uels(i));
0273 end
0274 [out.FaceVertexCData scl_data] = calc_colours(ed,varargin{:});
0275 try
0276 out.FaceVertexAlphaData = double(abs(scl_data) > varargin{1}.calc_colours.transparency_thresh);
0277 out.FaceAlpha = 'flat';
0278 end
0279 out.FaceColor = 'flat';
0280 out.CDataMapping = 'direct';
0281
0282
0283
0284 function res = intersection_test(A,B,C,D)
0285
0286
0287
0288
0289
0290 res = false;
0291
0292
0293 idx = 1:3;
0294 for i = 0:2
0295 id = circshift(idx',i)';
0296 id = id(1:2);
0297 res = res || ( sign(signed_area(A(id), B(id), C(id))) ~= ...
0298 sign(signed_area(A(id), B(id), D(id))) );
0299 end
0300
0301 function a = signed_area(A,B,C)
0302 a = ( B(1) - A(1) ) * ( C(2) - A(2) ) - ...
0303 ( C(1) - A(1) ) * ( B(2) - A(2) );
0304
0305
0306 function [nodeval dist] = nodes_above_or_below(mdl,level)
0307
0308
0309 tol = min(max(mdl.nodes) - min(mdl.nodes)) * 1e-10;
0310
0311 dist = mdl.nodes(:,3) - level;
0312 dist(abs(dist) < tol) = 0;
0313 nodeval = sign(dist);
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323
0324 function NODE= level_model( fwd_model, level )
0325
0326 vtx= fwd_model.nodes;
0327 [nn, dims] = size(vtx);
0328 if dims ==2
0329 NODE= vtx';
0330 return;
0331 end
0332
0333
0334
0335 level( isinf(level) | isnan(level) ) = realmax;
0336 level( level==0 ) = 1e-10;
0337
0338
0339
0340 invlev= 1./level;
0341 ctr= invlev / sum( invlev.^2 );
0342
0343
0344
0345 [jnk, s_ax]= sort( - abs(level - ctr) );
0346 v1= [0,0,0]; v1(s_ax(1))= level(s_ax(1));
0347 v1= v1 - ctr;
0348 v1= v1 / norm(v1);
0349
0350
0351 v2= [0,0,0]; v2(s_ax(2))= level(s_ax(2));
0352 v2= v2 - ctr;
0353 v2= v2 / norm(v2);
0354 v3= cross(v1,v2);
0355
0356
0357 v2= cross(v1,v3);
0358
0359
0360 v1= v1 * (1-2*(sum(v1)<0));
0361 v2= v2 * (1-2*(sum(v2)<0));
0362 v3= v3 * (1-2*(sum(v3)<0));
0363
0364 NODE= [v1;v2;v3] * (vtx' - ctr'*ones(1,nn) );
0365
0366
0367 function do_unit_test
0368 imdl = mk_common_model('n3r2',[16,2]);
0369 img = mk_image(imdl.fwd_model,1);
0370 load datacom.mat A B;
0371 img.elem_data(A) = 1.2;
0372 img.elem_data(B) = 0.8;
0373 subplot(131)
0374 show_fem(img);
0375 subplot(132)
0376 cla
0377
0378 hold on
0379
0380
0381
0382 slc = mdl_slice_mesher(img, [0 inf inf]);
0383 slc.calc_colours.transparency_thresh = -1;
0384 slc.fwd_model.boundary = slc.fwd_model.elems;
0385 show_fem(slc);
0386 slc = mdl_slice_mesher(img, [inf inf 2]);
0387 slc.calc_colours.transparency_thresh = -1;
0388 slc.fwd_model.boundary = slc.fwd_model.elems;
0389 show_fem(slc,[0 1 0]);
0390 slc = mdl_slice_mesher(img, [inf inf 2.5]);
0391 slc.calc_colours.transparency_thresh = -1;
0392 slc.fwd_model.boundary = slc.fwd_model.elems;
0393 show_fem(slc);
0394 slc = mdl_slice_mesher(img, [inf inf 1.3]);
0395 slc.calc_colours.transparency_thresh = -1;
0396 slc.fwd_model.boundary = slc.fwd_model.elems;
0397 show_fem(slc);
0398 zlim([0 3]);
0399 view(3)
0400 hold off
0401 subplot(133)
0402 hold on
0403 mdl_slice_mesher(img, [0 inf inf]);
0404 mdl_slice_mesher(img, [inf inf 2]);
0405 mdl_slice_mesher(img, [inf inf 2.5]);
0406 mdl_slice_mesher(img, [inf inf 1.3]);
0407 axis equal
0408 axis tight
0409 view(3)
0410
0411
0412 img.elem_data(:,2) = img.elem_data;
0413 slc = mdl_slice_mesher(img, [0 inf inf]);