/*
 * Geo.c -- Implementation of common geometric routines.
 *
 * Authors		: Patrick Lecoanet.
 * Creation date	:
 *
 * $Id$
 */

/*
 *  Copyright (c) 1993 - 1999 CENA, Patrick Lecoanet --
 *
 * This code is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This code is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this code; if not, write to the Free
 * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 */

/*
 * Much of the code here is inspired by (or copied from) the Tk code.
 */


#include "Geo.h"
#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.
 */
void
Anchor2Origin(ZnPoint	*position,
	      ZnDim	width,
	      ZnDim	height,
	      ZnAnchor	anchor,
	      ZnPoint	*origin)
{
  switch (anchor) {
  case ZnAnchorCenter:
    origin->x = position->x - width/2;
    origin->y = position->y - height/2;
    break;
  case ZnAnchorNW:
    *origin = *position;
    break;
  case ZnAnchorN:
    origin->x = position->x - width/2;
    origin->y = position->y;
    break;
  case ZnAnchorNE:
    origin->x = position->x - width;
    origin->y = position->y;
    break;
  case ZnAnchorE:
    origin->x = position->x - width;
    origin->y = position->y - height/2;
    break;
  case ZnAnchorSE:
    origin->x = position->x - width;
    origin->y = position->y - height;
    break;
  case ZnAnchorS:
    origin->x = position->x - width/2;
    origin->y = position->y - height;
    break;
  case ZnAnchorSW:
    origin->x = position->x;
    origin->y = position->y - height;
    break;
  case ZnAnchorW:
    origin->x = position->x;
    origin->y = position->y - height/2;
    break;
  }
}


/*
 * Compute the anchor position given the bbox origin, width,
 * height and the anchor.
 */
void
Origin2Anchor(ZnPoint	*origin,
	      ZnDim	width,
	      ZnDim	height,
	      ZnAnchor	anchor,
	      ZnPoint	*position)
{
  switch (anchor) {
  case ZnAnchorCenter:
    position->x = origin->x + width/2;
    position->y = origin->y + height/2;
    break;
  case ZnAnchorNW:
    *position = *origin;
    break;
  case ZnAnchorN:
    position->x = origin->x + width/2;
    position->y = origin->y;
    break;
  case ZnAnchorNE:
    position->x = origin->x + width;
    position->y = origin->y;
    break;
  case ZnAnchorE:
    position->x = origin->x + width;
    position->y = origin->y + height/2;
    break;
  case ZnAnchorSE:
    position->x = origin->x + width;
    position->y = origin->y + height;
    break;
  case ZnAnchorS:
    position->x = origin->x + width/2;
    position->y = origin->y + height;
    break;
  case ZnAnchorSW:
    position->x = origin->x;
    position->y = origin->y + height;
    break;
  case ZnAnchorW:
    position->x = origin->x;
    position->y = origin->y + height/2;
    break;
  }
}


void
BBox2XRect(ZnBBox	*bbox,
	   XRectangle	*r)
{
  r->x = REAL_TO_INT(bbox->orig.x);
  r->y = REAL_TO_INT(bbox->orig.y);
  r->width = REAL_TO_INT(bbox->corner.x) - r->x;
  r->height = REAL_TO_INT(bbox->corner.y) - r->y;
}


void
GetStringBBox(char	*str,
	      ZnFont	font,
	      ZnPos	x,
	      ZnPos	y,
	      ZnBBox	*str_bbox)
{
  Tk_FontMetrics fm;

  str_bbox->orig.x = x;
  str_bbox->corner.x  = x + ZnTextWidth(font, str, strlen(str));
  Tk_GetFontMetrics(font, &fm);
  str_bbox->orig.y = y - fm.ascent;
  str_bbox->corner.y = str_bbox->orig.y + fm.ascent + fm.descent;
}


void
ResetBBox(ZnBBox *bbox)
{
  bbox->orig.x = bbox->orig.y = 0;
  bbox->corner = bbox->orig;
}


void
CopyBBox(ZnBBox *bbox_from,
	 ZnBBox *bbox_to)
{
  bbox_to->orig = bbox_from->orig;
  bbox_to->corner = bbox_from->corner;
}


void
IntersectBBox(ZnBBox *bbox1,
	      ZnBBox *bbox2,
	      ZnBBox *bbox_inter)
{
  if ((bbox1->corner.x < bbox2->orig.x) ||
      (bbox1->corner.y < bbox2->orig.y) ||
      (bbox2->corner.x < bbox1->orig.x) ||
      (bbox2->corner.y < bbox1->orig.y)) {
    ResetBBox(bbox_inter);
  }
  else {
    bbox_inter->orig.x = MAX(bbox1->orig.x, bbox2->orig.x);
    bbox_inter->orig.y = MAX(bbox1->orig.y, bbox2->orig.y);
    bbox_inter->corner.x = MIN(bbox1->corner.x, bbox2->corner.x);
    bbox_inter->corner.y = MIN(bbox1->corner.y, bbox2->corner.y);
  }
}


ZnBool
IsEmptyBBox(ZnBBox *bbox)
{
  return (bbox->orig.x >= bbox->corner.x) || (bbox->orig.y >= bbox->corner.y);
}


void
AddBBoxToBBox(ZnBBox *bbox,
	      ZnBBox *bbox2)
{
  if (IsEmptyBBox(bbox2)) {
    return;
  }
  if (IsEmptyBBox(bbox)) {
    CopyBBox(bbox2, bbox);
  }
  else {
    bbox->orig.x = MIN(bbox->orig.x, bbox2->orig.x);
    bbox->orig.y = MIN(bbox->orig.y, bbox2->orig.y);
    bbox->corner.x = MAX(bbox->corner.x, bbox2->corner.x);
    bbox->corner.y = MAX(bbox->corner.y, bbox2->corner.y);
  }
}


void
AddPointToBBox(ZnBBox	*bbox,
	       ZnPos	px,
	       ZnPos	py)
{
  if (IsEmptyBBox(bbox)) {
    bbox->orig.x = px;
    bbox->orig.y = py;
    bbox->corner.x = bbox->orig.x + 1;
    bbox->corner.y = bbox->orig.y + 1;
  }
  else {
    bbox->orig.x = MIN(bbox->orig.x, px);
    bbox->orig.y = MIN(bbox->orig.y, py);
    bbox->corner.x = MAX(bbox->corner.x, px + 1);
    bbox->corner.y = MAX(bbox->corner.y, py + 1);
  }
}


void
AddPointsToBBox(ZnBBox	*bbox,
		ZnPoint	*points,
		int	num_points)
{
  int		x1, y1, x2, y2, cur;

  if (points == NULL) {
    return;
  }
  
  if (num_points == 0) {
    return;
  }
  
  if (IsEmptyBBox(bbox)) {
    x1 = points->x;
    y1 = points->y;
    x2 = x1 + 1;
    y2 = y1 + 1;
    num_points--;
    points++;
  }
  else {
    x1 = bbox->orig.x;
    y1 = bbox->orig.y;
    x2 = bbox->corner.x;
    y2 = bbox->corner.y;
  }

  for ( ; num_points > 0; num_points--, points++) {
    cur = points->x;
    if (cur < x1) {
      x1 = cur;
    }
    if (cur > x2) {
      x2 = cur;
    }
    cur = points->y;
    if (cur < y1) {
      y1 = cur;
    }
    if (cur > y2) {
      y2 = cur;
    }
  }
  bbox->orig.x = x1;
  bbox->orig.y = y1;
  bbox->corner.x = x2 + 1;
  bbox->corner.y = y2 + 1;
}


void
AddStringToBBox(ZnBBox	*bbox,
		char	*str,
		ZnFont	font,
		ZnPos	cx,
		ZnPos	cy)
{
  ZnBBox	str_bbox;
  
  GetStringBBox(str, font, cx, cy, &str_bbox);
  AddBBoxToBBox(bbox, &str_bbox);
}

ZnBool
PointInBBox(ZnBBox	*bbox,
	    ZnPos	x,
	    ZnPos	y)
{
  return ((x >= bbox->orig.x) && (x < bbox->corner.x) &&
	  (y >= bbox->orig.y) && (y < bbox->corner.y));
}


