merge metch toolbox with iso2mesh

Author: fangq

affinemap.m

@@ -0,0 +1,40 @@
+function [A,b]=affinemap(pfrom,pto)
+%  [A,b]=affinemap(pfrom,pto)
+%
+%  calculate an affine transform (A matrix and b vector) to map n
+%  vertices from one space to the other using least square solutions
+%
+%  author: Qianqian Fang <q.fang at neu.edu>
+%  date: 12/12/2008
+%
+% parameters: 
+%      pfrom: nx3 matrix, each row is a 3d point in original space
+%      pto: nx3 matrix, each row is a 3d point in the mapped space
+%
+% outputs:
+%      A: 3x3 matrix, the calculated affine A matrix
+%      b: 3x1 vector, the calculated affine b vector
+%
+% the solution will satisfy the following equation: A*pfrom'+b=pto
+%
+% Please find more information at http://iso2mesh.sf.net/cgi-bin/index.cgi?metch
+%
+% this function is part of "metch" toobox, see COPYING for license
+
+bsubmat=eye(3);
+ptnum=size(pfrom,1);
+if(size(pto,1)~=ptnum)
+        error('two inputs should have the same size');
+end
+amat=zeros(ptnum*3,9);
+for i=1:ptnum
+        amat(i*3-2:i*3,:)=kron(bsubmat,pfrom(i,:));
+end
+amat=[amat,repmat(bsubmat,ptnum,1)];
+
+bvec=pto';
+bvec=bvec(:);
+
+x=amat\bvec;
+A=reshape(x(1:9),3,3)';
+b=x(end-2:end);

dist2surf.m

@@ -0,0 +1,37 @@
+function [d2surf,cn]=dist2surf(node,nv,p,cn)
+%  [d2surf,cn]=dist2surf(node,nv,p)
+%
+%  calculate the distances from a point cloud to a surface, and return
+%  the indices of the closest surface node
+%
+%  author: Qianqian Fang <q.fang at neu.edu>
+%  date: 12/12/2008
+%
+% parameters: 
+%      node: node coordinate of the surface mesh (nn x 3)
+%      nv: nodal norms (vector) calculated from nodesurfnorm.m
+%          with dimensions of (size(node,1),3), this can be 
+%          calcuated from nodesurfnorm.m
+%      pt: points to be calculated, 3 columns for x,y and z respectively
+%
+% outputs:
+%      d2surf: a vector of length of p, the distances from p(i) to the surface
+%      cn: a integer vector with the length of p, the indices of the closest surface node
+%
+% Please find more information at http://iso2mesh.sf.net/cgi-bin/index.cgi?metch
+%
+% this function is part of "metch" toobox, see COPYING for license
+
+
+if(nargin<4)
+        nn=size(node,1);
+        pnum=size(p,1);
+        mindist=zeros(pnum,1);
+        cn=zeros(pnum,1);
+        for i=1:pnum
+                d0=node-repmat(p(i,:),nn,1);
+                d0=sum(d0.*d0,2);
+                [mindist(i),cn(i)]=min(d0);
+        end
+end
+d2surf=abs(sum(nv(cn,:).*(p-node(cn,:)),2));

finddisconnsurf.m

@@ -5,7 +5,7 @@
 % subroutine to extract disconnected surfaces from a 
 % cluster of surfaces
 % 
-% author: Qianqian Fang ([email protected].edu)
+% author: Qianqian Fang (q.fang at neu.edu)
 % Date: 2008/03/06
 %
 % input: 

linextriangle.m

@@ -0,0 +1,74 @@
+function [isinside,pt,coord]=linextriangle(p0,p1,plane)
+%  [isinside,pt,coord]=linextriangle(p0,p1,plane)
+%
+%  calculate the intersection of a 3d line (passing two points)
+%  with a plane (determined by 3 points)
+%
+%  author: Qianqian Fang <q.fang at neu.edu>
+%  date: 12/12/2008
+%
+% parameters: 
+%      p0: a 3d point in form of (x,y,z)
+%      p1: another 3d point in form of (x,y,z), p0 and p1 determins the line
+%      plane: a 3x3 matrix, each row is a 3d point in form of (x,y,z)
+%             this is used to define a plane
+% outputs:
+%      isinside: a boolean variable, 1 for the intersection is within the 
+%               3d triangle determined by the 3 points in plane; 0 is outside
+%      pt: the coordinates of the intersection pint
+%      coord: 1x3 vector, if isinside=1, coord will record the barycentric 
+%          coordinate for the intersection point within the triangle; 
+%          otherwise it will be all zeros.
+%
+% for degenerated lines or triangles, this will stop
+%
+% Please find more information at http://iso2mesh.sf.net/cgi-bin/index.cgi?metch
+%
+% this function is part of "metch" toobox, see COPYING for license
+
+[a,b,c,d]=getplanefrom3pt(plane);
+
+if(a*a+b*b+c*c==0.0)
+        error('degenerated plane');
+end
+
+dl_n=sum([a b c].*(p1-p0));
+
+if(dl_n==0.0)
+        error('degenerated line');
+end
+
+% solve for the intersection point
+t=-(a*p0(1)+b*p0(2)+c*p0(3)+d)/dl_n;
+pt=p0+(p1-p0)*t;
+
+
+dist=sum(abs(diff(plane)));
+[md,imax]=sort(dist);
+if(md(2)==0.0)
+        error('degenerated triangle');
+end
+goodidx=imax(2:end);
+
+ptproj=pt(goodidx);
+mat0=[plane(:,goodidx)',ptproj';1 1 1 1];
+
+isinside=0;
+coord=[0 0 0];
+
+det1=det(mat0(:,[4 2 3],:));
+det2=det(mat0(:,[1 4 3],:));
+if(det1*det2<0) 
+        return; 
+end
+det3=det(mat0(:,[1 2 4],:));
+if(det2*det3<0) 
+        return;
+end
+if(det1*det3<0) 
+        return;
+end
+isinside=1;
+det0=det(mat0(:,1:3));
+
+coord=[det1 det2 det3]/det0;