[cutpos,cutvalue,facedata,elemid]=qmeshcut(elem,node,value,cutat)
fast tetrahedral mesh slicer
author:Qianqian Fang, <fangq at nmr.mgh.harvard.edu>
input:
elem: integer array with dimensions of NE x 4, each row contains
the indices of all the nodes for each tetrahedron
node: node coordinates, 3 columns for x, y and z respectively
value: a scalar array with the length of node numbers, can have
multiple columns
cutat: cutat can have different forms:
if cutat is a 3x3 matrix, it defines a plane by 3 points:
cutat=[x1 y1 z1;x2 y2 z2;x3 y3 z3]
if cutat is a vector of 4 element, it defines a plane by
a*x+b*y+c*z+d=0 and cutat=[a b c d]
if cutat is a single scalar, it defines an isosurface
inside the mesh at value=cutat
if cutat is a string, it defines an implicit surface
at which the cut is made. it must has form expr1=expr2
where expr1 expr2 are expressions made of x,y,z,v and
constants
if cutat is a cell in the form of {'expression',scalar},
the expression will be evaluated at each node to
produce a new value, then an isosurface is produced
at the levelset where new value=scalar; the
expression can contain constants and x,y,z,v
output:
cutpos: all the intersections of mesh edges by the cutat
cutpos is similar to node, containing 3 columns for x/y/z
cutvalue: interpolated values at the intersections, with row number
being the num. of the intersections, column number being the
same as "value".
facedata: define the intersection polygons in the form of patch "Faces"
elemid: the index of the elem in which each intersection polygon locates
without any output, qmeshcut generates a cross-section plot
the outputs of this subroutine can be easily plotted using
% if value has a length of node:
patch('Vertices',cutpos,'Faces',facedata,'FaceVertexCData',cutvalue,'FaceColor','interp');
% if value has a length of elem:
patch('Vertices',cutpos,'Faces',facedata,'CData',cutvalue,'FaceColor','flat');
-- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)0001 function [cutpos,cutvalue,facedata,elemid]=qmeshcut(elem,node,value,cutat,varargin) 0002 % 0003 % [cutpos,cutvalue,facedata,elemid]=qmeshcut(elem,node,value,cutat) 0004 % 0005 % fast tetrahedral mesh slicer 0006 % 0007 % author:Qianqian Fang, <fangq at nmr.mgh.harvard.edu> 0008 % 0009 % input: 0010 % elem: integer array with dimensions of NE x 4, each row contains 0011 % the indices of all the nodes for each tetrahedron 0012 % node: node coordinates, 3 columns for x, y and z respectively 0013 % value: a scalar array with the length of node numbers, can have 0014 % multiple columns 0015 % cutat: cutat can have different forms: 0016 % if cutat is a 3x3 matrix, it defines a plane by 3 points: 0017 % cutat=[x1 y1 z1;x2 y2 z2;x3 y3 z3] 0018 % if cutat is a vector of 4 element, it defines a plane by 0019 % a*x+b*y+c*z+d=0 and cutat=[a b c d] 0020 % if cutat is a single scalar, it defines an isosurface 0021 % inside the mesh at value=cutat 0022 % if cutat is a string, it defines an implicit surface 0023 % at which the cut is made. it must has form expr1=expr2 0024 % where expr1 expr2 are expressions made of x,y,z,v and 0025 % constants 0026 % if cutat is a cell in the form of {'expression',scalar}, 0027 % the expression will be evaluated at each node to 0028 % produce a new value, then an isosurface is produced 0029 % at the levelset where new value=scalar; the 0030 % expression can contain constants and x,y,z,v 0031 % 0032 % output: 0033 % cutpos: all the intersections of mesh edges by the cutat 0034 % cutpos is similar to node, containing 3 columns for x/y/z 0035 % cutvalue: interpolated values at the intersections, with row number 0036 % being the num. of the intersections, column number being the 0037 % same as "value". 0038 % facedata: define the intersection polygons in the form of patch "Faces" 0039 % elemid: the index of the elem in which each intersection polygon locates 0040 % 0041 % without any output, qmeshcut generates a cross-section plot 0042 % 0043 % the outputs of this subroutine can be easily plotted using 0044 % 0045 % % if value has a length of node: 0046 % patch('Vertices',cutpos,'Faces',facedata,'FaceVertexCData',cutvalue,'FaceColor','interp'); 0047 % 0048 % % if value has a length of elem: 0049 % patch('Vertices',cutpos,'Faces',facedata,'CData',cutvalue,'FaceColor','flat'); 0050 % 0051 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net) 0052 % 0053 0054 % get the coefficients of the cutat equation: ax+by+cz+d=0 0055 if(nargin<4) 0056 error('qmeshcut requires at least 4 inputs'); 0057 end 0058 if(size(value,1)~=size(node,1) && size(value,1)~=size(elem,1) && ~isempty(value)) 0059 error('the length of value must be either that of node or elem'); 0060 end 0061 if(isempty(value)) 0062 cutvalue=[]; 0063 end 0064 if(ischar(cutat) || (iscell(cutat) && length(cutat)==2 && ischar(cutat{1}))) 0065 x=node(:,1); 0066 y=node(:,2); 0067 z=node(:,3); 0068 if(ischar(cutat)) 0069 expr=regexp(cutat,'(.+)=(.+)','tokens','once'); %regexp(cutat,'=','split'); 0070 if(length(expr)~=2) error('single expression must contain a single "=" sign'); end 0071 dist=eval(expr{1})-eval(expr{2}); 0072 else 0073 dist=eval(cutat{1})-cutat{2}; 0074 end 0075 if(all(dist<=0)) 0076 asign=-double(dist<0); 0077 asign(asign==0)=1; 0078 else 0079 asign=double(dist>0); 0080 asign(asign==0)=-1; 0081 end 0082 elseif(numel(cutat)==9 || numel(cutat)==4) 0083 if(numel(cutat)==9) 0084 [a,b,c,d]=getplanefrom3pt(cutat); 0085 else 0086 [a,b,c,d]=deal(cutat(:)); 0087 end 0088 0089 % compute which side of the cutat for all nodes in the mesh 0090 co=repmat([a b c],size(node,1),1); 0091 dist=sum( (co.*node)' )+d; 0092 asign=dist; 0093 asign(find(asign>=0))=1; 0094 asign(find(asign<0))=-1; 0095 else 0096 if(size(value,1)~=size(node,1)) 0097 error('must use nodal value list when cutting mesh at an isovalue'); 0098 end 0099 dist=value-cutat; 0100 if(all(dist<=0)) 0101 asign=-double(dist<0); 0102 asign(asign==0)=1; 0103 else 0104 asign=double(dist>0); 0105 asign(asign==0)=-1; 0106 end 0107 end 0108 0109 % get all the edges of the mesh 0110 esize=size(elem,2); 0111 if(esize==4) 0112 edges=[elem(:,[1,2]);elem(:,[1,3]);elem(:,[1,4]); 0113 elem(:,[2,3]);elem(:,[2,4]);elem(:,[3,4])]; 0114 elseif(esize==3) 0115 edges=[elem(:,[1,2]);elem(:,[1,3]);elem(:,[2,3])]; 0116 elseif(esize==10) 0117 edges=[elem(:,[1,5]);elem(:,[1,8]);elem(:,[1,7]); 0118 elem(:,[2,5]);elem(:,[2,6]);elem(:,[2,9]); 0119 elem(:,[3,6]);elem(:,[3,7]);elem(:,[3,10]); 0120 elem(:,[4,8]);elem(:,[4,9]);elem(:,[4,10])]; 0121 end 0122 0123 % find all edges with two ends at the both sides of the plane 0124 edgemask=sum(asign(edges),2); 0125 cutedges=find(edgemask==0); 0126 %edgemask=prod(asign(edges)'); 0127 %cutedges=find(edgemask<0); 0128 0129 % calculate the distances of the two nodes, and use them as interpolation weight 0130 cutweight=dist(edges(cutedges,:)); 0131 totalweight=diff(cutweight'); 0132 0133 %caveat: if an edge is co-planar to the cutat, then totalweight will be 0 0134 % and dividing zero will cause trouble for cutweight 0135 0136 cutweight=abs(cutweight./repmat(totalweight(:),1,2)); 0137 0138 % calculate the cross-cut position and the interpolated values 0139 0140 cutpos=node(edges(cutedges,1),:).*repmat(cutweight(:,2),[1 3])+... 0141 node(edges(cutedges,2),:).*repmat(cutweight(:,1),[1 3]); 0142 if(size(value,1)==size(node,1)) 0143 if(iscell(cutat) || ischar(cutat) || numel(cutat)==9 || numel(cutat)==4) 0144 cutvalue=value(edges(cutedges,1),:).*repmat(cutweight(:,2),[1 size(value,2)])+... 0145 value(edges(cutedges,2),:).*repmat(cutweight(:,1),[1 size(value,2)]); 0146 elseif(numel(cutat)==1) 0147 cutvalue=ones(size(cutpos,1),1)*cutat; 0148 end 0149 end 0150 % organize all cross-cuts into patch facedata format 0151 0152 emap=zeros(size(edges,1),1); 0153 emap(cutedges)=1:length(cutedges); 0154 if(esize==10) 0155 emap=reshape(emap,[size(elem,1),12]); % 10-node element 0156 else 0157 emap=reshape(emap,[size(elem,1),esize*(esize-1)/2]); % C^n_2 0158 end 0159 faceid=find(any(emap,2)==1); 0160 facelen=length(faceid); 0161 0162 % cross-cuts can only be triangles or quadrilaterals for tetrahedral mesh 0163 % (co-plannar mesh needs to be considered) 0164 0165 etag=sum(emap>0,2); % emap & etag are of length size(elem,1) 0166 0167 if(esize==3) % surface mesh cut by a plane 0168 linecut=find(etag==2); 0169 lineseg=emap(linecut,:)'; 0170 facedata=reshape(lineseg(find(lineseg)),[2,length(linecut)])'; 0171 elemid=linecut; 0172 if(size(value,1)==size(elem,1) && ~exist('cutvalue','var')) 0173 cutvalue=value(elemid,:); 0174 end 0175 return; 0176 end 0177 0178 tricut=find(etag==3); 0179 quadcut=find(etag==4); 0180 elemid=[tricut(:);quadcut(:)]; 0181 if(size(value,1)==size(elem,1) && ~exist('cutvalue','var')) 0182 cutvalue=value(elemid,:); 0183 end 0184 % fast way (vector-form) to get all triangles 0185 0186 tripatch=emap(tricut,:)'; 0187 tripatch=reshape(tripatch(find(tripatch)),[3,length(tricut)])'; 0188 0189 % fast way to get all quadrilaterals in convexhull form (avoid using 0190 % convhulln) 0191 0192 quadpatch=emap(quadcut,:)'; 0193 quadpatch=reshape(quadpatch(find(quadpatch)),[4,length(quadpatch)])'; 0194 0195 % combine the two sets to create the final facedata 0196 % using the matching-tetrahedra algorithm as shown in 0197 % https://visualization.hpc.mil/wiki/Marching_Tetrahedra 0198 0199 facedata=[tripatch(:,[1 2 3 3]); quadpatch(:,[1 2 4 3])]; 0200 0201 % plot your results with the following command 0202 0203 if(nargout==0) 0204 patch('Vertices',cutpos,'Faces',facedata,'FaceVertexCData',cutvalue,'facecolor','interp',varargin{:}); 0205 end