/*
 * Tell where aa area is with respect to another area.
 * Return -1 if the first is entirely outside the second,
 * 1 if it is entirely inside and 0 otherwise.
 */
int
BBoxInBBox(ZnBBox	*bbox1,
	   ZnBBox	*bbox2)
{
  if ((bbox1->corner.x <= bbox2->orig.x) ||
      (bbox1->orig.x >= bbox2->corner.x) ||
      (bbox1->corner.y <= bbox2->orig.y) ||
      (bbox1->orig.y >= bbox2->corner.y)) {
    return -1;
  }
  if ((bbox2->orig.x <= bbox1->orig.x) &&
      (bbox1->corner.x <= bbox2->corner.x) &&
      (bbox2->orig.y <= bbox1->orig.y) &&
      (bbox1->corner.y <= bbox2->corner.y)) {
    return 1;
  }
  return 0;
}

/*
 * Tell where a line is with respect to an area.
 * Return -1 if the line is entirely outside the bbox, 1
 * if it is entirely inside and 0 otherwise.
 */
int
LineInBBox(ZnPoint	*p1,
	   ZnPoint	*p2,
	   ZnBBox	*bbox)
{
  ZnBool	p1_inside = PointInBBox(bbox, p1->x, p1->y);
  ZnBool	p2_inside = PointInBBox(bbox, p2->x, p2->y);

  if (p1_inside != p2_inside) {
    return 0;
  }

  if (p1_inside && p2_inside) {
    return 1;
  }
  
  /*
   * Segment may intersect area, check it more thoroughly.
   */
  /* Vertical line */
  if (p1->x == p2->x) {
    if (((p1->y >= bbox->orig.y) ^ (p2->y >= bbox->orig.y)) &&
	(p1->x >= bbox->orig.x) &&
	(p1->x <= bbox->corner.x)) {
      return 0;
    }
  }
  /* Horizontal line */
  else if (p1->y == p2->y) {
    if (((p1->x >= bbox->orig.x) ^ (p2->x >= bbox->orig.x)) &&
	(p1->y >= bbox->orig.y) &&
	(p1->y <= bbox->corner.y)) {
      return 0;
    }
  }
  /* Diagonal, do it the hard way. */
  else {
    double	slope = ((double) p2->y - p1->y) / ((double) p2->x - p1->x);
    int		low, high, x, y;
    int		bbox_width = bbox->corner.x - bbox->orig.x;
    int		bbox_height = bbox->corner.y - bbox->orig.y;

    /* Check against left edge */
    if (p1->x < p2->x) {
      low = p1->x;
      high = p2->x;
    }
    else {
      low = p2->x;
      high = p1->x;
    }

    y = p1->y + (bbox->orig.x - p1->x) * slope;
    if ((bbox->orig.x >= low) && (bbox->orig.x <= high) &&
	(y >= bbox->orig.y) && (y <= bbox->corner.y))
      return 0;

    /* Check against right edge */
    y += bbox_width * slope;
    if ((y >= bbox->orig.y) && (y <= bbox->corner.y) &&
	(bbox->corner.x >= low) && (bbox->corner.x <= high))
      return 0;
	
    /* Check against bottom edge */
    if (p1->y < p2->y) {
      low = p1->y;
      high = p2->y;
    }
    else {
      low = p2->y;
      high = p1->y;
    }

    x = p1->x + (bbox->orig.y - p1->y) / slope;
    if ((x >= bbox->orig.x) && (x <= bbox->corner.x) &&
	(bbox->orig.y >= low) && (bbox->orig.y <= high))
      return 0;

    /* Check against top edge */
    x += bbox_height / slope;
    if ((x >= bbox->orig.x) && (x <= bbox->corner.x) &&
	(bbox->corner.y >= low) && (bbox->corner.y <= high))
      return 0;
  }

  return -1;
}


/*
 * ShiftLine --
 *	Given two points describing a line and a distance, return
 *	to points describing a line parallel to it at the given distance.
 *	When looking the line from p1 to p2 the new line will be dist away
 *	on its left. Negative values are allowed for dist, resulting in a line
 *	on the right.
 */
void
ShiftLine(ZnPoint	*p1,
	  ZnPoint	*p2,
	  ZnReal	dist,
	  ZnPoint	*p3,
	  ZnPoint	*p4)
{
  static int	shift_table[129];
  ZnBool	dx_neg, dy_neg;
  int		dx, dy;

  /*
   * Initialize the conversion table.
   */
  if (shift_table[0] == 0) {
    int		i;
    double	tangent, cosine;
    
    for (i = 0; i <= 128; i++) {
      tangent = i/128.0;
      cosine = 128/cos(atan(tangent)) + 0.5;
      shift_table[i] = (int) cosine;
    }
  }

  *p3 = *p1;
  dx = p2->x - p1->x;
  dy = p2->y - p1->y;
  if (dx < 0) {
    dx = -dx;
    dx_neg = True;
  }
  else {
    dx_neg = False;
  }
  if (dy < 0) {
    dy = -dy;
    dy_neg = True;
  }
  else {
    dy_neg = False;
  }
  if ((dy < PRECISION_LIMIT) && (dx < PRECISION_LIMIT)) {
    printf("ShiftLine: segment is a point\n");
    return;
  }

  if (dy <= dx) {
    dy = ((dist * shift_table[(dy*128)/dx]) + 64) / 128;
    if (!dx_neg) {
      dy = -dy;
    }
    p3->y += dy;
  }
  else {
    dx = ((dist * shift_table[(dx*128)/dy]) + 64) / 128;
    if (dy_neg) {
      dx = -dx;
    }
    p3->x += dx;
  }

  p4->x = p3->x + (p2->x - p1->x);
  p4->y = p3->y + (p2->y - p1->y);
}


/*
 * IntersectLines --
 *	Given two lines described by two points, compute their intersection.
 *	The function returns True if the lines are not parallel and False
 *	otherwise.
 */
ZnBool
IntersectLines(ZnPoint	*a1,
	       ZnPoint	*a2,
	       ZnPoint	*b1,
	       ZnPoint	*b2,
	       ZnPoint	*pi)
{
  ZnReal	dxadyb, dxbdya, dxadxb, dyadyb, p, q;

  dxadyb = (a2->x - a1->x)*(b2->y - b1->y);
  dxbdya = (b2->x - b1->x)*(a2->y - a1->y);
  dxadxb = (a2->x - a1->x)*(b2->x - b1->x);
  dyadyb = (a2->y - a1->y)*(b2->y - b1->y);
  
  if (dxadyb == dxbdya) {
    return False;
  }
  
  p = a1->x*dxbdya - b1->x*dxadyb + (b1->y - a1->y)*dxadxb;
  q = dxbdya - dxadyb;
  if (q < 0) {
    p = -p;
    q = -q;
  }
  if (p < 0) {
    pi->x = -(-2*p + q)/(2*q); /*-(-p + q/2)/q;*/
  }
  else {
    pi->x = (2*p + q)/(2*q); /*(p + q/2)/q;*/
  }
  
  p = a1->y*dxadyb - b1->y*dxbdya + (b1->x - a1->x)*dyadyb;
  q = dxadyb - dxbdya;
  if (q < 0) {
    p = -p;
    q = -q;
  }
  if (p < 0) {
    pi->y = -(-2*p + q)/(2*q); /*-(-p + q/2)/q;*/
  }
  else {
    pi->y = (2*p + q)/(2*q); /*(p + q/2)/q;*/
  }
  
  return True;
}


/*
 * InsetPolygon --
 *	Inset the given polygon by the given amount. The
 *     value can be negative, in this case the polygon will
 *     be outset.
 */
/**** A FINIR ****/
void
InsetPolygon(ZnPoint	*p,
	     int	num_points,
	     ZnDim	inset)
{
  ZnPoint	*p1, *p2;
  ZnPoint	new_p1, new_p2;
  /*  ZnPoint	shift1, shift2;*/
  int		i, processed_points;

  processed_points = 0;
  
  if ((p->x == p[num_points-1].x) && (p->y == p[num_points-1].y)) {
    num_points--;
  }
  for (p1 = p, p2 = p1+1, i = 0; i < num_points; i++, p1 = p2, p2++) {
    /*
     * Wrap to the first point.
     */
    if (i == num_points-1) {
      p2 = p;
    }
    /*
     * Skip duplicate vertices.
     */
    if ((p2->x == p1->x) && (p2->y == p1->y)) {
      continue;
    }
    
    ShiftLine(p1, p2, inset, &new_p1, &new_p2);

    if (processed_points >= 1) {
    }
  }
}


