Home > eidors > graphics > matlab > mdl_slice_mesher.m

mdl_slice_mesher

PURPOSE ^

MDL_SLICE_MESHER A slice of a 3D FEM as a 2D FEM

SYNOPSIS ^

function [nimg out] = mdl_slice_mesher(fmdl,level,varargin)

DESCRIPTION ^

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 5734 2018-03-29 23:57:07Z 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    end
0113 end
0114 e_edges = (sum(e_nodes(edges),2)/ 2)  .* (e_nodes(edges(:,1)) == e_nodes(edges(:,2)));
0115 %% crossed edges
0116 % exclude edges on plane (dealt with later)
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 % new nodes along the edges
0121 nodes = mdl.nodes(edges(idx,1),:) + ...
0122     repmat(t,1,3).*(mdl.nodes(edges(idx,2),:) - mdl.nodes(edges(idx,1),:));
0123 % nn indexes the just-created nodes, els indexes elements
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 %% crossed nodes
0133 idx = nodeval == 0;
0134 ln = length(nodes); %store the size
0135 nodes = [nodes; mdl.nodes(idx,:)]; % add the crossed nodes to the new model
0136 electrode_node = [electrode_node; e_nodes(idx)];
0137 % nnn indexes mdl.nodes(idx), eee indexes mdl.elems
0138 [nnn eee] = find(mdl.node2elem(idx,:));
0139 nnn = nnn + ln; %make a proper index into nodes
0140 % n1: eee = ueee(n1)
0141 [ueee jnk n1] = unique(eee,'last');
0142 nodes_per_elem = jnk;
0143 nodes_per_elem(2:end) = diff(jnk);
0144 % if an elem has 3 crossed nodes, there must be 2 of them, add both for now
0145 add = find(nodes_per_elem == 3);
0146 if ~isempty(add)
0147     % what to do with faces shared between elements?
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 % for elems with less than 4 crossed edges -> add crossed nodes if needed
0158 [uels jnk n] = unique(els_edge,'last');
0159 
0160 % only consider elements who have both a crossed node and edge
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 % TODO: Speed this up
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 % deal with double elements (from shared faces)
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 % put all source elements that contribute to destination element on one
0210 % column
0211 f2c = sparse(x,idx,1);
0212 % ensure correct number of columns (happens when the last source element
0213 % doesn't contribute
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 % add electrodes
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) % plane doesn't cut model
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 % show_fem(nimg);
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 %     patch(nodes(nn(id),1),nodes(nn(id),2),nodes(nn(id),3),img.elem_data(uels(i)));
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 % colormap(calc_colours('colourmap'));
0282 
0283 
0284 function res = intersection_test(A,B,C,D)
0285 % checks for intersection of segments AB and CD
0286 % if AB and CD intersect in 3D, then their projections on a 2D plane also
0287 % intersect (or are colliniar).
0288 
0289 % assume they don't
0290 res = false;
0291 
0292 % check for interesection on the 3 cartesian planes
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 % Set a model-dependent tolerance
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 % Level model: usage
0319 %   NODE= level_model( fwd_model, level );
0320 %
0321 % Level is a 1x3 vector specifying the x,y,z axis intercepts
0322 % NODE describes the vertices in this coord space
0323 
0324 function NODE= level_model( fwd_model, level )
0325 
0326    vtx= fwd_model.nodes;
0327    [nn, dims] = size(vtx);
0328    if dims ==2 % 2D case
0329        NODE= vtx';
0330        return;
0331    end
0332 
0333    % Infinities tend to cause issues -> replace with realmax
0334    % Don't need to worry about the sign of the inf
0335    level( isinf(level) | isnan(level) ) = realmax;
0336    level( level==0 ) =     1e-10; %eps;
0337 
0338    % Step 1: Choose a centre point in the plane
0339    %  Weight the point by it's inv axis coords
0340    invlev= 1./level;
0341    ctr= invlev / sum( invlev.^2 );
0342 
0343    % Step 2: Choose basis vectors in the plane
0344    %  First is the axis furthest from ctr
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    % Step 3: Get off-plane vector, by cross product
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    % Step 4: Get orthonormal basis. Replace v2
0357    v2= cross(v1,v3);
0358 
0359    % Step 5: Get bases to point in 'positive directions'
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 %     show_fem(img.fwd_model);
0378     hold on
0379 %     slc = mdl_slice_mesher(img, [3 -3 2]);
0380 %     slc.calc_colours.transparency_thresh = -1;
0381 %     show_fem(slc);
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     % test multi-column image
0412     img.elem_data(:,2) = img.elem_data;
0413     slc = mdl_slice_mesher(img, [0 inf inf]);

Generated on Fri 01-Jun-2018 15:59:55 by m2html © 2005