src/main/xspline.c

/*
 * Source code from Xfig 3.2.4 modified to work with arrays of doubles
 * instead linked lists of F_points and to remove some globals(!)
 * See copyright etc below.
 *
 * #included from engine.c
 * That manages the R_alloc stack.
 */

/*
 * From w_drawprim.h
 */

#define         MAXNUMPTS       25000

/*
 * From u_draw.c
 */

/*
 * FIG : Facility for Interactive Generation of figures
 * Copyright (c) 1985-1988 by Supoj Sutanthavibul
 * Parts Copyright (c) 1989-2002 by Brian V. Smith
 * Parts Copyright (c) 1991 by Paul King
 * Parts Copyright (c) 1992 by James Tough
 * Parts Copyright (c) 1998 by Georg Stemmer
 * Parts Copyright (c) 1995 by C. Blanc and C. Schlick
 *
 * Any party obtaining a copy of these files is granted, free of charge, a
 * full and unrestricted irrevocable, world-wide, paid up, royalty-free,
 * nonexclusive right and license to deal in this software and
 * documentation files (the "Software"), including without limitation the
 * rights to use, copy, modify, merge, publish and/or distribute copies of
 * the Software, and to permit persons who receive copies from any such
 * party to do so, with the only requirement being that this copyright
 * notice remain intact.
 *
 */

/************** POLYGON/CURVE DRAWING FACILITIES ****************/

static int npoints;
static int max_points;
static double *xpoints;
static double *ypoints;

/************* Code begins here *************/

/* R_allocs or mallocs global arrays */
static Rboolean
add_point(double x, double y, pGEDevDesc dd)
{
    if (npoints >= max_points) {
	int tmp_n;
	double *tmp_px;
	double *tmp_py;
	tmp_n = max_points + 200;
	/* too many points, return false */
	if (tmp_n > MAXNUMPTS) {
	    error(_("add_point - reached MAXNUMPTS (%d)"),tmp_n);
	}
	if (max_points == 0) {
	    tmp_px = (double *) R_alloc(tmp_n, sizeof(double));
	    tmp_py = (double *) R_alloc(tmp_n, sizeof(double));
	} else {
	    tmp_px = (double *) S_realloc((char *) xpoints,
					  tmp_n, max_points,
					  sizeof(double));
	    tmp_py = (double *) S_realloc((char *) ypoints,
					  tmp_n, max_points,
					  sizeof(double));
	}
	if (tmp_px == NULL || tmp_py == NULL) {
	    error(_("insufficient memory to allocate point array"));
	}
	xpoints = tmp_px;
	ypoints = tmp_py;
	max_points = tmp_n;
    }
    /* ignore identical points */
    if (npoints > 0 && xpoints[npoints-1] == x && ypoints[npoints-1] == y)
	return TRUE;
    /*
     * Convert back from 1200ppi to DEVICE coordinates
     */
    xpoints[npoints] = toDeviceX(x / 1200, GE_INCHES, dd);
    ypoints[npoints] = toDeviceY(y / 1200, GE_INCHES, dd);
    npoints = npoints + 1;
    return TRUE;
}

/*
 * From u_draw_spline.c
 */

/*
 * FIG : Facility for Interactive Generation of figures
 * Copyright (c) 1985-1988 by Supoj Sutanthavibul
 * Parts Copyright (c) 1989-2002 by Brian V. Smith
 * Parts Copyright (c) 1991 by Paul King
 *
 * Any party obtaining a copy of these files is granted, free of charge, a
 * full and unrestricted irrevocable, world-wide, paid up, royalty-free,
 * nonexclusive right and license to deal in this software and
 * documentation files (the "Software"), including without limitation the
 * rights to use, copy, modify, merge, publish and/or distribute copies of
 * the Software, and to permit persons who receive copies from any such
 * party to do so, with the only requirement being that this copyright
 * notice remain intact.
 *
 */

/* THIS FILE IS #included FROM u_draw.c and u_geom.c */

/********************* CURVES FOR SPLINES *****************************

 The following spline drawing routines are from

    "X-splines : A Spline Model Designed for the End User"

    by Carole BLANC and Christophe SCHLICK,
    in Proceedings of SIGGRAPH ' 95

***********************************************************************/

#define         HIGH_PRECISION    0.5
#define         LOW_PRECISION     1.0
#define         ZOOM_PRECISION    5.0
#define         ARROW_START       4
#define         MAX_SPLINE_STEP   0.2

/***********************************************************************/

#define Q(s)  (-(s))
#define EQN_NUMERATOR(dim) \
  (A_blend[0]*dim[0]+A_blend[1]*dim[1]+A_blend[2]*dim[2]+A_blend[3]*dim[3])