/*
 * Compute the two corner points of a thick line end.
 * Two points describing the line segment and the width
 * are given as input. If projecting is true this function
 * mimics the X11 line projecting behaviour. The computed
 * end is located around p2.
 */
void
GetButtPoints(ZnPoint	*p1,
	      ZnPoint	*p2,
	      int	width,
	      ZnBool	projecting,
	      ZnPoint	*c1,
	      ZnPoint	*c2)
{
  double	w_2 = width/2.0;
  double	length = hypot(p2->x - p1->x, p2->y - p1->y);
  double	delta_x, delta_y;
  
  if (length == 0.0) {
    c1->x = c2->x = p2->x;
    c1->y = c2->y = p2->y;
  }
  else {
    delta_x = -w_2 * (p2->y - p1->y) / length;
    delta_y = w_2 * (p2->x - p1->x) / length;
    c1->x = p2->x + delta_x;
    c2->x = p2->x - delta_x;
    c1->y = p2->y + delta_y;
    c2->y = p2->y - delta_y;
    if (projecting) {
      c1->x += delta_y;
      c2->x += delta_y;
      c1->y -= delta_x;
      c2->y -= delta_x;
    }
  }
}

/*
 * Compute the inside and outside points of the mitered
 * corner formed by a thick line going through 3 points.
 * If the angle formed by the three points is less than
 * 11 degrees, False is returned an no points are computed.
 * Else True is returned and the points are in c1, c2.
 *
 * If someday the code is switched to REAL coordinates, we
 * must round each coordinate to the nearer integer to mimic
 * the way pixels are drawn. Sample code: floor(p->x+0.5);
 *
 * Hmmm, the switch has been done but not the rounding ;-)
 */
ZnBool
GetMiterPoints(ZnPoint	*p1,
	       ZnPoint	*p2,
	       ZnPoint	*p3,
	       int	width,
	       ZnPoint	*c1,
	       ZnPoint	*c2)
{
  static double	deg11 = (11.0*2.0*M_PI)/360.0;
  double	theta1;		/* angle of p2-p1 segment. */
  double	theta2;		/* angle of p2-p3 segment. */
  double	theta;		/* angle of the joint */
  double	theta3;		/* angle of bisector of the joint toward
				 * the external point of the joint. */
  double	dist;		/* distance of the external points
				 * of the corner from the mid point
				 * p2. */
  double	delta_x, delta_y; /* projection of (dist,theta3) on x
				   * and y. */
  
  if (p2->y == p1->y) {
    theta1 = (p2->x < p1->x) ? 0.0 : M_PI;
  }
  else if (p2->x == p1->x) {
    theta1 = (p2->y < p1->y) ? M_PI/2.0 : -M_PI/2.0;
  }
  else {
    theta1 = atan2(p1->y - p2->y, p1->x - p2->x);
  }
  if (p3->y == p2->y) {
    theta2 = (p3->x > p2->x) ? 0.0 : M_PI;
  }
  else if (p3->x == p2->x) {
    theta2 = (p3->y > p2->y) ? M_PI/2.0 : -M_PI/2.0;
  }
  else {
    theta2 = atan2(p3->y - p2->y, p3->x - p2->x);
  }
  theta = theta1 - theta2;
  if (theta > M_PI) {
    theta -= 2.0*M_PI;
  }
  else if (theta < -M_PI) {
    theta += 2*M_PI;
  }
  if ((theta < deg11) && (theta > -deg11)) {
    return False;
  }
  /*
   * Compute the distance of the internal and external
   * corner points from the intersection p2 (considered
   * at 0,0).
   */
  dist = 0.5*width / sin(0.5*theta);
  dist = ABS(dist);

  /*
   * Compute the angle of the bisector of the joint that
   * goes toward the outside of the joint (the left hand
   * when looking from p1-p2).
   */
  theta3 = (theta1 + theta2)/2.0;
  if (sin(theta3 - (theta1 + M_PI)) < 0.0) {
    theta3 += M_PI;
  }
  
  delta_x = dist * cos(theta3);
  c1->x = p2->x + delta_x;
  c2->x = p2->x - delta_x;
  delta_y = dist * sin(theta3);
  c1->y = p2->y + delta_y;
  c2->y = p2->y - delta_y;
  
  return True;
}

/*
 * Tell where a thick polyline is with respect to an area.
 * Return -1 if the polyline is entirely outside the bbox, 1
 * if it is entirely inside and 0 otherwise. The joints can
 * be specified as JoinMiter, JoinRound, JoinBevel. The ends
 * can be: CapRound, CapButt, CapProjecting.
 */
int
PolylineInBBox(ZnPoint	*points,
	       int	num_points,
	       int	width,
	       int	cap_style,
	       int	join_style,
	       ZnBBox	*bbox)
{
  int		count, inside = -1;
  ZnBool	do_miter_as_bevel;
  ZnPoint	poly[4];
  
  /*
   * If the first point is inside the area, change inside
   * accordingly.
   */
  if ((points[0].x >= bbox->orig.x) && (points[0].x <= bbox->corner.x) &&
      (points[0].y >= bbox->orig.y) && (points[0].y <= bbox->corner.y)) {
    inside = 1;
  }

  /*
   * Now iterate through all the edges. Compute a polygon for
   * each and test it against the area. At each vertex an oval
   * of radius width/2 is also tested to account for round ends
   * and joints.
   */
  do_miter_as_bevel = False;
  for (count = num_points; count >= 2; count--, points++) {

    /*
     * Test a circle around the first point if CapRound or
     * around every joint if JoinRound.
     */
    if (((cap_style == CapRound) && (count == num_points)) ||
	((join_style == JoinRound) && (count != num_points))) {
      if (OvalInBBox(points, width, width, bbox) != inside) {
	return 0;
      }
    }
    /*
     * Build a polygon to represent an edge from the current
     * point to the next. Special cases for the first and
     * last edges to implement the line ends.
     */
    /*
     * First vertex of the edge
     */
    if (count == num_points) {
      GetButtPoints(&points[1], points, width,
		    cap_style == CapProjecting, poly, &poly[1]);
    }
    /*
     * Here we are at a joint starting a new edge. If the
     * joint is mitered, start by carrying over the points
     * from the previous edge. Otherwise compute new points
     * for a butt end.
     */
    else if ((join_style == JoinMiter) && !do_miter_as_bevel) {
      poly[0] = poly[3];
      poly[1] = poly[2];
    }
    else {
      GetButtPoints(&points[1], points, width, 0, poly, &poly[1]);
      /*
       * If the previous joint was beveled (or considered so),
       * check the polygon that fill the bevel. It has more or
       * less an X shape, i.e, it's self intersecting. If this
       * is not ok, it may be necessary to permutte poly[1] &
       * poly[2].
       */
      if ((join_style == JoinBevel) || do_miter_as_bevel) {
	if (PolygonInBBox(poly, 4, bbox, NULL) != inside) {
	  return 0;
	}
	do_miter_as_bevel = False;
      }
    }

    /*
     * Opposite vertex of the edge.
     */
    if (count == 2) {
      GetButtPoints(points, &points[1], width, cap_style == CapProjecting,
		    &poly[2], &poly[3]);
    }
    else if (join_style == JoinMiter) {
      if (GetMiterPoints(points, &points[1], &points[2], width,
			 &poly[2], &poly[3]) == False) {
	do_miter_as_bevel = True;
	GetButtPoints(points, &points[1], width, 0, &poly[2], &poly[3]);
      }
    }
    else {
      GetButtPoints(points, &points[1], width, 0, &poly[2], &poly[3]);
    }

    if (PolygonInBBox(poly, 4, bbox, NULL) != inside) {
      return 0;
    }
  }

  /*
   * Test a circle around the last point if CapRound.
   */
  if (cap_style == CapRound) {
    if (OvalInBBox(points, width, width, bbox) != inside) {
      return 0;
    }
  }

  return inside;
}


