Ajout des fonctions de manipulation de la structure de polygone.

Exportation de la fonction de lissage par des Beziers.
Ajout de la fonction d'interpolation par des Beziers.

1 files changed, 608 insertions, 10 deletions

diff --git a/generic/Geo.c b/generic/Geo.c
index c564be7..aca1590 100644
--- a/generic/Geo.c
+++ b/generic/Geo.c
@@ -35,12 +35,76 @@
 #include "WidgetInfo.h"
 
 #include <memory.h>
+#include <malloc.h>
 
 
 static const char rcsid[] = "$Id$";
 static const char compile_id[]="$Compile: " __FILE__ " " __DATE__ " " __TIME__ " $";
 
 
+void
+POLY_INIT(ZnPoly	*poly) {
+  poly->num_contours = 0;
+  poly->contours = NULL;
+  poly->holes = NULL;
+  poly->cw = NULL;
+}
+
+void
+POLY_CONTOUR1(ZnPoly	*poly,
+	      ZnPoint	*pts,
+	      int	num_pts) {
+  poly->num_contours = 1;
+  poly->holes = &poly->hole1;
+  poly->hole1 = False;
+  poly->cw = &poly->cw1;
+  poly->contours = &poly->contour1;
+  poly->contour1.num_points = num_pts;
+  poly->contour1.points = pts;
+}
+
+void
+POLY_SET(ZnPoly	*poly1,
+	 ZnPoly	*poly2) {
+  POLY_FREE(poly1);
+  if (poly2->num_contours == 1) {
+    POLY_CONTOUR1(poly1, poly2->contours[0].points, poly2->contours[0].num_points);
+    if (poly2->contours != &poly2->contour1) {
+      ZnFree(poly2->contours);
+    }
+  }
+  else {
+    poly1->num_contours = poly2->num_contours;
+    poly1->contours = poly2->contours;
+    poly1->holes = poly2->holes;
+    poly1->cw = poly2->cw;
+  }
+}
+
+void
+POLY_FREE(ZnPoly	*poly) {
+  if (poly->num_contours) {
+    int i;
+    for (i = 0; i < poly->num_contours; i++) {
+      ZnFree(poly->contours[i].points);
+    }
+    if ((poly->contours != &poly->contour1) &&
+        (poly->num_contours > 1)) {
+      ZnFree(poly->contours);
+      if (poly->holes) {
+        ZnFree(poly->holes);
+      }
+      if (poly->cw) {
+        ZnFree(poly->cw);
+      }
+    }
+    poly->num_contours = 0;
+    poly->contours = NULL;
+    poly->holes = NULL;
+    poly->cw = NULL;
+  }
+}
+
 /*
  * Compute the origin of the rectangle given
  * by position, anchor, width and height.
@@ -828,7 +892,7 @@ PolylineInBBox(ZnPoint	*points,
        * poly[2].
        */
       if ((join_style == JoinBevel) || do_miter_as_bevel) {
-	if (PolygonInBBox(poly, 4, bbox) != inside) {
+	if (PolygonInBBox(poly, 4, bbox, NULL) != inside) {
 	  return 0;
 	}
 	do_miter_as_bevel = False;
@@ -853,7 +917,7 @@ PolylineInBBox(ZnPoint	*points,
       GetButtPoints(points, &points[1], width, 0, &poly[2], &poly[3]);
     }
 
-    if (PolygonInBBox(poly, 4, bbox) != inside) {
+    if (PolygonInBBox(poly, 4, bbox, NULL) != inside) {
       return 0;
     }
   }
@@ -874,17 +938,23 @@ PolylineInBBox(ZnPoint	*points,
 /*
  * Tell where a polygon is with respect to an area.
  * Return -1 if the polygon is entirely outside the bbox, 1
- * if it is entirely inside and 0 otherwise.
+ * if it is entirely inside and 0 otherwise. If area_enclosed
+ * is not NULL it tells whether the area is enclosed by the
+ * polygon or not.
  */
 int
 PolygonInBBox(ZnPoint	*points,
 	      int	num_points,
-	      ZnBBox	*bbox)
+	      ZnBBox	*bbox,
+	      ZnBool	*area_enclosed)
 {
   int		inside, count;
   ZnPoint	*p, *head, *first, *second;
   ZnBool	closed;
 
+  if (area_enclosed) {
+    *area_enclosed = False;
+  }
   p = head = points;
   closed = True;
   count = num_points-2;
@@ -937,6 +1007,9 @@ PolygonInBBox(ZnPoint	*points,
 	 num_points, bbox->orig.x, bbox->orig.y,
 	 PolygonToPointDist(points, num_points, &bbox->orig));*/
   if (PolygonToPointDist(points, num_points, &bbox->orig) <= 0.0) {
+    if (area_enclosed) {
+      *area_enclosed = True;
+    }
     return 0;
   }
   
@@ -2058,18 +2131,14 @@ GetArcPath(ZnReal	start_angle,
  **********************************************************************************
  */
 void
-SmoothPathWithBezier(ZnList	from_points,
+SmoothPathWithBezier(ZnPoint	*fp,
+		     int	num_fp,
 		     ZnList	to_points)
 {
-  ZnPoint	*fp;
-  int		num_fp;
   ZnBool	closed;
   ZnPoint	s[4];
   int		i;
 
-  fp = (ZnPoint *) ZnListArray(from_points);
-  num_fp = ZnListSize(from_points);
-
   /*
    * make sure the output vector is empty
    */
@@ -2150,6 +2219,535 @@ SmoothPathWithBezier(ZnList	from_points,
 /*
  **********************************************************************************
  *
+ * FitBezier --
+ *  	Fit a Bezier curve to a (sub)set of digitized points.
+ *
+ *	From: An Algorithm for Automatically Fitting Digitized Curves
+ *	      by Philip J. Schneider in "Graphics Gems", Academic Press, 1990
+ *
+ **********************************************************************************
+ */
+
+static ZnReal
+V2DistanceBetween2Points(ZnPoint	*a,
+			 ZnPoint	*b)
+{
+  ZnReal dx = a->x - b->x;
+  ZnReal dy = a->y - b->y;
+  return sqrt((dx*dx)+(dy*dy));
+}
+
+static ZnReal
+V2SquaredLength(ZnPoint	*a)
+{	
+  return (a->x * a->x)+(a->y * a->y);
+}
+
+static ZnReal
+V2Length(ZnPoint	*a)
+{
+  return sqrt(V2SquaredLength(a));
+}
+	
+static ZnPoint *
+V2Scale(ZnPoint	*v,
+	ZnReal	newlen)
+{
+  ZnReal len = V2Length(v);
+  if (len != 0.0) {
+    v->x *= newlen/len;
+    v->y *= newlen/len;
+  }
+  return v;
+}
+
+static ZnPoint *
+V2Negate(ZnPoint *v)
+{
+  v->x = -v->x;  v->y = -v->y;
+  return v;
+}
+
+static ZnPoint *
+V2Normalize(ZnPoint *v)
+{
+  ZnReal len = sqrt(V2Length(v));
+  if (len != 0.0) {
+    v->x /= len;
+    v->y /= len;
+  }
+  return v;
+}
+static ZnPoint *
+V2Add(ZnPoint	*a,
+      ZnPoint	*b,
+      ZnPoint	*c)
+{
+  c->x = a->x + b->x;
+  c->y = a->y + b->y;
+  return c;
+}
+	
+static ZnReal
+V2Dot(ZnPoint	*a,
+      ZnPoint	*b) 
+{
+  return (a->x*b->x) + (a->y*b->y);
+}
+
+static ZnPoint
+V2AddII(ZnPoint	a,
+	ZnPoint	b)
+{
+  ZnPoint c;
+  c.x = a.x + b.x;
+  c.y = a.y + b.y;
+  return c;
+}
+
+static ZnPoint
+V2ScaleIII(ZnPoint	v,
+	   ZnReal	s)
+{
+  ZnPoint result;
+  result.x = v.x * s;
+  result.y = v.y * s;
+  return result;
+}
+
+static ZnPoint
+V2SubII(ZnPoint	a,
+	ZnPoint	b)
+{
+  ZnPoint c;
+  c.x = a.x - b.x;
+  c.y = a.y - b.y;
+  return c;
+}
+
+/*
+ * B0, B1, B2, B3, Bezier multipliers.
+ */
+static ZnReal
+B0(ZnReal	u)
+{
+  ZnReal tmp = 1.0 - u;
+  return tmp * tmp * tmp;
+}
+
+static ZnReal
+B1(ZnReal	u)
+{
+  ZnReal tmp = 1.0 - u;
+  return 3 * u * (tmp * tmp);
+}
+
+static ZnReal
+B2(ZnReal	u)
+{
+  ZnReal tmp = 1.0 - u;
+  return 3 * u * u * tmp;
+}
+
+static ZnReal
+B3(ZnReal	u)
+{
+  return u * u * u;
+}
+
+/*
+ * ChordLengthParameterize  --
+ *	Assign parameter values to digitized points 
+ *	using relative distances between points.
+ */
+static ZnReal *
+ChordLengthParameterize(ZnPoint	*d,
+			int	first,
+			int	last)
+{
+  int	i;
+  ZnReal	*u;
+
+  u = (ZnReal *) ZnMalloc((unsigned) (last-first+1) * sizeof(ZnReal));
+  
+  u[0] = 0.0;
+  for (i = first+1; i <= last; i++) {
+    u[i-first] = u[i-first-1] + V2DistanceBetween2Points(&d[i], &d[i-1]);
+  }
+  
+  for (i = first + 1; i <= last; i++) {
+    u[i-first] = u[i-first] / u[last-first];
+  }
+  
+  return u;
+}
+
+/*
+ * Bezier --
+ *	Evaluate a Bezier curve at a particular parameter value
+ * 
+ */
+static ZnPoint
+BezierII(int		degree,
+	 ZnPoint	*V,
+	 ZnReal		t)
+{
+  int		i, j;		
+  ZnPoint 	Q;	        /* Point on curve at parameter t	*/
+  ZnPoint 	*Vtemp;		/* Local copy of control points		*/
+  
+  /* Copy array	*/
+  Vtemp = (ZnPoint *) ZnMalloc((unsigned)((degree+1) * sizeof (ZnPoint)));
+  for (i = 0; i <= degree; i++) {
+    Vtemp[i] = V[i];
+  }
+  
+  /* Triangle computation */
+  for (i = 1; i <= degree; i++) {	
+    for (j = 0; j <= degree-i; j++) {
+      Vtemp[j].x = (1.0 - t) * Vtemp[j].x + t * Vtemp[j+1].x;
+      Vtemp[j].y = (1.0 - t) * Vtemp[j].y + t * Vtemp[j+1].y;
+    }
+  }
+  
+  Q = Vtemp[0];
+  ZnFree(Vtemp);
+  return Q;
+}
+
+/*
+ * NewtonRaphsonRootFind --
+ *	Use Newton-Raphson iteration to find better root.
+ */
+static ZnReal
+NewtonRaphsonRootFind(ZnPoint	*Q,
+		      ZnPoint	P,
+		      ZnReal	u)
+{
+  ZnReal 	numerator, denominator;
+  ZnPoint 	Q1[3], Q2[2];		/*  Q' and Q''			*/
+  ZnPoint	Q_u, Q1_u, Q2_u;	/*u evaluated at Q, Q', & Q''	*/
+  ZnReal 	uPrime;			/*  Improved u			*/
+  int 		i;
+    
+  /* Compute Q(u)	*/
+  Q_u = BezierII(3, Q, u);
+    
+  /* Generate control vertices for Q'	*/
+  for (i = 0; i <= 2; i++) {
+    Q1[i].x = (Q[i+1].x - Q[i].x) * 3.0;
+    Q1[i].y = (Q[i+1].y - Q[i].y) * 3.0;
+  }
+    
+  /* Generate control vertices for Q'' */
+  for (i = 0; i <= 1; i++) {
+    Q2[i].x = (Q1[i+1].x - Q1[i].x) * 2.0;
+    Q2[i].y = (Q1[i+1].y - Q1[i].y) * 2.0;
+  }
+  
+  /* Compute Q'(u) and Q''(u)	*/
+  Q1_u = BezierII(2, Q1, u);
+  Q2_u = BezierII(1, Q2, u);
+    
+  /* Compute f(u)/f'(u) */
+  numerator = (Q_u.x - P.x) * (Q1_u.x) + (Q_u.y - P.y) * (Q1_u.y);
+  denominator = (Q1_u.x) * (Q1_u.x) + (Q1_u.y) * (Q1_u.y) +
+                (Q_u.x - P.x) * (Q2_u.x) + (Q_u.y - P.y) * (Q2_u.y);
+    
+  /* u = u - f(u)/f'(u) */
+  uPrime = u - (numerator/denominator);
+  return (uPrime);
+}
+
+/*
+ * Reparameterize --
+ *	Given set of points and their parameterization, try to find
+ *	a better parameterization.
+ */
+static ZnReal *
+Reparameterize(ZnPoint	*d,
+	       int	first,
+	       int	last, 
+	       ZnReal	*u,
+	       ZnPoint	*bezCurve)
+{
+  int		nPts = last-first+1;	
+  int		i;
+  ZnReal	*uPrime;	/*  New parameter values	*/
+
+  uPrime = (ZnReal *) ZnMalloc(nPts * sizeof(ZnReal));
+  for (i = first; i <= last; i++) {
+    uPrime[i-first] = NewtonRaphsonRootFind(bezCurve, d[i], u[i-first]);
+  }
+  return (uPrime);
+}
+
+/*
+ * GenerateBezier --
+ *	Use least-squares method to find Bezier control
+ *	points for region.
+ */
+static void
+GenerateBezier(ZnPoint	*d,
+	       int	first,
+	       int	last, 
+	       ZnReal	*uPrime, 
+	       ZnPoint	tHat1,
+	       ZnPoint	tHat2,
+	       ZnPoint	*bez_curve)
+{
+  int		i;
+  ZnPoint 	*A0, *A1;	/* Precomputed rhs for eqn	*/
+  int		num_points;	/* Number of pts in sub-curve */
+  ZnReal 	C[2][2];	/* Matrix C		*/
+  ZnReal 	X[2];		/* Matrix X			*/
+  ZnReal 	det_C0_C1;	/* Determinants of matrices	*/
+  ZnReal	det_C0_X, det_X_C1;
+  ZnReal 	alpha_l;	/* Alpha values, left and right	*/
+  ZnReal	alpha_r;
+  ZnPoint 	tmp;		/* Utility variable		*/
+  
+  num_points = last - first + 1;
+  A0 = (ZnPoint *) ZnMalloc(num_points * sizeof(ZnPoint));
+  A1 = (ZnPoint *) ZnMalloc(num_points * sizeof(ZnPoint));
+  
+  /* Compute the A's	*/
+  for (i = 0; i < num_points; i++) {
+    ZnPoint	v1, v2;
+    v1 = tHat1;
+    v2 = tHat2;
+    V2Scale(&v1, B1(uPrime[i]));
+    V2Scale(&v2, B2(uPrime[i]));
+    A0[i] = v1;
+    A1[i] = v2;
+  }
+
+  /* Create the C and X matrices	*/
+  C[0][0] = 0.0;
+  C[0][1] = 0.0;
+  C[1][0] = 0.0;
+  C[1][1] = 0.0;
+  X[0]    = 0.0;
+  X[1]    = 0.0;
+
+  for (i = 0; i < num_points; i++) {
+    C[0][0] += V2Dot(&A0[i], &A0[i]);
+    C[0][1] += V2Dot(&A0[i], &A1[i]);
+    C[1][0] = C[0][1];
+    C[1][1] += V2Dot(&A1[i], &A1[i]);
+
+    tmp = V2SubII(d[first + i],
+		  V2AddII(V2ScaleIII(d[first], B0(uPrime[i])),
+			  V2AddII(V2ScaleIII(d[first], B1(uPrime[i])),
+				  V2AddII(V2ScaleIII(d[last], B2(uPrime[i])),
+					  V2ScaleIII(d[last], B3(uPrime[i]))))));
+
+    X[0] += V2Dot(&A0[i], &tmp);
+    X[1] += V2Dot(&A1[i], &tmp);
+  }
+
+  /* Compute the determinants of C and X	*/
+  det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
+  det_C0_X  = C[0][0] * X[1]    - C[0][1] * X[0];
+  det_X_C1  = X[0]    * C[1][1] - X[1]    * C[0][1];
+
+  /* Finally, derive alpha values	*/
+  if (det_C0_C1 == 0.0) {
+    det_C0_C1 = (C[0][0] * C[1][1]) * 10e-12;
+  }
+  alpha_l = det_X_C1 / det_C0_C1;
+  alpha_r = det_C0_X / det_C0_C1;
+
+  /*  If alpha negative, use the Wu/Barsky heuristic (see text) */
+  if (alpha_l < 0.0 || alpha_r < 0.0) {
+    ZnReal dist = V2DistanceBetween2Points(&d[last], &d[first]) / 3.0;
+    
+    bez_curve[0] = d[first];
+    bez_curve[3] = d[last];
+    V2Add(&bez_curve[0], V2Scale(&tHat1, dist), &bez_curve[1]);
+    V2Add(&bez_curve[3], V2Scale(&tHat2, dist), &bez_curve[2]);
+  }
+  else {
+    /*  First and last control points of the Bezier curve are */
+    /*  positioned exactly at the first and last data points */
+    /*  Control points 1 and 2 are positioned an alpha distance out */
+    /*  on the tangent vectors, left and right, respectively */
+    bez_curve[0] = d[first];
+    bez_curve[3] = d[last];
+    V2Add(&bez_curve[0], V2Scale(&tHat1, alpha_l), &bez_curve[1]);
+    V2Add(&bez_curve[3], V2Scale(&tHat2, alpha_r), &bez_curve[2]);
+  }
+  ZnFree(A0);
+  ZnFree(A1);
+}
+
+/*
+ * ComputeMaxError --
+ *	Find the maximum squared distance of digitized points
+ *	to fitted curve.
+*/
+static ZnReal
+ComputeMaxError(ZnPoint	*d,
+		int	first,
+		int	last, 
+		ZnPoint	*bez_curve,
+		ZnReal	*u, 
+		int	*splitPoint)
+{
+  int		i;
+  ZnReal	maxDist;	/*  Maximum error		*/
+  ZnReal	dist;		/*  Current error		*/
+  ZnPoint	P;		/*  Point on curve		*/
+  ZnPoint	v;		/*  Vector from point to curve	*/
+  
+  *splitPoint = (last - first + 1)/2;
+  maxDist = 0.0;
+  for (i = first + 1; i < last; i++) {
+    P = BezierII(3, bez_curve, u[i-first]);
+    v = V2SubII(P, d[i]);
+    dist = V2SquaredLength(&v);
+    if (dist >= maxDist) {
+      maxDist = dist;
+      *splitPoint = i;
+    }
+  }
+  return (maxDist);
+}
+
+/*
+ * ComputeLeftTangent,
+ * ComputeRightTangent,
+ * ComputeCenterTangent --
+ *	Approximate unit tangents at endpoints and
+ *	center of digitized curve.
+ */
+static ZnPoint
+ComputeLeftTangent(ZnPoint	*d,
+		   int		end)
+{
+  ZnPoint tHat1;
+  tHat1 = V2SubII(d[end+1], d[end]);
+  tHat1 = *V2Normalize(&tHat1);
+  return tHat1;
+}
+
+static ZnPoint
+ComputeRightTangent(ZnPoint	*d,
+		    int		end)
+{
+  ZnPoint tHat2;
+  tHat2 = V2SubII(d[end-1], d[end]);
+  tHat2 = *V2Normalize(&tHat2);
+  return tHat2;
+}
+
+
+static ZnPoint
+ComputeCenterTangent(ZnPoint	*d,
+		     int	center)
+{
+  ZnPoint	V1, V2, tHatCenter;
+
+  V1 = V2SubII(d[center-1], d[center]);
+  V2 = V2SubII(d[center], d[center+1]);
+  tHatCenter.x = (V1.x + V2.x)/2.0;
+  tHatCenter.y = (V1.y + V2.y)/2.0;
+  tHatCenter = *V2Normalize(&tHatCenter);
+  return tHatCenter;
+}
+
+static void
+FitCubic(ZnPoint	*d,
+	 int		first,
+	 int		last,
+	 ZnPoint	tHat1, 
+	 ZnPoint	tHat2,
+	 ZnReal		error,
+	 ZnList		controls)
+{
+  ZnPoint	*bez_curve;	/* Control points of fitted Bezier curve*/
+  ZnReal	*u;		/* Parameter values for point  */
+  ZnReal	*uPrime;	/* Improved parameter values */
+  ZnReal	max_err;	/* Maximum fitting error	 */
+  int		splitPoint;	/* Point to split point set at	 */
+  int		num_points;	/* Number of points in subset  */
+  ZnReal	iteration_err;	/* Error below which you try iterating  */
+  int		max_iter = 4;	/* Max times to try iterating  */
+  ZnPoint	tHatCenter;   	/* Unit tangent vector at splitPoint */
+  int		i;		
+
+  iteration_err = error * error;
+  num_points = last - first + 1;
+  ZnListAssertSize(controls, ZnListSize(controls)+4);
+  bez_curve = ZnListAt(controls, ZnListSize(controls)-4);
+  
+  /*  Use heuristic if region only has two points in it */
+  if (num_points == 2) {
+    ZnReal dist = V2DistanceBetween2Points(&d[last], &d[first]) / 3.0;
+
+    bez_curve[0] = d[first];
+    bez_curve[3] = d[last];
+    V2Add(&bez_curve[0], V2Scale(&tHat1, dist), &bez_curve[1]);
+    V2Add(&bez_curve[3], V2Scale(&tHat2, dist), &bez_curve[2]);
+    return;
+  }
+
+  /*  Parameterize points, and attempt to fit curve */
+  u = ChordLengthParameterize(d, first, last);
+  GenerateBezier(d, first, last, u, tHat1, tHat2, bez_curve);
+  
+  /*  Find max deviation of points to fitted curve */
+  max_err = ComputeMaxError(d, first, last, bez_curve, u, &splitPoint);
+  if (max_err < error) {
+    ZnFree(u);
+    return;
+  }
+  
+  /*
+   * If error not too large, try some reparameterization
+   *  and iteration.
+   */
+  if (max_err < iteration_err) {
+    for (i = 0; i < max_iter; i++) {
+      uPrime = Reparameterize(d, first, last, u, bez_curve);
+      GenerateBezier(d, first, last, uPrime, tHat1, tHat2, bez_curve);
+      max_err = ComputeMaxError(d, first, last,
+				bez_curve, uPrime, &splitPoint);
+      if (max_err < error) {
+	ZnFree(u);
+	return;
+      }
+      ZnFree(u);
+      u = uPrime;
+    }
+  }
+  
+  /* Fitting failed -- split at max error point and fit recursively */
+  ZnFree(u);
+  ZnListAssertSize(controls, ZnListSize(controls)-4);
+  tHatCenter = ComputeCenterTangent(d, splitPoint);
+  FitCubic(d, first, splitPoint, tHat1, tHatCenter, error, controls);
+  V2Negate(&tHatCenter);
+  FitCubic(d, splitPoint, last, tHatCenter, tHat2, error, controls);
+}
+
+void
+FitBezier(ZnPoint	*pts,
+	  int		num_points,
+	  ZnReal	error,
+	  ZnList	controls)
+{
+  ZnPoint	tHat1, tHat2;	/*  Unit tangent vectors at endpoints */
+
+  tHat1 = ComputeLeftTangent(pts, 0);
+  tHat2 = ComputeRightTangent(pts, num_points-1);
+  FitCubic(pts, 0, num_points-1, tHat1, tHat2, error, controls);
+}
+
+
+/*
+ **********************************************************************************
+ *
  * GetLineEnd --
  *	Compute the points describing the given line end style at point p1 for
  *	the line p1,p2. Point p1 is adjusted to fit the line end.