static double
f_blend(double numerator, double denominator)
{
  double p = 2 * denominator * denominator;
  double u = numerator / denominator;
  double u2 = u * u;

  return (u * u2 * (10 - p + (2*p - 15)*u + (6 - p)*u2));
}

static double
g_blend(double u, double q)             /* p equals 2 */
{
  return(u*(q + u*(2*q + u*(8 - 12*q + u*(14*q - 11 + u*(4 - 5*q))))));
}

static double
h_blend(double u, double q)
{
    double u2=u*u;
    return (u * (q + u * (2 * q + u2 * (-2*q - u*q))));
}

static void
negative_s1_influence(double t, double s1, double *A0, double *A2)
{
  *A0 = h_blend(-t, Q(s1));
  *A2 = g_blend(t, Q(s1));
}

static void
negative_s2_influence(double t, double s2, double *A1, double *A3)
{
  *A1 = g_blend(1-t, Q(s2));
  *A3 = h_blend(t-1, Q(s2));
}

static void
positive_s1_influence(double k, double t, double s1, double *A0, double *A2)
{
  double Tk;

  Tk = k+1+s1;
  *A0 = (t+k+1<Tk) ? f_blend(t+k+1-Tk, k-Tk) : 0.0;

  Tk = k+1-s1;
  *A2 = f_blend(t+k+1-Tk, k+2-Tk);
}

static void
positive_s2_influence(double k, double t, double s2, double *A1, double *A3)
{
  double Tk;

  Tk = k+2+s2;
  *A1 = f_blend(t+k+1-Tk, k+1-Tk);

  Tk = k+2-s2;
  *A3 = (t+k+1>Tk) ? f_blend(t+k+1-Tk, k+3-Tk) : 0.0;
}

static void
point_adding(double *A_blend, double *px, double *py,
	     pGEDevDesc dd)
{
  double weights_sum;

  weights_sum = A_blend[0] + A_blend[1] + A_blend[2] + A_blend[3];
  add_point(EQN_NUMERATOR(px) / (weights_sum),
	    EQN_NUMERATOR(py) / (weights_sum),
	    dd);
}

static void
point_computing(double *A_blend,
		double *px, double *py,
		double *x, double *y)
{
  double weights_sum;

  weights_sum = A_blend[0] + A_blend[1] + A_blend[2] + A_blend[3];

  *x = EQN_NUMERATOR(px) / (weights_sum);
  *y = EQN_NUMERATOR(py) / (weights_sum);
}