/*
 * 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 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,
	      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;
  /*
   * Check to see if closed. If not adjust num_points and
   * record this.
   */
  if ((points[0].x != points[num_points-1].x) ||
      (points[0].y != points[num_points-1].y)) {
    closed = False;
    count = num_points-1;
  }

  /*
   * Get the status of the first edge.
   */
  inside = LineInBBox(&p[0], &p[1], bbox);
  p++;
  if (inside == 0) {
    return 0;
  }
  for (; count > 0; p++, count--) {
    first = &p[0];
    /*
     * Pretend the polygon is closed if this is not the case.
     */
    if (!closed && (count == 1)) {
      second = head;
    }
    else {
      second = &p[1];
    }
    
    if (LineInBBox(first, second, bbox) != inside) {
      return 0;
    }
  }

  /*
   * If all the edges are inside the area, the polygon is
   * inside the area. If all the edges are outside, the polygon
   * may completely enclose the area. Test if the origin of
   * the area is inside the polygon to decide this.
   */
  if (inside == 1) {
    return 1;
  }

  /*printf("PolygonInBBox, np = %d, x = %g, y = %g, dist = %g\n",
	 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;
  }
  
  return -1;
}


/*
 * Tell where an oval is with respect to an area.
 * Return -1 if the oval is entirely outside the bbox, 1
 * if it is entirely inside and 0 otherwise.
 */
int
OvalInBBox(ZnPoint	*center,
	   int		width,
	   int		height,
	   ZnBBox	*bbox)
{
  ZnPoint	origin, corner;
  int		w_2, h_2;
  double	delta_x, delta_y;
  
  w_2 = (width+1)/2;
  h_2 = (height+1)/2;
  
  origin.x = center->x - w_2;
  origin.y = center->y - h_2;
  corner.x = center->x + w_2;
  corner.y = center->y + h_2;

  /*
   * This if the oval bbox is completely inside or outside
   * of the area. Trivial case first.
   */
  if ((bbox->orig.x <= origin.x) && (bbox->corner.x >= corner.x) &&
      (bbox->orig.y <= origin.y) && (bbox->corner.y >= corner.y)) {
    return 1;
  }
  if ((bbox->corner.x < origin.x) || (bbox->orig.x > corner.x) ||
      (bbox->corner.y < origin.y) || (bbox->orig.y > corner.y)) {
    return -1;
  }

  /*
   * Then test all sides of the area against the oval center.
   * If the point of a side closest to the center is inside
   * the oval, then the oval intersects the area.
   */

  /*
   * Compute the square of the Y axis distance, reducing
   * the oval to a unit circle to take into account the
   * shape factor.
   */
  delta_y = bbox->orig.y - center->y;
  if (delta_y < 0.0) {
    delta_y = center->y - bbox->corner.y;
    if (delta_y < 0.0) {
      delta_y = 0.0;
    }
  }
  delta_y /= h_2;
  delta_y *= delta_y;

  /*
   * Test left and then right edges.
   */
  delta_x = (bbox->orig.x - center->x) / w_2;
  delta_x *= delta_x;
  if ((delta_x + delta_y) <= 1.0) {
    return 0;
  }
  delta_x = (bbox->corner.x - center->x) / w_2;
  delta_x *= delta_x;
  if ((delta_x + delta_y) <= 1.0) {
    return 0;
  }
  
  /*
   * Compute the square of the X axis distance, reducing
   * the oval to a unit circle to take into account the
   * shape factor.
   */
  delta_x = bbox->orig.x - center->x;
  if (delta_x < 0.0) {
    delta_x = center->x - bbox->corner.x;
    if (delta_x < 0.0) {
      delta_x = 0.0;
    }
  }
  delta_x /= w_2;
  delta_x *= delta_x;
  
  /*
   * Test top and then bottom edges.
   */
  delta_y = (bbox->orig.y - center->y) / h_2;
  delta_y *= delta_y;
  if ((delta_x + delta_y) <= 1.0) {
    return 0;
  }
  delta_y = (bbox->corner.y - center->y) / h_2;
  delta_y *= delta_y;
  if ((delta_x + delta_y) <= 1.0) {
    return 0;
  }

  return -1;
}


/*
 * Tell if a point is in an angular range whose center is 0,0.
 * The range is specified by a starting angle and an angle extent.
 * The use of a double here is important, don't change it. In some
 * case we need to normalize a figure to take care of its shape and
 * the result needs precision.
 */
ZnBool
PointInAngle(int	start_angle,
	     int	angle_extent,
	     ZnPoint	*p)
{
  double	point_angle;
  int		angle_diff;

  if ((p->x == 0) && (p->y == 0)) {
    point_angle = 0.0;
  }
  else {
    point_angle = atan2(p->y, p->x) * 180.0 / M_PI;
  }
  angle_diff = (REAL_TO_INT(point_angle) - start_angle) % 360;
  if (angle_diff < 0) {
    angle_diff += 360;
  }
  return ((angle_diff <= angle_extent) ||
	  ((angle_extent < 0) && ((angle_diff - 360) >= angle_extent)));
}

/*
 * PointPolarToCartesian --
 *	Convert a point in polar coordinates (rho, theta)
 *	in a reference system described by angle heading
 *	(angles running clockwise) to a point in cartesian
 *	coordinates (delta_x, delta_y).
 *
 */
void
PointPolarToCartesian(ZnReal	heading,
		      int	rho,
		      int	theta,
		      int	*delta_x,
		      int	*delta_y)
{
  double	to_angle;

  /* Compute angle in trigonometric system */
  to_angle = DegreesToRadian(theta) + heading - M_PI_2;
  /*  to_angle = heading - DegreesToRadian(theta);*/
  /* Compute cartesian coordinates */
  *delta_x = (int) (rho * cos(to_angle));
  *delta_y = (int) (rho * sin(to_angle));
}

/*
 * Return a vector angle given its projections
 */
ZnReal
ProjectionToAngle(ZnDim	dx,
		  ZnDim	dy)
{
  if (dx == 0) {
    if (dy < 0) {
      return -M_PI_2;
    }
    else if (dy > 0) {
      return M_PI_2;
    }
    else {
      return 0.0;
    }
  }
  else if (dx < 0) {
    return atan((double) dy / (double) dx) - M_PI;
  }
  else {
    return atan((double) dy / (double) dx);
  }
  return 0.0;
}


/*
 * Tell if an horizontal line intersect an axis aligned
 * elliptical arc.
 *
 * Returns True if the line described by (x1, x2, y) intersects
 * the arc described by (r1, r2, start_angle and angle_extent).
 * This arc is origin centered.
 */
ZnBool
HorizLineToArc(ZnReal	x1,
	       ZnReal	x2,
	       ZnReal	y,
	       ZnReal	rx,
	       ZnReal	ry,
	       ZnReal	start_angle,
	       ZnReal	angle_extent)
{
  ZnReal	tmp, x;
  ZnPoint	t;
  
  /*
   * Compute the x-coordinate of one possible intersection point
   * between the arc and the line.  Use a transformed coordinate
   * system where the oval is a unit circle centered at the origin.
   * Then scale back to get actual x-coordinate.
   */
  t.y = y/ry;
  tmp = 1 - t.y*t.y;
  if (tmp < 0.0) {
    return False;
  }
  t.x = sqrt(tmp);
  x = t.x*rx;
  
  /*
   * Test both intersection points.
   */  
  if ((x >= x1) && (x <= x2) && PointInAngle(start_angle, angle_extent, &t)) {
    return True;
  }
  t.x = -t.x;
  if ((-x >= x1) && (-x <= x2) && PointInAngle(start_angle, angle_extent, &t)) {
    return True;
  }
  return False;
}


/*
 * Tell if an vertical line intersect an axis aligned
 * elliptical arc.
 *
 * Returns True if the line described by (x1, x2, y) intersects
 * the arc described by (r1, r2, start_angle and angle_extent).
 * This arc is origin centered.
 */
