mdl_slice_mesher

MDL_SLICE_MESHER A slice of a 3D FEM as a 2D FEM 
 img2d = mdl_slice_mesher(mdl3d,level) returns a 2D FEM model MDL2D 
 suitable for viewing with SHOW_FEM representing a cut through MDL3D at
 LEVEL. 
 Note that where the intersection of an element of MDL3D and the LEVEL is
 a quadrangle, this will be represented as two triangles in IMG2D. 
 Faces of MDL3D co-planar with LEVEL will be assigned an avarage of the
 values of the two elements that share them. 

 [img2d ptch] = mdl_slice_mesher(...) also returns a struct PTCH suitable
 for use with the PATCH function. It offers the advantage of displaying
 both quad and tri elements. Colors can be controlled by adding
 additional arguments passed to calc_colours:
 [img2d ptch] = mdl_slice_mesher(mdl3d,level, ... )
 
 Inputs:
   MDL3D  - an EIDORS fwd_model or img struct with elem_data
   LEVEL  - any single-slice specification accepted by LEVEL_MODEL_SLICE

 Additional options can be specified as fields of MDL3D.mdl_slice_mesher
 or MDL3D.fwd_model.mdl_slice_mesher:
       .interp_elems (defualt: true) : average element values for faces
                                       co-planar with LEVEL

 To control the transparency use transparency_tresh (see CALC_COLOURS for
 details), e.g.:
    img2d.calc_colours.transparency_thresh = -1; (no transperency)
    calc_colours('transparency_thresh', 0.25); (some transparency)

 See also: LEVEL_MODEL_SLICE, SHOW_FEM, MDL_SLICE_MAPPER, SHOW_3D_SLICES,
           CROP_MODEL, CALC_COLOURS, PATCH

