Transformation Quad_XYZ --> Quad_XY + Find the radius of a circle tangent to three edges of a Quad

Numeric/Numeric.cpp

@@ -321,6 +321,7 @@ void circumCenterXY(double *p1, double *p2, double *p3, double *res)
   res[1] = (double)((a1 * (x2 - x3) + a2 * (x3 - x1) + a3 * (x1 - x2)) / d);
 }
 
+
 void circumCenterXYZ(double *p1, double *p2, double *p3, double *res, double *uv)
 {
   double v1[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
@@ -348,6 +349,81 @@ void circumCenterXYZ(double *p1, double *p2, double *p3, double *res, double *uv
   res[2] = p1[2] + resP[0] * vx[2] + resP[1] * vy[2];
 }
 
+void planarQuad_xyz2xy(double *x, double *y, double *z, double *xn, double *yn)
+{
+	double v1[3] = {x[1] - x[0], y[1] - y[0], z[1] - z[0]};
+	double v2[3] = {x[2] - x[0], y[2] - y[0], z[2] - z[0]};
+	double v3[3] = {x[3] - x[0], y[3] - y[0], z[3] - z[0]};
+	
+	double vx[3] = {x[1] - x[0], y[1] - y[0], z[1] - z[0]};
+	double vy[3] = {x[2] - x[0], y[2] - y[0], z[2] - z[0]};
+	double vz[3]; prodve(vx, vy, vz); prodve(vz, vx, vy);
+	
+	norme(vx); norme(vy); norme(vz);
+	
+	double p1P[2] = {0.0, 0.0};
+	double p2P[2]; prosca(v1, vx, &p2P[0]); prosca(v1, vy, &p2P[1]);
+	double p3P[2]; prosca(v2, vx, &p3P[0]); prosca(v2, vy, &p3P[1]);
+	double p4P[2]; prosca(v3, vx, &p4P[0]); prosca(v3, vy, &p4P[1]);
+	
+	xn[0] = p1P[0];
+	xn[1] = p2P[0];
+	xn[2] = p3P[0];
+	xn[3] = p4P[0];
+	yn[0] = p1P[1];
+	yn[1] = p2P[1];
+	yn[2] = p3P[1];
+	yn[3] = p4P[1];
+}
+
+double computeInnerRadiusForQuad(double *x, double *y, int i){
+	
+	// parameters of the equations of the 3 edges	
+	double a1 = y[(4+i)%4]-y[(5+i)%4];
+	double a2 = y[(5+i)%4]-y[(6+i)%4];
+	double a3 = y[(6+i)%4]-y[(7+i)%4];
+
+	double b1 = x[(5+i)%4]-x[(4+i)%4];
+	double b2 = x[(6+i)%4]-x[(5+i)%4];
+	double b3 = x[(7+i)%4]-x[(6+i)%4];
+
+	double c1 = y[(5+i)%4]*x[(4+i)%4]-y[(4+i)%4]*x[(5+i)%4];
+	double c2 = y[(6+i)%4]*x[(5+i)%4]-y[(5+i)%4]*x[(6+i)%4];
+	double c3 = y[(7+i)%4]*x[(6+i)%4]-y[(6+i)%4]*x[(7+i)%4];
+
+	// length of the 3 edges
+	double l1 = sqrt(a1*a1+b1*b1);
+	double l2 = sqrt(a2*a2+b2*b2);
+	double l3 = sqrt(a3*a3+b3*b3);
+	
+	// parameters of the 2 bisectors
+	double a12 = a1/l1-a2/l2;
+	double a23 = a2/l2-a3/l3;
+	
+	double b12 = b1/l1-b2/l2;
+	double b23 = b2/l2-b3/l3;
+	
+	double c12 = c1/l1-c2/l2;
+	double c23 = c2/l2-c3/l3;
+	
+	// compute the coordinates of the center of the incircle, 
+	// that is the point where the 2 bisectors meet
+	double x_s = (c12*b23-c23*b12)/(a23*b12-a12*b23);
+	double y_s = 0.;
+	if (b12 != 0) {
+		y_s = -a12/b12*x_s-c12/b12;
+	}
+	else {
+		y_s = -a23/b23*x_s-c23/b23;
+	}
+ 
+	// finally get the radius of the circle 
+	double r = (a1*x_s+b1*y_s+c1)/l1;
+
+	return r;
+}
+
+
 char float2char(float f)
 {
   // float normalized in [-1, 1], char in [-127, 127]

Numeric/Numeric.h

@@ -66,6 +66,11 @@ double triangle_area2d(double p0[2], double p1[2], double p2[2]);
 double triangle_polar_inertia(double p0[2], double p1[2], double p2[2]);
 void circumCenterXY(double *p1, double *p2, double *p3, double *res);
 void circumCenterXYZ(double *p1, double *p2, double *p3, double *res, double *uv=0);
+// operate a transformation on the 4 points of a Quad in 3D, to have an equivalent Quad in 2D
+void planarQuad_xyz2xy(double *x, double *y, double *z, double *xn, double *yn);
+// compute the radius of the circle that is tangent to the 3 edges defined by 4 points
+// edge_1=(x0,y0)->(x1,y1); edge_2=(x1,y1)->(x2,y2); edge_3=(x2,y2)->(x3,y3) 
+double computeInnerRadiusForQuad(double *x, double *y, int i);
 char float2char(float f);
 float char2float(char c);
 void eigenvalue(double mat[3][3], double re[3]);