ZnBool
VertLineToArc(ZnReal	x,
	      ZnReal	y1,
	      ZnReal	y2,
	      ZnReal	rx,
	      ZnReal	ry,
	      ZnReal	start_angle,
	      ZnReal	angle_extent)
{
  ZnReal	tmp, y;
  ZnPoint	t;
  
  /*
   * Compute the y-coordinate of one possible intersection point
   * between the arc and the line.  Use a transformed coordinate
   * system where the oval is a unit circle centered at the origin.
   * Then scale back to get actual y-coordinate.
   */
  t.x = x/rx;
  tmp = 1 - t.x*t.x;
  if (tmp < 0.0) {
    return False;
  }
  t.y = sqrt(tmp);
  y = t.y*ry;

  /*
   * Test both intersection points.
   */
  if ((y > y1) && (y < y2) && PointInAngle(start_angle, angle_extent, &t)) {
    return True;
  }
  t.y = -t.y;
  if ((-y > y1) && (-y < y2) && PointInAngle(start_angle, angle_extent, &t)) {
    return True;
  }
  return False;
}


/*
 * Return the distance of the given point to the rectangle
 * described by rect. Return negative values for points in
 * the rectangle.
 */
double
RectangleToPointDist(ZnBBox	*bbox,
		     ZnPoint	*p)
{
  double	new_dist, dist;
  ZnPoint	p1, p2;

  p1.x = bbox->orig.x;
  p1.y = p2.y = bbox->orig.y;
  p2.x = bbox->corner.x;
  dist = LineToPointDist(&p1, &p2, p);
  if (dist == 0.0) {
    return 0.0;
  }

  p1 = p2;
  p2.y = bbox->corner.y;
  new_dist = LineToPointDist(&p1, &p2, p);
  dist = MIN(dist, new_dist);
  if (dist == 0.0) {
    return 0.0;
  }

  p1 = p2;
  p2.x = bbox->orig.x;
  new_dist = LineToPointDist(&p1, &p2, p);
  dist = MIN(dist, new_dist);
  if (dist == 0.0) {
    return 0.0;
  }
  
  p1 = p2;
  p2.y = bbox->orig.y;
  new_dist = LineToPointDist(&p1, &p2, p);
  dist = MIN(dist, new_dist);

  if (PointInBBox(bbox, p->x, p->y)) {
    return -dist;
  }
  else {
    return dist;
  }
}


/*
 * Return the distance of the given point to the line
 * described by <xl1,yl1>, <xl2,yl2>..
 */
double
LineToPointDist(ZnPoint	*p1,
		ZnPoint	*p2,
		ZnPoint	*p)
{
  double	x, y;
  int		x_int, y_int;

  /*
   * First compute the closest point on the line. This is done
   * separately for vertical, horizontal, other lines.
   */

  /* Vertical */
  if (p1->x == p2->x) {
    x = p1->x;
    if (p1->y >= p2->y) {
      y_int = MIN(p1->y, p->y);
      y_int = MAX(y_int, p2->y);
    }
    else {
      y_int = MIN(p2->y, p->y);
      y_int = MAX(y_int, p1->y);
    }
    y = y_int;
  }

  /* Horizontal */
  else if (p1->y == p2->y) {
    y = p1->y;
    if (p1->x >= p2->x) {
      x_int = MIN(p1->x, p->x);
      x_int = MAX(x_int, p2->x);
    }
    else {
      x_int = MIN(p2->x, p->x);
      x_int = MAX(x_int, p1->x);
    }
    x = x_int;
  }

  /*
   * Other. Compute its parameters of form y = a1*x + b1 and
   * then compute the parameters of the perpendicular passing
   * through the point y = a2*x + b2, last find the closest point
   * on the segment.
   */
  else {
    double	a1, a2, b1, b2;

    a1 = ((double) (p2->y - p1->y)) / ((double) (p2->x - p1->x));
    b1 = p1->y - a1*p1->x;

    a2 = -1.0/a1;
    b2 = p->y - a2*p->x;

    x = (b2 - b1) / (a1 - a2);
    y = a1*x + b1;

    if (p1->x > p2->x) {
      if (x > p1->x) {
	x = p1->x;
	y = p1->y;
      }
      else if (x < p2->x) {
	x = p2->x;
	y = p2->y;
      }
    }
    else {
      if (x > p2->x) {
	x = p2->x;
	y = p2->y;
      }
      else if (x < p1->x) {
	x = p1->x;
	y = p1->y;
      }
    }
  }

  /* Return the distance */
  return hypot(p->x - x, p->y - y);
}


/*
 * Return the distance of the polygon described by
 * points, to the given point. If the point is
 * inside return values are negative.
 */
double
PolygonToPointDist(ZnPoint	*points,
		   int		num_points,
		   ZnPoint	*p)
{
  double	best_distance;
  int		intersections;
  int		x_int, y_int;
  ZnPoint	*first_point;
  double	x, y, dist;
  ZnPoint	p1, p2;

  /*
   * The algorithm iterates through all the edges of the polygon
   * computing for each the distance to the point and whether a vertical
   * ray starting at the point, intersects the edge. The smallest
   * distance of all edges is stored in best_distance while intersections
   * hold the count of edges to rays intersections. For more informations
   * on how the distance is computed see LineToPointDist.
   */
  best_distance = 1.0e40;
  intersections = 0;

  first_point = points;

  /*
   * Check to see if closed. Adjust num_points to open it (the
   * algorithm always consider a set of points as a closed polygon).
   */
  if ((points[0].x == points[num_points-1].x) &&
      (points[0].y == points[num_points-1].y)) {
    num_points--;
  }

  for ( ; num_points >= 1; num_points--, points++) {
    p1 = points[0];
    /*
     * Wrap over to the first point.
     */
    if (num_points == 1) {
      p2 = *first_point;
    }
    else {
      p2 = points[1];
    }
    
    /*
     * First try to find the closest point on this edge.
     */

    /* Vertical edge */
    if (p1.x == p2.x) {
      x = p1.x;
      if (p1.y >= p2.y) {
	y_int = MIN(p1.y, p->y);
	y_int = MAX(y_int, p2.y);
      }
      else {
	y_int = MIN(p2.y, p->y);
	y_int = MAX(y_int, p1.y);
      }
      y = y_int;
    }

    /* Horizontal edge */
    else if (p1.y == p2.y) {
      y = p1.y;
      if (p1.x >= p2.x) {
	x_int = MIN(p1.x, p->x);
	x_int = MAX(x_int, p2.x);
	if ((p->y < y) && (p->x < p1.x) && (p->x >= p2.x)) {
	  intersections++;
	}
      }
      else {
	x_int = MIN(p2.x, p->x);
	x_int = MAX(x_int, p1.x);
	if ((p->y < y) && (p->x < p2.x) && (p->x >= p1.x)) {
	  intersections++;
	}
      }
      x = x_int;
    }

    /* Other */
    else {
      double	a1, b1, a2, b2;

      a1 = ((double) (p2.y - p1.y)) / ((double) (p2.x - p1.x));
      b1 = p1.y - a1 * p1.x;

      a2 = -1.0/a1;
      b2 = p->y - a2 * p->x;

      x = (b2 - b1)/(a1 - a2);
      y = a1 * x + b1; 

      if (p1.x > p2.x) {
	if (x > p1.x) {
	  x = p1.x;
	  y = p1.y;
	}
	else if (x < p2.x) {
	  x = p2.x;
	  y = p2.y;
	}
      }
      else {
	if (x > p2.x) {
	  x = p2.x;
	  y = p2.y;
	}
	else if (x < p1.x) {
	  x = p1.x;
	  y = p1.y;
	}
      }

      if (((a1 * p->x + b1) > p->y) &&	/* True if point is lower */
	  (p->x >= MIN(p1.x, p2.x)) &&
	  (p->x < MAX(p1.x, p2.x))) {
	intersections++;
      }
    }

    /*
     * Now compute the distance to the closest point and
     * keep it if it is the shortest.
     */
    dist = hypot(p->x - x, p->y - y);
    best_distance = MIN(best_distance, dist);
    /*
     * We can safely escape here if distance is zero.
     */
    if (best_distance == 0.0) {
      return 0.0;
    }
  }

  /*
   * Well, all the edges are processed, if the
   * intersection count is odd the point is inside.
   */
  if (intersections & 0x1) {
    return -best_distance;
  }
  else {
    return best_distance;
  }
}


