aboutsummaryrefslogtreecommitdiff
path: root/generic/Geo.c
diff options
context:
space:
mode:
Diffstat (limited to 'generic/Geo.c')
-rw-r--r--generic/Geo.c618
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.
generated by cgit v1.1 at 2024-11-03 02:44:47 +0000