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  - Vector [1x3] of intercepts
          of the slice on the x, y, z axis. To specify a z=2 plane
          parallel to the x,y: use levels= [inf,inf,2]

 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: 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  - Vector [1x3] of intercepts
0020 %          of the slice on the x, y, z axis. To specify a z=2 plane
0021 %          parallel to the x,y: use levels= [inf,inf,2]
0022 %
0023 % To control the transparency use transparency_tresh (see CALC_COLOURS for
0024 % details), e.g.:
0025 %    img2d.calc_colours.transparency_thresh = -1; (no transperency)
0026 %    calc_colours('transparency_thresh', 0.25); (some transparency)
0027 %
0028 % See also: SHOW_FEM, MDL_SLICE_MAPPER, SHOW_3D_SLICES, CROP_MODEL,
0029 %           CALC_COLOURS. PATCH
0030 
0031 % (C) 2012-2015 Bartlomiej Grychtol.
0032 % License: GPL version 2 or version 3
0033 % $Id: mdl_slice_mesher.m 5991 2019-06-27 12:48:02Z aadler $
0034 
0035 % TODO:
0036 %  1. More intuitive cut plane specification
0037 %  2. Support node_data
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 % This is a start  of a function to extrude_3d_if_reqd, so that
0082 %   we can show 3d shapes. Currently not working
0083 % fmdl = extrude_3d_if_reqd( fmdl );
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; %10 is arbitrary == a guess
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 % find which edges are on electrodes
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       if isfield(mdl.electrode(i),'faces')
0113          e_nodes(unique(mdl.electrode(i).faces)) = i;
0114       end
0115    end
0116 end
0117 e_edges = (sum(e_nodes(edges),2)/ 2)  .* (e_nodes(edges(:,1)) == e_nodes(edges(:,2)));
0118 %% crossed edges
0119 % exclude edges on plane (dealt with later)
0120 idx = (sum(nodeval(edges),2) == 0) & (nodeval(edges(:,1)) ~= 0) ; 
0121 dist = (nodedist(edges(idx,2)) - nodedist(edges(idx,1)));
0122 t = -nodedist(edges(idx,1))./dist;
0123 % new nodes along the edges
0124 nodes = mdl.nodes(edges(idx,1),:) + ...
0125     repmat(t,1,3).*(mdl.nodes(edges(idx,2),:) - mdl.nodes(edges(idx,1),:));
0126 % nn indexes the just-created nodes, els indexes elements
0127 if any(idx)
0128     [nn els] = find(edge2elem(idx,:));
0129 else
0130     nn = []; els = [];
0131 end
0132 els_edge = els;
0133 
0134 electrode_node = e_edges(idx);
0135 %% crossed nodes
0136 idx = nodeval == 0;
0137 ln = length(nodes); %store the size
0138 nodes = [nodes; mdl.nodes(idx,:)]; % add the crossed nodes to the new model
0139 electrode_node = [electrode_node; e_nodes(idx)];
0140 % nnn indexes mdl.nodes(idx), eee indexes mdl.elems
0141 [nnn eee] = find(mdl.node2elem(idx,:));
0142 nnn = nnn + ln; %make a proper index into nodes
0143 % n1: eee = ueee(n1)
0144 [ueee jnk n1] = unique(eee,'last');
0145 nodes_per_elem = jnk;
0146 nodes_per_elem(2:end) = diff(jnk);
0147 % if an elem has 3 crossed nodes, there must be 2 of them, add both for now
0148 add = find(nodes_per_elem == 3);
0149 if ~isempty(add)
0150     % what to do with faces shared between elements?
0151     addel = ueee(add*[1 1 1])';
0152     els = [els; addel(:)];
0153     for i = 1:length(add)
0154        addnd = nnn(n1 == add(i));
0155        nn = [nn; addnd(:)];
0156     end
0157     [els idx] = sort(els);
0158     nn = nn(idx);
0159 end
0160 % for elems with less than 4 crossed edges -> add crossed nodes if needed
0161 [uels jnk n] = unique(els_edge,'last');
0162 
0163 % only consider elements who have both a crossed node and edge
0164 [idx ia ib] = intersect(ueee, uels);
0165 for i = 1:length(ia)
0166     newnodes = nnn(n1==ia(i));
0167     nn = [nn; newnodes];
0168     els = [els; repmat(uels(ib(i)),length(newnodes),1)];
0169 end
0170 [els idx] = sort(els);
0171 nn = nn(idx);
0172 [uels jnk n] = unique(els,'last');
0173 nodes_per_elem = jnk;
0174 nodes_per_elem(2:end) = diff(jnk);
0175 
0176 n_tri = length(uels) + sum(nodes_per_elem==4);
0177 
0178 nmdl.type = 'fwd_model';
0179 nmdl.nodes = nodes;
0180 nmdl.elems = zeros(n_tri,3);
0181 
0182 if n_tri == 0 
0183     error('EIDORS:NoIntersection',... 
0184         'No intersection found between the cut plane [%.2f %.2f %.2f] and the model.', ...
0185         level(1),level(2),level(3));
0186 end
0187 n_el_data = size(fmdl.elems,1);
0188 f2c = sparse(n_el_data,length(uels));
0189 c = 1;
0190 % TODO: Speed this up
0191 for i = 1:length(uels)
0192     switch nodes_per_elem(i)
0193         case 3
0194             nmdl.elems(c,:) = nn(n==i);
0195             f2c(uels(i),c) = 1;
0196             c = c + 1;
0197         case 4
0198             nds = nn(n==i);
0199             nmdl.elems(c,:) = nds(1:3);
0200             f2c(uels(i),c) = 1;
0201             nmdl.elems(c+1,:) = nds(2:4);
0202             f2c(uels(i),c+1) = 1;
0203             c = c + 2;
0204     end
0205 end
0206 % deal with double elements (from shared faces)
0207 nmdl.elems = sort(nmdl.elems,2);
0208 [nmdl.elems idx] = sortrows(nmdl.elems);
0209 f2c = f2c(:,idx);
0210 [nmdl.elems n idx] = unique(nmdl.elems, 'rows');
0211 [x y] = find(f2c);
0212 % put all source elements that contribute to destination element on one
0213 % column
0214 f2c = sparse(x,idx,1);
0215 % ensure correct number of columns (happens when the last source element
0216 % doesn't contribute
0217 if size(f2c,1) < n_el_data
0218    f2c(n_el_data,end) = 0;
0219 end
0220 n_src_els = sum(f2c,1);
0221 f2c = f2c * spdiag(1./n_src_els);
0222 
0223 
0224 % add electrodes
0225 try
0226    for i = 1:length(mdl.electrode)
0227       nmdl.electrode(i) = mdl.electrode(i);
0228       nmdl.electrode(i).nodes = find(electrode_node == i);
0229    end
0230 end
0231 
0232 out.uels = uels;
0233 out.els  = els;
0234 out.nn   = nn;
0235 
0236 
0237 function nimg = build_image(nmdl, f2c, img)
0238 nimg = mk_image(nmdl,1);
0239 if isempty(nimg.elem_data) % plane doesn't cut model
0240     return
0241 end
0242 nimg.elem_data = (get_img_data(img)' * f2c)';
0243 try
0244    nimg.calc_colours = img.calc_colours;
0245 end
0246 
0247 
0248 function out = draw_patch(in, nmdl, elem_data, varargin)
0249 
0250 uels = in.uels;
0251 els  = in.els;
0252 nn   = in.nn;
0253 nodes = nmdl.nodes;
0254 
0255 img.elem_data = elem_data;
0256 
0257 out.Vertices = nodes;
0258 out.Faces    = NaN(length(uels),4);
0259 ed = zeros(length(uels),1);
0260 
0261 % show_fem(nimg);
0262 for i = 1:length(uels)
0263     idx = els == uels(i);
0264     id = find(idx);
0265     nnn = nodes(nn(idx),:);
0266     if nnz(idx)==4
0267         if intersection_test(nnn(1,:),nnn(4,:),nnn(2,:),nnn(3,:))
0268             id(3:4) = id([4 3]);
0269         elseif intersection_test(nnn(1,:),nnn(2,:),nnn(3,:),nnn(4,:))
0270             id(2:3) = id([3 2]);
0271         end
0272     end
0273 %     patch(nodes(nn(id),1),nodes(nn(id),2),nodes(nn(id),3),img.elem_data(uels(i)));
0274     out.Faces(i,1:length(id)) = nn(id);
0275     ed(i) = img.elem_data(uels(i));
0276 end
0277 [out.FaceVertexCData scl_data] = calc_colours(ed,varargin{:});
0278 try
0279    out.FaceVertexAlphaData = double(abs(scl_data) > varargin{1}.calc_colours.transparency_thresh);
0280    out.FaceAlpha = 'flat';
0281 end
0282 out.FaceColor = 'flat';
0283 out.CDataMapping = 'direct';
0284 % colormap(calc_colours('colourmap'));
0285 
0286 
0287 function res = intersection_test(A,B,C,D)
0288 % checks for intersection of segments AB and CD
0289 % if AB and CD intersect in 3D, then their projections on a 2D plane also
0290 % intersect (or are colliniar).
0291 
0292 % assume they don't
0293 res = false;
0294 
0295 % check for interesection on the 3 cartesian planes
0296 idx = 1:3;
0297 for i = 0:2
0298     id = circshift(idx',i)';
0299     id = id(1:2);
0300     res = res || ( sign(signed_area(A(id), B(id), C(id))) ~= ...
0301                    sign(signed_area(A(id), B(id), D(id)))        );
0302 end
0303 
0304 function a = signed_area(A,B,C)
0305     a = ( B(1) - A(1) ) * ( C(2) - A(2) ) - ...
0306         ( C(1) - A(1) ) * ( B(2) - A(2) );
0307 
0308 
0309 function [nodeval dist] = nodes_above_or_below(mdl,level)
0310 
0311 % Set a model-dependent tolerance
0312 tol = min(max(mdl.nodes) - min(mdl.nodes)) * 1e-10;
0313 
0314 dist = mdl.nodes(:,3) - level;
0315 dist(abs(dist) < tol) = 0;
0316 nodeval = sign(dist);
0317 
0318 
0319 
0320 
0321 % Level model: usage
0322 %   NODE= level_model( fwd_model, level );
0323 %
0324 % Level is a 1x3 vector specifying the x,y,z axis intercepts
0325 % NODE describes the vertices in this coord space
0326 
0327 function NODE= level_model( fwd_model, level )
0328 
0329    vtx= fwd_model.nodes;
0330    [nn, dims] = size(vtx);
0331    if dims ==2 % 2D case
0332        NODE= vtx';
0333        return;
0334    end
0335 
0336    % Infinities tend to cause issues -> replace with realmax
0337    % Don't need to worry about the sign of the inf
0338    level( isinf(level) | isnan(level) ) = realmax;
0339    level( level==0 ) =     1e-10; %eps;
0340 
0341    % Step 1: Choose a centre point in the plane
0342    %  Weight the point by it's inv axis coords
0343    invlev= 1./level;
0344    ctr= invlev / sum( invlev.^2 );
0345 
0346    % Step 2: Choose basis vectors in the plane
0347    %  First is the axis furthest from ctr
0348    [jnk, s_ax]= sort( - abs(level - ctr) );
0349    v1= [0,0,0]; v1(s_ax(1))= level(s_ax(1));
0350    v1= v1 - ctr;
0351    v1= v1 / norm(v1);
0352 
0353    % Step 3: Get off-plane vector, by cross product
0354    v2= [0,0,0]; v2(s_ax(2))= level(s_ax(2));
0355    v2= v2 - ctr;
0356    v2= v2 / norm(v2);
0357    v3= cross(v1,v2);
0358 
0359    % Step 4: Get orthonormal basis. Replace v2
0360    v2= cross(v1,v3);
0361 
0362    % Step 5: Get bases to point in 'positive directions'
0363    v1= v1 * (1-2*(sum(v1)<0));
0364    v2= v2 * (1-2*(sum(v2)<0));
0365    v3= v3 * (1-2*(sum(v3)<0));
0366 
0367    NODE= [v1;v2;v3] * (vtx' - ctr'*ones(1,nn) );
0368 
0369    
0370 function do_unit_test
0371     imdl = mk_common_model('n3r2',[16,2]);
0372     img = mk_image(imdl.fwd_model,1);
0373     load datacom.mat A B;
0374     img.elem_data(A) = 1.2;
0375     img.elem_data(B) = 0.8;
0376     subplot(131)
0377     show_fem(img);
0378     subplot(132)
0379     cla
0380 %     show_fem(img.fwd_model);
0381     hold on
0382 %     slc = mdl_slice_mesher(img, [3 -3 2]);
0383 %     slc.calc_colours.transparency_thresh = -1;
0384 %     show_fem(slc);
0385     slc = mdl_slice_mesher(img, [0 inf inf]);
0386     slc.calc_colours.transparency_thresh = -1;
0387     slc.fwd_model.boundary = slc.fwd_model.elems;
0388     show_fem(slc);
0389     slc = mdl_slice_mesher(img, [inf inf 2]);
0390     slc.calc_colours.transparency_thresh = -1;
0391     slc.fwd_model.boundary = slc.fwd_model.elems;
0392     show_fem(slc,[0 1 0]);
0393     slc = mdl_slice_mesher(img, [inf inf 2.5]);
0394     slc.calc_colours.transparency_thresh = -1;
0395     slc.fwd_model.boundary = slc.fwd_model.elems;
0396     show_fem(slc);
0397     slc = mdl_slice_mesher(img, [inf inf 1.3]);
0398     slc.calc_colours.transparency_thresh = -1;
0399     slc.fwd_model.boundary = slc.fwd_model.elems;
0400     show_fem(slc);
0401     zlim([0 3]);
0402     view(3)
0403     hold off
0404     subplot(133)
0405     hold on
0406     mdl_slice_mesher(img, [0 inf inf]);
0407     mdl_slice_mesher(img, [inf inf 2]);
0408     mdl_slice_mesher(img, [inf inf 2.5]);
0409     mdl_slice_mesher(img, [inf inf 1.3]);
0410     axis equal
0411     axis tight
0412     view(3)
0413     
0414     % test multi-column image
0415     img.elem_data(:,2) = img.elem_data;
0416     slc = mdl_slice_mesher(img, [0 inf inf]);