/*
 * Return the distance of a thick polyline to the
 * given point. Cap and Join parameters are considered
 * in the process.
 */
double
PolylineToPointDist(ZnPoint	*points,
		    int		num_points,
		    int		width,
		    int		cap_style,
		    int		join_style,		   
		    ZnPoint	*p)
{
  ZnBool	miter2bevel = False;
  int		count;
  ZnPoint	*ptr;
  ZnPoint	outline[5];
  double	dist, best_dist, h_width;

  best_dist = 1.0e36;
  h_width = width/2.0;

  for (count = num_points, ptr = points; count >= 2; count--, ptr++) {
    if (((cap_style == CapRound) && (count == num_points)) ||
	((join_style == JoinRound) && (count != num_points))) {
      dist = hypot(ptr->x - p->x, ptr->y - p->y) - h_width;
      if (dist <= 0.0) {
	best_dist = 0.0;
	goto done;
      }
      else if (dist < best_dist) {
	best_dist = dist;
      }
    }
    /*
     * Build the polygonal outline of the current edge.
     */
    if (count == num_points) {
      GetButtPoints(&ptr[1], ptr, width, cap_style==CapProjecting, outline, &outline[1]);
    }
    else if ((join_style == JoinMiter) && !miter2bevel) {
      outline[0] = outline[3];
      outline[1] = outline[2];
    }
    else {
      GetButtPoints(&ptr[1], ptr, width, 0, outline, &outline[1]);
      /*
       * If joints are beveled, check the distance to the polygon
       * that fills the joint.
       */
      if ((join_style == JoinBevel) || miter2bevel) {
	outline[4] = outline[0];
	dist = PolygonToPointDist(outline, 5, p);
	if (dist <= 0.0) {
	  best_dist = 0.0;
	  goto done;
	}
	else if (dist < best_dist) {
	  best_dist = dist;
	}
	miter2bevel = False;
      }
    }
    if (count == 2) {
      GetButtPoints(ptr, &ptr[1], width, cap_style==CapProjecting,
		    &outline[2], &outline[3]);
    }
    else if (join_style == JoinMiter) {
      if (GetMiterPoints(ptr, &ptr[1], &ptr[2], width,
			 &outline[2], &outline[3]) == False) {
	miter2bevel = True;
	GetButtPoints(ptr, &ptr[1], width, 0, &outline[2], &outline[3]);
      }
      /*printf("2=%g+%g, 3=%g+%g\n",
	outline[2].x, outline[2].y, outline[3].x, outline[3].y);*/
    }
    else {
      GetButtPoints(ptr, &ptr[1], width, 0, &outline[2], &outline[3]);
    }
    outline[4] = outline[0];
    /*printf("0=%g+%g, 1=%g+%g, 2=%g+%g, 3=%g+%g, 4=%g+%g\n",
	   outline[0].x, outline[0].y, outline[1].x, outline[1].y,
	   outline[2].x, outline[2].y, outline[3].x, outline[3].y,
	   outline[4].x, outline[4].y);*/
    dist = PolygonToPointDist(outline, 5, p);
    if (dist <= 0.0) {
      best_dist = 0.0;
      goto done;
    }
    else if (dist < best_dist) {
      best_dist = dist;
    }
  }

  /*
   * Test the final point if cap style is round. The code so far
   * has only handled the butt and projecting cases.
   */
  if (cap_style == CapRound) {
    dist = hypot(ptr->x - p->x, ptr->y - p->y) - h_width;
    if (dist <= 0.0) {
      best_dist = 0.0;
      goto done;
    }
    else if (dist < best_dist) {
      best_dist = dist;
    }
  }
  
 done:
  return best_dist;
}


/*
 * Return the distance of the given oval to the point given.
 * The oval is described by its bounding box <xbb,ybb,wbb,hbb>,
 * the thickness of its outline <width>. Return values are negative
 * if the point is inside.
 */
double
OvalToPointDist(ZnPoint		*center,
		int		width,
		int		height,
		unsigned int	line_width,
		ZnPoint		*p)
{
  double	x_delta, y_delta;
  /*  double	x_diameter, y_diameter;*/
  double	scaled_distance;
  double	distance_to_outline;
  double	distance_to_center;

  /*
   * Compute the distance from the point given to the center
   * of the oval. Then compute the same distance in a coordinate
   * system where the oval is a circle with unit radius.
   */

  x_delta = p->x - center->x;
  y_delta = p->y - center->y;
  distance_to_center = hypot(x_delta, y_delta);
  scaled_distance = hypot(x_delta / ((width + line_width) / 2.0),
			  y_delta / ((height + line_width) / 2.0));

  /*
   * If the scaled distance is greater than 1.0 the point is outside
   * the oval. Compute the distance to the edge and convert it back
   * to the original coordinate system. This distance is not much
   * accurate and can overestimate the real distance if the oval is
   * very eccentric.
   */
  if (scaled_distance > 1.0) {
    distance_to_outline = (distance_to_center / scaled_distance) * (scaled_distance - 1.0);
    return distance_to_outline;
  }

  /*
   * The point is inside the oval. Compute the distance as above and check
   * if the point is within the outline.
   */
  if (scaled_distance > 1.0e-10) {
    distance_to_outline = (distance_to_center / scaled_distance) * (1.0 - scaled_distance) - line_width;
  }
  else {
    /*
     * If the point is very close to the center avoid dividing by a
     * very small number, take another method.
     */
    if (width < height)
      distance_to_outline = ((double) (width - line_width)) / 2;
    else
      distance_to_outline = ((double) (height - line_width)) / 2;
  }

  if (distance_to_outline < 0.0)
    return 0.0;
  else
    return -distance_to_outline;
}


static int bezier_basis[][4] =
{
    {	-1,	3,     -3,	1},
    {	 3,    -6,	3,	0},
    {	-3,	3,	0,	0},
    {	 1,	0,	0,	0}
};

/*
 **********************************************************************************
 *
 * Arc2Param --
 *
 *	Given a Bezier curve describing an arc and an angle return the parameter
 *	value for the intersection point between the arc and the ray at angle.
 *
 **********************************************************************************
 */
#define EVAL(coeff, t) (((coeff[0]*t + coeff[1])*t + coeff[2]) * t + coeff[3])
static ZnReal
Arc2Param(ZnPoint	*controls,
	  ZnReal	angle)
{
  ZnReal	coeff_x[4], coeff_y[4];
  ZnReal	min_angle, min_t, max_angle, max_t, cur_angle, cur_t;
  int		i, j, depth = 0;

  /* assume angle >= 0 */
  while (angle > M_PI) {
    angle -= 2 * M_PI;
  }

  for (i = 0; i < 4; i++) {
    coeff_x[i] = coeff_y[i] = 0.0;
    for (j = 0; j < 4; j++) {
      coeff_x[i] += bezier_basis[i][j] * controls[j].x;
      coeff_y[i] += bezier_basis[i][j] * controls[j].y;
    }
  }

  min_angle = atan2(controls[0].y, controls[0].x);
  max_angle = atan2(controls[3].y, controls[3].x);
  if (max_angle < min_angle) {
    min_angle -= M_PI+M_PI;
  }
  if (angle > max_angle) {
    angle -= M_PI+M_PI;
  }

  min_t = 0.0; max_t = 1.0;

  while (depth < 15) {
    cur_t = (max_t + min_t) / 2.0;
    cur_angle = atan2(EVAL(coeff_y, cur_t), EVAL(coeff_x, cur_t));
    if (angle > cur_angle) {
      min_t = cur_t;
      min_angle = cur_angle;
    }
    else {
      max_t = cur_t;
      max_angle = cur_angle;
    }
    depth += 1;
  }

  if ((max_angle-angle) < (angle-min_angle)) {
    return max_t;
  }

  return min_t;
}
#undef EVAL


