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

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005