static double
step_computing(int k,
	       double *px, double *py,
	       double s1, double s2,
	       double precision,
               pGEDevDesc dd)
{
  double A_blend[4];
  double xstart, ystart, xend, yend, xmid, ymid, xlength, ylength;
  double start_to_end_dist, devWidth, devHeight, devDiag, number_of_steps;
  double  step, angle_cos, scal_prod, xv1, xv2, yv1, yv2, sides_length_prod;

  /* This function computes the step used to draw the segment (p1, p2)
     (xv1, yv1) : coordinates of the vector from middle to origin
     (xv2, yv2) : coordinates of the vector from middle to extremity */

  if ((s1 == 0) && (s2 == 0))
    return(1.0);              /* only one step in case of linear segment */

  /* compute coordinates of the origin */
  if (s1>0) {
      if (s2<0) {
	  positive_s1_influence(k, 0.0, s1, &A_blend[0], &A_blend[2]);
	  negative_s2_influence(0.0, s2, &A_blend[1], &A_blend[3]);
      } else {
	  positive_s1_influence(k, 0.0, s1, &A_blend[0], &A_blend[2]);
	  positive_s2_influence(k, 0.0, s2, &A_blend[1], &A_blend[3]);
      }
      point_computing(A_blend, px, py, &xstart, &ystart);
  } else {
      xstart = px[1];
      ystart = py[1];
  }

  /* compute coordinates  of the extremity */
  if (s2>0) {
      if (s1<0) {
	  negative_s1_influence(1.0, s1, &A_blend[0], &A_blend[2]);
	  positive_s2_influence(k, 1.0, s2, &A_blend[1], &A_blend[3]);
      } else {
	  positive_s1_influence(k, 1.0, s1, &A_blend[0], &A_blend[2]);
	  positive_s2_influence(k, 1.0, s2, &A_blend[1], &A_blend[3]);
      }
      point_computing(A_blend, px, py, &xend, &yend);
  } else {
      xend = px[2];
      yend = py[2];
  }

  /* compute coordinates  of the middle */
  if (s2>0) {
      if (s1<0) {
	  negative_s1_influence(0.5, s1, &A_blend[0], &A_blend[2]);
	  positive_s2_influence(k, 0.5, s2, &A_blend[1], &A_blend[3]);
      } else {
	  positive_s1_influence(k, 0.5, s1, &A_blend[0], &A_blend[2]);
	  positive_s2_influence(k, 0.5, s2, &A_blend[1], &A_blend[3]);
	}
  } else if (s1<0) {
      negative_s1_influence(0.5, s1, &A_blend[0], &A_blend[2]);
      negative_s2_influence(0.5, s2, &A_blend[1], &A_blend[3]);
  } else {
      positive_s1_influence(k, 0.5, s1, &A_blend[0], &A_blend[2]);
      negative_s2_influence(0.5, s2, &A_blend[1], &A_blend[3]);
  }
  point_computing(A_blend, px, py, &xmid, &ymid);

  xv1 = xstart - xmid;
  yv1 = ystart - ymid;
  xv2 = xend - xmid;
  yv2 = yend - ymid;

  scal_prod = xv1*xv2 + yv1*yv2;

  sides_length_prod = sqrt((xv1*xv1 + yv1*yv1)*(xv2*xv2 + yv2*yv2));

  /* compute cosinus of origin-middle-extremity angle, which approximates the
     curve of the spline segment */
  if (sides_length_prod == 0.0)
    angle_cos = 0.0;
  else
    angle_cos = scal_prod/sides_length_prod;

  xlength = xend - xstart;
  ylength = yend - ystart;

  start_to_end_dist = sqrt(xlength*xlength + ylength*ylength);

  /* Paul 2009-01-25
   * It is possible for origin and extremity to be very remote
   * indeed (if the control points are located WAY off the device).
   * In order to avoid having ridiculously many steps, limit
   * the start_to_end_dist to being the length of the diagonal of the 
   * device.
   */
  devWidth = fromDeviceWidth(toDeviceWidth(1, GE_NDC, dd),
                             GE_INCHES, dd)*1200;
  devHeight = fromDeviceHeight(toDeviceHeight(1, GE_NDC, dd),
                               GE_INCHES, dd)*1200;
  devDiag = sqrt(devWidth* devWidth + devHeight*devHeight);
  if (start_to_end_dist > devDiag)
      start_to_end_dist = devDiag;

  /* more steps if segment's origin and extremity are remote */
  number_of_steps = sqrt(start_to_end_dist)/2;

  /* more steps if the curve is high */
  number_of_steps += (int)((1 + angle_cos)*10);

  if (number_of_steps == 0)
    step = 1;
  else
    step = precision/number_of_steps;

  if ((step > MAX_SPLINE_STEP) || (step == 0))
    step = MAX_SPLINE_STEP;

  return (step);
}

static void
spline_segment_computing(double step, int k,
			 double *px, double *py,
			 double s1, double s2,
			 pGEDevDesc dd)
{
  double A_blend[4];
  double t;

  if (s1<0) {
     if (s2<0) {
	 for (t=0.0 ; t<1 ; t+=step) {
	     negative_s1_influence(t, s1, &A_blend[0], &A_blend[2]);
	     negative_s2_influence(t, s2, &A_blend[1], &A_blend[3]);

	     point_adding(A_blend, px, py, dd);
	 }
     } else {
	 for (t=0.0 ; t<1 ; t+=step) {
	     negative_s1_influence(t, s1, &A_blend[0], &A_blend[2]);
	     positive_s2_influence(k, t, s2, &A_blend[1], &A_blend[3]);

	     point_adding(A_blend, px, py, dd);
	 }
     }
  } else if (s2<0) {
      for (t=0.0 ; t<1 ; t+=step) {
	     positive_s1_influence(k, t, s1, &A_blend[0], &A_blend[2]);
	     negative_s2_influence(t, s2, &A_blend[1], &A_blend[3]);

	     point_adding(A_blend, px, py, dd);
      }
  } else {
      for (t=0.0 ; t<1 ; t+=step) {
	     positive_s1_influence(k, t, s1, &A_blend[0], &A_blend[2]);
	     positive_s2_influence(k, t, s2, &A_blend[1], &A_blend[3]);

	     point_adding(A_blend, px, py, dd);
      }
  }
}

/*
 * For adding last line segment when computing open spline
 * WITHOUT end control points repeated
 * (i.e., can't just connect to last control point)
 */
static void
spline_last_segment_computing(double step, int k,
			      double *px, double *py,
			      double s1, double s2,
			      pGEDevDesc dd)
{
  double A_blend[4];
  double t = 1;

  if (s1<0) {
      if (s2<0) {
	  negative_s1_influence(t, s1, &A_blend[0], &A_blend[2]);
	  negative_s2_influence(t, s2, &A_blend[1], &A_blend[3]);

	  point_adding(A_blend, px, py, dd);
      } else {
	  negative_s1_influence(t, s1, &A_blend[0], &A_blend[2]);
	  positive_s2_influence(k, t, s2, &A_blend[1], &A_blend[3]);

	  point_adding(A_blend, px, py, dd);
      }
  } else if (s2<0) {
      positive_s1_influence(k, t, s1, &A_blend[0], &A_blend[2]);
      negative_s2_influence(t, s2, &A_blend[1], &A_blend[3]);

      point_adding(A_blend, px, py, dd);
  } else {
      positive_s1_influence(k, t, s1, &A_blend[0], &A_blend[2]);
      positive_s2_influence(k, t, s2, &A_blend[1], &A_blend[3]);

      point_adding(A_blend, px, py, dd);
  }
}