/*
 **********************************************************************************
 *
 * BezierSubdivide --
 *
 *	Subdivide a Bezier curve given by controls at parameter t. Return
 *	in controls, the first or the last part depending on boolean first.
 *
 **********************************************************************************
 */
void
BezierSubdivide(ZnPoint	*controls,
		ZnReal	t,
		ZnBool	first)
{
  ZnReal	s = 1.0 - t;
  ZnPoint	r[7];
  ZnPoint	a;

  r[0] = controls[0];
  r[6] = controls[3];
  a.x = s*controls[1].x + t*controls[2].x;
  a.y = s*controls[1].y + t*controls[2].y;
  r[1].x = s*r[0].x + t*controls[1].x;
  r[1].y = s*r[0].y + t*controls[1].y;
  r[2].x = s*r[1].x + t*a.x;
  r[2].y = s*r[1].y + t*a.y;
  r[5].x = s*controls[2].x + t*r[6].x;
  r[5].y = s*controls[2].y + t*r[6].y;
  r[4].x = s*a.x + t*r[5].x;
  r[4].y = s*a.y + t*r[5].y;
  r[3].x = s*r[2].x + t*r[4].x;
  r[3].y = s*r[2].y + t*r[4].y;

  if (first) {
    memcpy(controls, r, 4*sizeof(ZnPoint));
  }
  else {
    memcpy(controls, &r[3], 4*sizeof(ZnPoint));
  }      
}


/*
 **********************************************************************************
 *
 * GetBezierPoints --
 *	Use recursive subdivision to approximate the curve. The subdivision stops
 *	when the error is under eps.
 *	This algorithm is adaptive, meaning that it computes the minimum number
 *	of segments needed to render each curve part.
 *
 **********************************************************************************
 */
static void
GetBezierPoints(ZnPoint	*controls,
		ZnList	to_points,
		double	eps)
{
  ZnReal	dist2;
  ZnPoint	mid_segm, mid_cord, delta;

  /*
   * Compute distance between cord center and curve at t = 0.5
   */
  mid_segm.x = (controls[0].x + 3*controls[1].x + 3*controls[2].x + controls[3].x) / 8.0;
  mid_segm.y = (controls[0].y + 3*controls[1].y + 3*controls[2].y + controls[3].y) / 8.0;
  mid_cord.x = (controls[0].x + controls[3].x) / 2.0;
  mid_cord.y = (controls[0].y + controls[3].y) / 2.0;
  delta.x = mid_segm.x - mid_cord.x;
  delta.y = mid_segm.y - mid_cord.y;
  dist2 = delta.x*delta.x + delta.y*delta.y;
  
  if (dist2 > eps) {
    ZnPoint	new_controls[4];
    /*
     * Subdivide the curve at t = 0.5
     * and compute each new curve.
     */
    new_controls[0] = controls[0];
    new_controls[1].x = (controls[0].x + controls[1].x) / 2.0;
    new_controls[1].y = (controls[0].y + controls[1].y) / 2.0;
    new_controls[2].x = (controls[0].x + 2*controls[1].x + controls[2].x) / 4.0;
    new_controls[2].y = (controls[0].y + 2*controls[1].y + controls[2].y) / 4.0;
    new_controls[3] = mid_segm;
    GetBezierPoints(new_controls, to_points, eps);
    
    new_controls[0] = mid_segm;
    new_controls[1].x = (controls[1].x + 2*controls[2].x + controls[3].x) / 4.0;
    new_controls[1].y = (controls[1].y + 2*controls[2].y + controls[3].y) / 4.0;
    new_controls[2].x = (controls[2].x + (controls[3].x)) / 2.0;
    new_controls[2].y = (controls[2].y + (controls[3].y)) / 2.0;
    new_controls[3] = controls[3];
    GetBezierPoints(new_controls, to_points, eps);
  }
  else {
    /*
     * Flat enough add the end to the current path.
     * The start should already be there.
     */
    ZnListAdd(to_points, &controls[3], ZnListTail);
  }
}


/*
 **********************************************************************************
 *
 * GetBezierPath --
 *	Compute in to_points a new set of points describing a Bezier path based
 *	on the control points given in from_points.
 *	If more than four points are given, the algorithm iterate over the
 *	set using the last point of a segment as the first point of the next.
 *	If 3 points are left, they are interpreted as a Bezier segment with
 *	coincident internal control points. If 2 points are left a straight
 *	is emitted.
 *
 **********************************************************************************
 */
void
GetBezierPath(ZnList	from_points,
	      ZnList	to_points)
{
  ZnPoint	*fp;
  int		num_fp, i;
  ZnPoint	s[4];
  
  fp = (ZnPoint *) ZnListArray(from_points);
  num_fp = ZnListSize(from_points);

  /*
   * make sure the output vector is empty, then add the first point.
   */
  ZnListEmpty(to_points);
  ZnListAdd(to_points, &fp[0], ZnListTail);

  for (i = 0; i < num_fp; ) {
    if (i < (num_fp-3)) {
      GetBezierPoints(fp, to_points, 1.0);
      if (i < (num_fp-4)) {
	fp += 3;
	i += 3;
      }
      else {
	break;
      }
    }
    else if (i == (num_fp-3)) {
      s[0] = fp[0];
      s[1] = s[2] = fp[1];
      s[3] = fp[2];
      GetBezierPoints(s, to_points, 1.0);
      break;
    }
    else if (i == (num_fp-2)) {
      ZnListAdd(to_points, &fp[1], ZnListTail);
      break;
    }
  }
}

/*
 **********************************************************************************
 *
 * GetArcCurve --
 *	Compute in to_points a set of Bezier control points describing an arc
 *	path given the start angle, the stop angle and the type: 0 for arc,
 *	1 for chord, 2 for pie slice.
 *	To obtain the actual polygonal shape, the client chould use GetBezierPath
 *	on the returned controls (after applying transform). The returned arc
 *	is circular and centered on 0,0.
 *
 **********************************************************************************
 */
static double arc_nodes_x[4] = { 1.0, 0.0, -1.0, 0.0 };
static double arc_nodes_y[4] = { 0.0, 1.0,  0.0, -1.0 };
static double arc_controls_x[8] = { 1.0, 0.55197, -0.55197, -1.0, -1.0, -0.55197, 0.55197, 1.0 };
static double arc_controls_y[8] = { 0.55197, 1.0, 1.0, 0.55197, -0.55197, -1.0, -1.0, -0.55197 };
void
GetArcPath(ZnReal	start_angle,
	   ZnReal	end_angle,
	   int		type,
	   ZnList	to_points)
{
  int		start_quad, end_quad, quadrant;
  ZnPoint	center_p = { 0.0, 0.0 };
  ZnPoint	start_p = center_p;
  
  /*
   * make sure the output vector is empty.
   */
  ZnListEmpty(to_points);
  
  /*
   * normalize start_angle and end_angle.
   */
  start_angle = fmod(start_angle, 2.0*M_PI);
  if (start_angle < 0.0) {
    start_angle += 2.0*M_PI;
  }
  end_angle = fmod(end_angle, 2.0*M_PI);
  if (end_angle < 0.0) {
    end_angle += 2.0*M_PI;
  }
  if (start_angle >= end_angle) {
    if (start_angle == end_angle) {
      type = 3;
    }
    end_angle += 2.0*M_PI;
  }
  
  /*
   * Now 0 <= start_angle < 2 * M_PI and start_angle <= end_angle.
   */

  start_quad = start_angle / (M_PI / 2.0);
  end_quad = end_angle / (M_PI / 2.0);

  for (quadrant = start_quad; quadrant <= end_quad; quadrant++) {
    ZnPoint controls[4];
    ZnReal  t;
    
    controls[0].x = arc_nodes_x[quadrant % 4];
    controls[0].y = arc_nodes_y[quadrant % 4];
    controls[1].x = arc_controls_x[2 * (quadrant % 4)];
    controls[1].y = arc_controls_y[2 * (quadrant % 4)];
    controls[2].x = arc_controls_x[2 * (quadrant % 4) + 1];
    controls[2].y = arc_controls_y[2 * (quadrant % 4) + 1];
    controls[3].x = arc_nodes_x[(quadrant + 1) % 4];
    controls[3].y = arc_nodes_y[(quadrant + 1) % 4];
    
    if (quadrant == start_quad)	{
      t = Arc2Param(controls, start_angle);
      BezierSubdivide(controls, t, False);
      /*
       * The path is still empty and we have to create the first
       * vertex.
       */
      start_p = controls[0];
      ZnListAdd(to_points, &controls[0], ZnListTail);
    }
    if (quadrant == end_quad) {
      t = Arc2Param(controls, end_angle);
      if (!t) {
	break;
      }
      BezierSubdivide(controls, t, True);
    }
    ZnListAdd(to_points, &controls[1], ZnListTail);
    ZnListAdd(to_points, &controls[2], ZnListTail);
    ZnListAdd(to_points, &controls[3], ZnListTail);
  }

  if (type == 2) {
    ZnListAdd(to_points, &center_p, ZnListTail);
    /*
     * Can't add this one, it would be interpreted by GetBezierPath
     * as an off-curve control. The path should be closed by the client
     * after expansion by GetBezierPath.
     *
      ZnListAdd(to_points, &start_p, ZnListTail); 
     */
  }
  else if (type == 1) {
    ZnListAdd(to_points, &start_p, ZnListTail); 
  }
}


