From 3e85ffd0ff6065666f1b06981f3d59ed69cdc1c8 Mon Sep 17 00:00:00 2001
From: cignoni <paolo.cignoni@isti.cnr.it>
Date: Thu, 8 Jan 2009 11:24:54 +0000
Subject: [PATCH] New function for computing baric coords in a more robust way

---
 vcg/space/triangle3.h | 55 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 55 insertions(+)

diff --git a/vcg/space/triangle3.h b/vcg/space/triangle3.h
index 2b4ae6d2..e8547df9 100644
--- a/vcg/space/triangle3.h
+++ b/vcg/space/triangle3.h
@@ -94,6 +94,7 @@ Initial commit
 #define __VCG_TRIANGLE3
 
 #include <vcg/space/box3.h>
+#include <vcg/space/point2.h>
 #include <vcg/space/point3.h>
 #include <vcg/space/plane3.h>
 #include <vcg/space/segment3.h>
@@ -163,6 +164,60 @@ typename TriangleType::ScalarType QualityFace(const TriangleType &t)
   return Quality(t.cP(0), t.cP(1), t.cP(2));
 }
 
+// More robust function to computing barycentric coords of a point inside a triangle.
+// it requires the knowledge of what is the direction that is more orthogonal to the face plane. 
+// Usually this info can be stored in a bit of the face flags (see updateFlags class members)
+// This direction is used to project the triangle in 2D and solve the problem in 2D where it is well defined.
+
+template<class TriangleType, class ScalarType>
+bool InterpolationParameters(const TriangleType t, const int Axis, const Point3<ScalarType> & P,  Point3<ScalarType> & L)
+{
+  Point2<ScalarType> test;
+	typedef Point2<ScalarType> P2;
+	if(Axis==0) return InterpolationParameters2( P2(t.P(0)[1],t.P(0)[2]), P2(t.P(1)[1],t.P(1)[2]), P2(t.P(2)[1],t.P(2)[2]), P2(P[1],P[2]), L);
+	if(Axis==1) return InterpolationParameters2( P2(t.P(0)[0],t.P(0)[2]), P2(t.P(1)[0],t.P(1)[2]), P2(t.P(2)[0],t.P(2)[2]), P2(P[0],P[2]), L);
+	if(Axis==2) return InterpolationParameters2( P2(t.P(0)[0],t.P(0)[1]), P2(t.P(1)[0],t.P(1)[1]), P2(t.P(2)[0],t.P(2)[1]), P2(P[0],P[1]), L);
+	return false; 
+}
+
+
+// Function that computes the barycentric coords of a 2D triangle. Used by the above function. 
+// Algorithm: simply find a base for the frame of the triangle, assuming v3 as origin (matrix T) invert it and apply to P-v3.
+
+template<class ScalarType>
+bool InterpolationParameters2(const Point2<ScalarType> &V1,
+													   const Point2<ScalarType> &V2,
+														 const Point2<ScalarType> &V3,
+														 const Point2<ScalarType> &P, Point3<ScalarType> &L)
+{
+	ScalarType T00 = V1[0]-V3[0];  ScalarType T01 = V2[0]-V3[0];  
+	ScalarType T10 = V1[1]-V3[1];  ScalarType T11 = V2[1]-V3[1]; 
+	ScalarType Det = T00 * T11 - T01*T10;
+	if(fabs(Det) < 0.0000001) 
+				return false;
+	
+	ScalarType IT00 =  T11/Det;	ScalarType IT01 = -T01/Det;	
+	ScalarType IT10 = -T10/Det;	ScalarType IT11 =  T00/Det;	
+
+	Point2<ScalarType> Delta = P-V3;
+	
+	L[0] = IT00*Delta[0] + IT01*Delta[1];
+	L[1] = IT10*Delta[0] + IT11*Delta[1]; 
+ 
+	if(L[0]<0) L[0]=0;
+	if(L[1]<0) L[1]=0;
+	if(L[0]>1.) L[0]=1;
+	if(L[1]>1.) L[1]=1;
+	
+	L[2] = 1. - L[1] - L[0];
+	if(L[2]<0) L[2]=0;
+	
+	assert(L[2] >= -0.00001);
+
+	return true;
+}
+
+
 /** Calcola i coefficienti della combinazione convessa.
 	@param bq Punto appartenente alla faccia
 	@param a Valore di ritorno per il vertice V(0)