0001 function [nimg out] = mdl_slice_mesher(fmdl,level,varargin)
0002 %MDL_SLICE_MESHER A slice of a 3D FEM as a 2D FEM
0003 % img2d = mdl_slice_mesher(mdl3d,level) returns a 2D FEM model MDL2D
0004 % suitable for viewing with SHOW_FEM representing a cut through MDL3D at
0005 % LEVEL.
0006 % Note that where the intersection of an element of MDL3D and the LEVEL is
0007 % a quadrangle, this will be represented as two triangles in IMG2D.
0008 % Faces of MDL3D co-planar with LEVEL will be assigned an avarage of the
0009 % values of the two elements that share them.
0010 %
0011 % [img2d ptch] = mdl_slice_mesher(...) also returns a struct PTCH suitable
0012 % for use with the PATCH function. It offers the advantage of displaying
0013 % both quad and tri elements. Colors can be controlled by adding
0014 % additional arguments passed to calc_colours:
0015 % [img2d ptch] = mdl_slice_mesher(mdl3d,level, ... )
0016 %
0017 % Inputs:
0018 %   MDL3D  - an EIDORS fwd_model or img struct with elem_data
0019 %   LEVEL  - any single-slice specification accepted by LEVEL_MODEL_SLICE
0020 %
0021 % Additional options can be specified as fields of MDL3D.mdl_slice_mesher
0022 % or MDL3D.fwd_model.mdl_slice_mesher:
0023 %       .interp_elems (defualt: true) : average element values for faces
0024 %                                       co-planar with LEVEL
0025 %
0026 % To control the transparency use transparency_tresh (see CALC_COLOURS for
0027 % details), e.g.:
0028 %    img2d.calc_colours.transparency_thresh = -1; (no transperency)
0029 %    calc_colours('transparency_thresh', 0.25); (some transparency)
0030 %
0031 % See also: LEVEL_MODEL_SLICE, SHOW_FEM, MDL_SLICE_MAPPER, SHOW_3D_SLICES,
0032 %           CROP_MODEL, CALC_COLOURS, PATCH
0033 
0034 % (C) 2012-2021 Bartlomiej Grychtol.
0035 % License: GPL version 2 or version 3
0036 % $Id: mdl_slice_mesher.m 6372 2022-05-09 13:37:47Z bgrychtol $
0037 
0038 
0039 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return, end;
0040 
0041 if isempty(varargin)
0042    try
0043       varargin{1}.calc_colours = fmdl.calc_colours;
0044    end
0045 end
0046 
0047 switch fmdl.type
0048     case 'image'  
0049        img = fmdl;
0050        fmdl = fmdl.fwd_model;
0051     case 'fwd_model'
0052        img = mk_image(fmdl,1);
0053     otherwise; error('Unknown object type');
0054 end
0055 
0056 fmdl.mdl_slice_mesher = parse_opt(fmdl);
0057 
0058 opt.cache_obj = {fmdl.nodes, fmdl.elems, fmdl.mdl_slice_mesher, level};
0059 if isfield(fmdl,'electrode');
0060     opt.cache_obj(end+1) = {fmdl.electrode};
0061 end
0062 opt.fstr      = 'mdl_slice_mesher';
0063 opt.log_level = 4;
0064 
0065 [nmdl, f2c, n2n, p_struct] = eidors_cache(@do_mdl_slice_mesher,{fmdl, level},opt);
0066 nimg = build_image(nmdl, f2c, n2n, img);
0067 
0068 switch nargout
0069    case 2
0070       out = draw_patch(p_struct, nimg.fwd_model, get_img_data(img), varargin{:});
0071    case 0
0072       out = draw_patch(p_struct, nimg.fwd_model, get_img_data(img), varargin{:});
0073       cmap_type = calc_colours('cmap_type');
0074       try 
0075          calc_colours('cmap_type',varargin{1}.calc_colours.cmap_type);
0076       end
0077       colormap(calc_colours('colourmap'));
0078       patch(out);
0079       calc_colours('cmap_type',cmap_type);
0080       clear nimg;
0081 end
0082 
0083 
0084 function opt = parse_opt(fmdl)
0085 opt.interp_elems = true;
0086 
0087 if isfield(fmdl,'mdl_slice_mesher')
0088     flds = fieldnames(fmdl.mdl_slice_mesher);
0089     for i = 1:numel(flds)
0090         opt.(flds{i}) = fmdl.mdl_slice_mesher.(flds{i});
0091     end
0092 end
0093 
0094 
0095 % This is a start  of a function to extrude_3d_if_reqd, so that
0096 %   we can show 3d shapes. Currently not working
0097 % fmdl = extrude_3d_if_reqd( fmdl );
0098 function fmdl = extrude_3d_if_reqd( fmdl );
0099    if size(fmdl.nodes,2)==3; return; end
0100    nn = fmdl.nodes; N= size(nn,1);
0101    ee = fmdl.elems; E= size(ee,1);
0102    oN= ones(N,1);
0103    oE= ones(E,1) + N;
0104    fmdl.nodes = [nn,-10*oN; %10 is arbitrary == a guess
0105                  nn,+10*oN]; 
0106    fmdl.elems = [ee(:,[1,2,3,1]) + oE*[0,0,0,1];
0107                  ee(:,[1,2,3,2]) + oE*[1,0,0,1];
0108                  ee(:,[1,2,3,3]) + oE*[1,1,0,1]];
0109 
0110 function [nmdl, f2c, n2n, out] = do_mdl_slice_mesher(fmdl,level)
0111 
0112 mdl = fmdl;
0113 opt.edge2elem = true;
0114 opt.node2elem = true;
0115 mdl = fix_model(mdl,opt);
0116 edges = mdl.edges;
0117 edge2elem = mdl.edge2elem;
0118 tmp = mdl; 
0119 tmp = level_model_slice( tmp, level, 1 );
0120 [nodeval nodedist] = nodes_above_or_below(tmp,0);
0121 % find which edges are on electrodes
0122 e_nodes = zeros(length(mdl.nodes),1);
0123 try
0124    for i = 1:length(mdl.electrode)
0125       e_nodes(mdl.electrode(i).nodes) = i;
0126       if isfield(mdl.electrode(i),'faces')
0127          e_nodes(unique(mdl.electrode(i).faces)) = i;
0128       end
0129    end
0130 end
0131 e_edges = (sum(e_nodes(edges),2)/ 2)  .* (e_nodes(edges(:,1)) == e_nodes(edges(:,2)));
0132 %% crossed edges
0133 % exclude edges on plane (dealt with later)
0134 idx = (sum(nodeval(edges),2) == 0) & (nodeval(edges(:,1)) ~= 0) ; 
0135 dist = (nodedist(edges(idx,2)) - nodedist(edges(idx,1)));
0136 t = -nodedist(edges(idx,1))./dist;
0137 % new nodes along the edges
0138 nodes = mdl.nodes(edges(idx,1),:) + ...
0139     repmat(t,1,3).*(mdl.nodes(edges(idx,2),:) - mdl.nodes(edges(idx,1),:));
0140 
0141 n2n = sparse(edges(idx,:),uint32((1:length(t))' * ones(1,2)),[1-t,t],size(mdl.nodes,1),size(nodes,1));
0142 
0143 % nn indexes the just-created nodes, els indexes elements
0144 if any(idx)
0145     [nn els] = find(edge2elem(idx,:));
0146 else
0147     nn = []; els = [];
0148 end
0149 els_edge = els;
0150 
0151 electrode_node = e_edges(idx);
0152 %% crossed nodes
0153 idx = find(nodeval == 0);
0154 ln = length(nodes); %store the size
0155 nodes = [nodes; mdl.nodes(idx,:)]; % add the crossed nodes to the new model
0156 n2n = [n2n, sparse(idx,(1:length(idx)),1,size(mdl.nodes,1),length(idx))];
0157 electrode_node = [electrode_node; e_nodes(idx)];
0158 % nnn indexes mdl.nodes(idx), eee indexes mdl.elems
0159 [nnn eee] = find(mdl.node2elem(idx,:));
0160 nnn = nnn + ln; %make a proper index into nodes
0161 % n1: eee = ueee(n1)
0162 [ueee jnk n1] = unique(eee,'last');
0163 nodes_per_elem = jnk;
0164 nodes_per_elem(2:end) = diff(jnk);
0165 % if an elem has 3 crossed nodes, there must be 2 of them, add both for now
0166 add = find(nodes_per_elem == 3);
0167 if ~isempty(add)
0168     % what to do with faces shared between elements?
0169     addel = ueee(add*[1 1 1])';
0170     els = [els; addel(:)];
0171     for i = 1:length(add)
0172        addnd = nnn(n1 == add(i));
0173        nn = [nn; addnd(:)];
0174     end
0175     [els idx] = sort(els);
0176     nn = nn(idx);
0177 end
0178 % for elems with less than 4 crossed edges -> add crossed nodes if needed
0179 [uels jnk n] = unique(els_edge,'last');
0180 
0181 % only consider elements who have both a crossed node and edge
0182 [idx ia ib] = intersect(ueee, uels);
0183 for i = 1:length(ia)
0184     newnodes = nnn(n1==ia(i));
0185     nn = [nn; newnodes];
0186     els = [els; repmat(uels(ib(i)),length(newnodes),1)];
0187 end
0188 [els idx] = sort(els);
0189 nn = nn(idx);
0190 [uels jnk n] = unique(els,'last');
0191 nodes_per_elem = jnk;
0192 nodes_per_elem(2:end) = diff(jnk);
0193 
0194 n_tri = length(uels) + sum(nodes_per_elem==4);
0195 
0196 nmdl.type = 'fwd_model';
0197 nmdl.nodes = nodes;
0198 nmdl.elems = zeros(n_tri,3);
0199 
0200 if n_tri == 0 
0201     error('EIDORS:NoIntersection',... 
0202         'No intersection found between the cut plane [%.2f %.2f %.2f] and the model.', ...
0203         level(1),level(2),level(3));
0204 end
0205 n_el_data = size(fmdl.elems,1);
0206 f2c = sparse(n_el_data,length(uels));
0207 c = 1;
0208 % TODO: Speed this up
0209 for i = 1:length(uels)
0210     switch nodes_per_elem(i)
0211         case 3
0212             nmdl.elems(c,:) = nn(n==i);
0213             f2c(uels(i),c) = 1;
0214             c = c + 1;
0215         case 4
0216             nds = nn(n==i);
0217             nmdl.elems(c,:) = nds(1:3);
0218             f2c(uels(i),c) = 1;
0219             nmdl.elems(c+1,:) = nds(2:4);
0220             f2c(uels(i),c+1) = 1;
0221             c = c + 2;
0222     end
0223 end
0224 % deal with double elements (from shared faces)
0225 nmdl.elems = sort(nmdl.elems,2);
0226 [nmdl.elems, idx] = sortrows(nmdl.elems);
0227 f2c = f2c(:,idx);
0228 [nmdl.elems, n, idx] = unique(nmdl.elems, 'rows');
0229 if ~fmdl.mdl_slice_mesher.interp_elems
0230     f2c = f2c(:,n);
0231 else
0232     [x y] = find(f2c);
0233     % put all source elements that contribute to destination element on one
0234     % column
0235     f2c = sparse(x,idx,1);
0236     % ensure correct number of columns (happens when the last source element
0237     % doesn't contribute
0238     if size(f2c,1) < n_el_data
0239        f2c(n_el_data,end) = 0;
0240     end
0241     n_src_els = sum(f2c,1);
0242     f2c = f2c * spdiag(1./n_src_els);
0243 end
0244 
0245 % add electrodes
0246 try
0247    for i = 1:length(mdl.electrode)
0248       nmdl.electrode(i) = mdl.electrode(i);
0249       nmdl.electrode(i).nodes = find(electrode_node == i);
0250    end
0251 end
0252 out.elem_map = false(size(mdl.elems,1),1);
0253 out.elem_map(uels) = true;
0254 out.node_map = n2n;
0255 out.els  = els;
0256 out.nn   = nn;
0257 
0258 
0259 function nimg = build_image(nmdl, f2c, n2n, img)
0260 nimg = mk_image(nmdl,1);
0261 if isempty(nimg.elem_data) % plane doesn't cut model
0262     return
0263 end
0264 
0265 img_data = get_img_data(img);
0266 if size(img_data,1) == num_nodes(img.fwd_model)
0267     nimg.node_data = n2n' * img_data;
0268     nimg = rmfield(nimg, 'elem_data');
0269 else
0270     nimg.elem_data = f2c' * img_data;
0271 end
0272 try
0273    nimg.calc_colours = img.calc_colours;
0274 end
0275 
0276 
0277 function out = draw_patch(in, nmdl, img_data, varargin)
0278 
0279 
0280 els  = in.els;
0281 nn   = in.nn;
0282 nodes = nmdl.nodes;
0283 uels = find(in.elem_map);
0284 
0285 out.Vertices = nodes;
0286 out.Faces    = NaN(length(uels),4);
0287 
0288 for i = 1:length(uels)
0289     idx = els == uels(i);
0290     id = find(idx);
0291     nnn = nodes(nn(idx),:);
0292     if nnz(idx)==4
0293         if intersection_test(nnn(1,:),nnn(4,:),nnn(2,:),nnn(3,:))
0294             id(3:4) = id([4 3]);
0295         elseif intersection_test(nnn(1,:),nnn(2,:),nnn(3,:),nnn(4,:))
0296             id(2:3) = id([3 2]);
0297         end
0298     end
0299     out.Faces(i,1:length(id)) = nn(id);   
0300 end
0301 switch length(img_data)
0302     case size(in.elem_map,1)
0303         data = img_data(uels);
0304         out.FaceColor = 'flat';
0305     case size(in.node_map,1)
0306         data = in.node_map' * img_data ;
0307         out.FaceColor = 'interp';
0308     otherwise
0309         error('wrong size of data');
0310 end
0311 [out.FaceVertexCData, scl_data] = calc_colours(data,varargin{:});
0312 try
0313    out.FaceVertexAlphaData = double(abs(scl_data) > varargin{1}.calc_colours.transparency_thresh);
0314    out.FaceAlpha = 'flat';
0315 end
0316 
0317 out.CDataMapping = 'direct';
0318 % colormap(calc_colours('colourmap'));
0319 
0320 
0321 
0322 function res = intersection_test(A,B,C,D)
0323 % checks for intersection of segments AB and CD
0324 % if AB and CD intersect in 3D, then their projections on a 2D plane also
0325 % intersect (or are colliniar).
0326 
0327 % assume they don't
0328 res = false;
0329 
0330 % check for interesection on the 3 cartesian planes
0331 idx = 1:3;
0332 for i = 0:2
0333     id = circshift(idx',i)';
0334     id = id(1:2);
0335     res = res || ( sign(signed_area(A(id), B(id), C(id))) ~= ...
0336                    sign(signed_area(A(id), B(id), D(id)))        );
0337 end
0338 
0339 function a = signed_area(A,B,C)
0340     a = ( B(1) - A(1) ) * ( C(2) - A(2) ) - ...
0341         ( C(1) - A(1) ) * ( B(2) - A(2) );
0342 
0343 
0344 function [nodeval dist] = nodes_above_or_below(mdl,level)
0345 
0346 % Set a model-dependent tolerance
0347 tol = min(max(mdl.nodes) - min(mdl.nodes)) * 1e-5;
0348 tol = max(tol, eps);
0349 dist = mdl.nodes(:,3) - level;
0350 dist(abs(dist) < tol) = 0;
0351 nodeval = sign(dist);
0352 
0353 
0354 function do_unit_test
0355     imdl = mk_common_model('n3r2',[16,2]);
0356     img = mk_image(imdl.fwd_model,1);
0357     load datacom.mat A B;
0358     img.elem_data(A) = 1.2;
0359     img.elem_data(B) = 0.8;
0360     subplot(131)
0361     show_fem(img);
0362     subplot(132)
0363     cla
0364 %     show_fem(img.fwd_model);
0365     hold on
0366 %     slc = mdl_slice_mesher(img, [3 -3 2]);
0367 %     slc.calc_colours.transparency_thresh = -1;
0368 %     show_fem(slc);
0369     slc = mdl_slice_mesher(img, [0 inf inf]);
0370     slc.calc_colours.transparency_thresh = -1;
0371     slc.fwd_model.boundary = slc.fwd_model.elems;
0372     show_fem(slc);
0373     slc = mdl_slice_mesher(img, [inf inf 2]);
0374     slc.calc_colours.transparency_thresh = -1;
0375     slc.fwd_model.boundary = slc.fwd_model.elems;
0376     show_fem(slc,[0 1 0]);
0377     slc = mdl_slice_mesher(img, [inf inf 2.5]);
0378     slc.calc_colours.transparency_thresh = -1;
0379     slc.fwd_model.boundary = slc.fwd_model.elems;
0380     show_fem(slc);
0381     slc = mdl_slice_mesher(img, [inf inf 1.3]);
0382     slc.calc_colours.transparency_thresh = -1;
0383     slc.fwd_model.boundary = slc.fwd_model.elems;
0384     show_fem(slc);
0385     zlim([0 3]);
0386     view(3)
0387     hold off
0388     subplot(133)
0389     hold on
0390     mdl_slice_mesher(img, [0 inf inf]);
0391     mdl_slice_mesher(img, [inf inf 2]);
0392     mdl_slice_mesher(img, [inf inf 2.5]);
0393     mdl_slice_mesher(img, [inf inf 1.3]);
0394     mdl_slice_mesher(img, struct('centre',[-.5,-.5,2],'normal_angle',[1 1 .2]))
0395     axis equal
0396     axis tight
0397     view(3)
0398     
0399     % test multi-column image
0400     img.elem_data(:,2) = img.elem_data;
0401     slc = mdl_slice_mesher(img, [0 inf inf]);