/*
 **********************************************************************************
 *
 * SmoothPathWithBezier --
 *	Compute in to_points a new set of points describing a smoothed path based
 *	on the path given in from_points. The algorithm use Bezier cubic curves.
 *
 **********************************************************************************
 */
void
SmoothPathWithBezier(ZnPoint	*fp,
		     int	num_fp,
		     ZnList	to_points)
{
  ZnBool	closed;
  ZnPoint	s[4];
  int		i;

  /*
   * make sure the output vector is empty
   */
  ZnListEmpty(to_points);

  /*
   * If the curve is closed, generates a Bezier curve that
   * spans the closure. Else simply add the first point to
   * the path.
   */
  if ((fp[0].x == fp[num_fp-1].x) && (fp[0].y == fp[num_fp-1].y)) {
    closed = True;
    s[0].x = REAL_TO_INT(0.5*fp[num_fp-2].x + 0.5*fp[0].x);
    s[0].y = REAL_TO_INT(0.5*fp[num_fp-2].y + 0.5*fp[0].y);
    s[1].x = REAL_TO_INT(0.167*fp[num_fp-2].x + 0.833*fp[0].x);
    s[1].y = REAL_TO_INT(0.167*fp[num_fp-2].y + 0.833*fp[0].y);
    s[2].x = REAL_TO_INT(0.833*fp[0].x + 0.167*fp[1].x);
    s[2].y = REAL_TO_INT(0.833*fp[0].y + 0.167*fp[1].y);
    s[3].x = REAL_TO_INT(0.5*fp[0].x + 0.5*fp[1].x);
    s[3].y = REAL_TO_INT(0.5*fp[0].y + 0.5*fp[1].y);
    ZnListAdd(to_points, s, ZnListTail);
    GetBezierPoints(s, to_points, 1.0);
  }
  else {
    closed = False;
    ZnListAdd(to_points, &fp[0], ZnListTail);
  }

  for (i = 2; i < num_fp; i++, fp++) {
    /*
     * Setup the first two control points. This differ
     * for first segment of open curves.
     */
    if ((i == 2) && !closed) {
      s[0] = fp[0];
      s[1].x = REAL_TO_INT(0.333*fp[0].x + 0.667*fp[1].x);
      s[1].y = REAL_TO_INT(0.333*fp[0].y + 0.667*fp[1].y);
    }
    else {
      s[0].x = REAL_TO_INT(0.5*fp[0].x + 0.5*fp[1].x);
      s[0].y = REAL_TO_INT(0.5*fp[0].y + 0.5*fp[1].y);
      s[1].x = REAL_TO_INT(0.167*fp[0].x + 0.833*fp[1].x);
      s[1].y = REAL_TO_INT(0.167*fp[0].y + 0.833*fp[1].y);
    }

    /*
     * Setup the last two control points. This also differ
     * for last segment of open curves.
     */
    if ((i == num_fp-1) && !closed) {
      s[2].x = REAL_TO_INT(0.667*fp[1].x + 0.333*fp[2].x);
      s[2].y = REAL_TO_INT(0.667*fp[1].y + 0.333*fp[2].y);
      s[3] = fp[2];
    }
    else {
      s[2].x = REAL_TO_INT(0.833*fp[1].x + 0.167*fp[2].x);
      s[2].y = REAL_TO_INT(0.833*fp[1].y + 0.167*fp[2].y);
      s[3].x = REAL_TO_INT(0.5*fp[1].x + 0.5*fp[2].x);
      s[3].y = REAL_TO_INT(0.5*fp[1].y + 0.5*fp[2].y);
    }

    /*
     * If the first two points or the last two are equal
     * output the last control point. Else generate the
     * Bezier curve.
     */
    if (((fp[0].x == fp[1].x) && (fp[0].y == fp[1].y)) ||
	((fp[1].x == fp[2].x) && (fp[1].y == fp[2].y))) {
      ZnListAdd(to_points, &s[3], ZnListTail);
    }
    else {
      GetBezierPoints(s, to_points, 1.0);
    }
  }
}


/*
 **********************************************************************************
 *
 * 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.
 *	If bbox is non null, it is filled with the bounding box of the end.
 *
 * For the time being this procedure handles open/filled arrows.
 *
 * Here are the parameters describing arrows.
 *
 *              *  | ARROW_SHAPE_C
 *            **   |
 *          * ***************************
 *        *  *
 *      *   *  +p1                      +p2
 *      | * |*
 *      |   * ***************************
 *      |   | **
 *      |   |   *
 *      |   |   |
 *      |---|   | ARROW_SHAPE_A
 *      |       |
 *      |-------| ARROW_SHAPE_B
 *
 **********************************************************************************
 */
void
GetLineEnd(ZnPoint	*p1,
	   ZnPoint	*p2,
	   unsigned int	line_width,
	   int		cap_style,
	   ZnLineEnd	end_style,
	   ZnPoint	*points)
{
  ZnReal	dx, dy, length, temp, backup;
  ZnReal	frac_height, sin_theta, cos_theta;
  ZnReal	vert_x, vert_y;
  ZnReal	shape_a, shape_b, shape_c;
  
  if (end_style != NULL) {
    shape_a = end_style->shape_a + 0.001;
    shape_b = end_style->shape_b + 0.001;
    shape_c = end_style->shape_c + line_width/2.0 + 0.001;

    frac_height = (line_width/2.0) / shape_c;
    dx = p1->x - p2->x;
    dy = p1->y - p2->y;
    length = hypot(dx, dy);
    if (length == 0) {
      sin_theta = cos_theta = 0.0;
    }
    else {
      sin_theta = dy/length;
      cos_theta = dx/length;
    }

    if (cap_style != CapProjecting) {
      temp = frac_height;
    }
    else {
      temp = line_width / shape_c;
    }
    backup = temp * shape_b + shape_a * (1.0 - temp) / 2.0;
    points[0].x = points[5].x = p1->x + backup * cos_theta;
    points[0].y = points[5].y = p1->y + backup * sin_theta;
    
    vert_x = points[0].x - shape_a*cos_theta;
    vert_y = points[0].y - shape_a*sin_theta;
    temp = shape_c*sin_theta;
    points[1].x = REAL_TO_INT(points[0].x - shape_b*cos_theta + temp);
    points[4].x = REAL_TO_INT(points[1].x - 2*temp);
    temp = shape_c*cos_theta;
    points[1].y = REAL_TO_INT(points[0].y - shape_b*sin_theta - temp);
    points[4].y = REAL_TO_INT(points[1].y + 2*temp);
    points[2].x = REAL_TO_INT(points[1].x*frac_height + vert_x*(1.0-frac_height));
    points[2].y = REAL_TO_INT(points[1].y*frac_height + vert_y*(1.0-frac_height));
    points[3].x = REAL_TO_INT(points[4].x*frac_height + vert_x*(1.0-frac_height));
    points[3].y = REAL_TO_INT(points[4].y*frac_height + vert_y*(1.0-frac_height));
  }
}