/********************* MAIN METHODS *************************************/

/*
 * x and y are in DEVICE coordinates
 * xfig works in 1200ppi
 *   (http://www.csit.fsu.edu/~burkardt/data/fig/fig_format.html)
 * so convert to 1200ppi so that step calculations are correct
 */
#define COPY_CONTROL_POINT(PI, I, N) \
      px[PI] = fromDeviceX(x[(I) % N], GE_INCHES, dd) * 1200; \
      py[PI] = fromDeviceY(y[(I) % N], GE_INCHES, dd) * 1200; \
      ps[PI] = s[(I) % N]

#define NEXT_CONTROL_POINTS(K, N) \
      COPY_CONTROL_POINT(0, K, N); \
      COPY_CONTROL_POINT(1, K + 1, N); \
      COPY_CONTROL_POINT(2, K + 2, N); \
      COPY_CONTROL_POINT(3, K + 3, N)

#define INIT_CONTROL_POINTS(N) \
      COPY_CONTROL_POINT(0, N - 1, N); \
      COPY_CONTROL_POINT(1, 0, N); \
      COPY_CONTROL_POINT(2, 1, N); \
      COPY_CONTROL_POINT(3, 2, N)

#define SPLINE_SEGMENT_LOOP(K, PX, PY, S1, S2, PREC) \
      step = step_computing(K, PX, PY, S1, S2, PREC, dd);    \
      spline_segment_computing(step, K, PX, PY, S1, S2, dd)

static Rboolean
compute_open_spline(int n, double *x, double *y, double *s,
		    Rboolean repEnds,
		    double precision,
		    pGEDevDesc dd)
{
  int       k;
  double     step = 0.0 /* -Wall */;
  double px[4];
  double py[4];
  double ps[4]={0.,0.,0.,0.};

  max_points = 0;
  npoints = 0;
  xpoints = NULL;
  ypoints = NULL;

  if (repEnds && n < 2)
      error(_("there must be at least two control points"));
  if (!repEnds && n < 4)
      error(_("there must be at least four control points"));

  if (repEnds) {
      /* first control point is needed twice for the first segment */
      COPY_CONTROL_POINT(0, 0, n);
      COPY_CONTROL_POINT(1, 0, n);
      COPY_CONTROL_POINT(2, 1, n);

      if (n == 2) {
	COPY_CONTROL_POINT(3, 1, n);
      } else {
	COPY_CONTROL_POINT(3, 2, n);
      }

      for (k = 0 ; ; k++) {
	  SPLINE_SEGMENT_LOOP(k, px, py, ps[1], ps[2], precision);
          /* Paul 2008-12-05
           * >= rather than == to handle special case of n == 2
           */
	  if (k + 3 >= n)
	      break;
	  NEXT_CONTROL_POINTS(k, n);
      }

      /* last control point is needed twice for the last segment */
      COPY_CONTROL_POINT(0, n - 3, n);
      COPY_CONTROL_POINT(1, n - 2, n);
      COPY_CONTROL_POINT(2, n - 1, n);
      COPY_CONTROL_POINT(3, n - 1, n);
      SPLINE_SEGMENT_LOOP(k, px, py, ps[1], ps[2], precision);

      add_point(px[3], py[3], dd);
  } else {
      for (k = 0 ; k + 3 < n ; k++) {
	  NEXT_CONTROL_POINTS(k, n);
	  SPLINE_SEGMENT_LOOP(k, px, py, ps[1], ps[2], precision);
      }
      spline_last_segment_computing(step, n - 4, px, py, ps[1], ps[2], dd);
  }

  return TRUE;
}

static Rboolean
compute_closed_spline(int n, double *x, double *y, double *s,
		      double precision,
		      pGEDevDesc dd)
{
  int k;
  double step;
  double px[4];
  double py[4];
  double ps[4];

  max_points = 0;
  npoints = 0;
  xpoints = NULL;
  ypoints = NULL;

  if (n < 3)
      error(_("There must be at least three control points"));

  INIT_CONTROL_POINTS(n);

  for (k = 0 ; k < n ; k++) {
      SPLINE_SEGMENT_LOOP(k, px, py, ps[1], ps[2], precision);
      NEXT_CONTROL_POINTS(k, n);
  }

  return TRUE;
}