1/***********************************************************************
2 *                                                                     *
3 * $Id: hpgsbezier.c 270 2006-01-29 21:12:23Z softadm $
4 *                                                                     *
5 * hpgs - HPGl Script, a hpgl/2 interpreter, which uses a Postscript   *
6 *        API for rendering a scene and thus renders to a variety of   *
7 *        devices and fileformats.                                     *
8 *                                                                     *
9 * (C) 2004-2006 ev-i Informationstechnologie GmbH  http://www.ev-i.at *
10 *                                                                     *
11 * Author: Wolfgang Glas                                               *
12 *                                                                     *
13 *  hpgs is free software; you can redistribute it and/or              *
14 * modify it under the terms of the GNU Lesser General Public          *
15 * License as published by the Free Software Foundation; either        *
16 * version 2.1 of the License, or (at your option) any later version.  *
17 *                                                                     *
18 * hpgs is distributed in the hope that it will be useful,             *
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of      *
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU   *
21 * Lesser General Public License for more details.                     *
22 *                                                                     *
23 * You should have received a copy of the GNU Lesser General Public    *
24 * License along with this library; if not, write to the               *
25 * Free Software  Foundation, Inc., 59 Temple Place, Suite 330,        *
26 * Boston, MA  02111-1307  USA                                         *
27 *                                                                     *
28 ***********************************************************************
29 *                                                                     *
30 * Bezier spline length and scanline cut implementation.               *
31 *                                                                     *
32 ***********************************************************************/
33
34#include <hpgspaint.h>
35#include <math.h>
36#include <assert.h>
37
38// #define HPGS_BEZIER_DEBUG
39
40/*! Calculates the derivation of the bezier curve starting
41    at point \c i of \c path at curve parameter \c t.
42*/
43void hpgs_bezier_path_delta(const hpgs_paint_path *path, int i,
44			    double t, hpgs_point *p)
45{
46  double tp,tm;
47
48  tp = t + 0.5;
49  tm = 0.5 - t;
50
51  p->x = 3.0* ((path->points[i+1].p.x-path->points[i].p.x)*tm*tm+
52	       2.0 * (path->points[i+2].p.x-path->points[i+1].p.x) * tm*tp+
53	       (path->points[i+3].p.x-path->points[i+2].p.x) * tp*tp);
54
55  p->y = 3.0 * ((path->points[i+1].p.y-path->points[i].p.y)*tm*tm+
56		2.0 * (path->points[i+2].p.y-path->points[i+1].p.y) * tm*tp+
57		(path->points[i+3].p.y-path->points[i+2].p.y) * tp*tp);
58}
59
60/*! Calculates the derivation of the bezier curve starting
61    at point \c i of \c path at curve parameter \c t.
62*/
63double hpgs_bezier_path_delta_x(const hpgs_paint_path *path, int i, double t)
64{
65  double tp,tm;
66
67  tp = t + 0.5;
68  tm = 0.5 - t;
69
70  return 3.0* ((path->points[i+1].p.x-path->points[i].p.x)*tm*tm+
71	       2.0 * (path->points[i+2].p.x-path->points[i+1].p.x) * tm*tp+
72	       (path->points[i+3].p.x-path->points[i+2].p.x) * tp*tp);
73}
74
75/*! Calculates the derivation of the bezier curve starting
76    at point \c i of \c path at curve parameter \c t.
77*/
78double hpgs_bezier_path_delta_y(const hpgs_paint_path *path, int i, double t)
79{
80  double tp,tm;
81
82  tp = t + 0.5;
83  tm = 0.5 - t;
84
85  return 3.0* ((path->points[i+1].p.y-path->points[i].p.y)*tm*tm+
86	       2.0 * (path->points[i+2].p.y-path->points[i+1].p.y) * tm*tp+
87	       (path->points[i+3].p.y-path->points[i+2].p.y) * tp*tp);
88}
89
90/*! Calculates the tangent to the bezier curve starting
91    at point \c i of \c path at curve parameter \c t.
92
93    If \c orientation is \c 1, the tangent is calculated with increasing
94    curve parameter, if we meet a singularity (right tangent).
95
96    If \c orientation is \c -1, the tangent is calculated with decreasing
97    curve parameter, if we meet a singularity (left tangent).
98
99    If \c orientation is \c 0, the tangent is calculated with decreased
100    and increased curve parameter, if we meet a singularity (mid tangent).
101*/
102
103void hpgs_bezier_path_tangent(const hpgs_paint_path *path, int i,
104                              double t, int orientation,
105                              double ytol, hpgs_point *p)
106{
107  hpgs_bezier_path_delta(path,i,t,p);
108
109  double l = hypot(p->x,p->y);
110
111  if (l < ytol)
112    {
113      hpgs_point p1;
114
115      if (orientation <= 0)
116        hpgs_bezier_path_point(path,i,t-1.0e-4,&p1);
117      else
118        hpgs_bezier_path_point(path,i,t,&p1);
119
120      if (orientation >= 0)
121        hpgs_bezier_path_point(path,i,t+1.0e-4,p);
122      else
123        hpgs_bezier_path_point(path,i,t,p);
124
125      p->x -= p1.x;
126      p->y -= p1.y;
127
128      l = hypot(p->x,p->y);
129      if (l < ytol*1.0e-4) return;
130    }
131
132  p->x /= l;
133  p->y /= l;
134}
135
136static int quad_roots (double a, double b, double c,
137                       double t0, double t1, double *t)
138{
139  if (fabs(a) < (fabs(b)>fabs(c) ? fabs(b) : fabs(c)) * 1.0e-8 || fabs(a) < 1.0e-16)
140    {
141      if (fabs(b) < fabs(c) * 1.0e-8)
142        return 0;
143
144      t[0] = -c/b;
145      return t[0] > t0 && t[0] < t1;
146    }
147
148  double p = b/a;
149  double q = c/a;
150
151  double det = 0.25*p-q;
152
153  if (det < 0) return 0;
154
155  det = sqrt(det);
156
157  if (p > 0.0)
158    t[0] = -0.5*p - det;
159  else
160    t[0] = -0.5*p + det;
161
162  if (fabs(t[0]) < 1.0e-8)
163    return t[0] > t0 && t[0] < t1;
164
165  t[1] = q / t[0];
166
167  if (t[0] <= t0 || t[0] >= t1)
168    {
169      t[0] = t[1];
170      return t[0] > t0 && t[0] < t1;
171    }
172
173  if (t[1] <= t0 || t[1] >= t1)
174    return 1;
175
176  if (t[0] > t[1])
177    {
178      double tmp = t[1];
179      t[1] = t[0];
180      t[0] = tmp;
181    }
182
183  return fabs(t[0]-t[1]) < 1.0e-8 ? 1 : 2;
184}
185
186
187/*!
188    Calculates the singularities of the bezier curve in the
189    parameter interval from \c t0 to \c t1.
190
191    We calculate points, where the cross product of the first derivative
192    and the second derivative vanishes. These points include points, where
193    the first or second derivative vanishes. Additionally, turning points
194    of the curve are also included.
195
196    \c tx is a pointer to an array of at least two double values,
197    in which the position of the singularities is returned.
198 */
199void hpgs_bezier_path_singularities(const hpgs_paint_path *path, int i,
200                                    double t0, double t1,
201                                    int *nx, double *tx)
202{
203  // The second order spline p'(t) cross p''(t)
204  double k0 =
205    (path->points[i+1].p.x-path->points[i].p.x)*(path->points[i+2].p.y-path->points[i+1].p.y) -
206    (path->points[i+1].p.y-path->points[i].p.y)*(path->points[i+2].p.x-path->points[i+1].p.x);
207
208  double k1 =
209    (path->points[i+1].p.x-path->points[i].p.x)*(path->points[i+3].p.y-path->points[i+2].p.y) -
210    (path->points[i+1].p.y-path->points[i].p.y)*(path->points[i+3].p.x-path->points[i+2].p.x);
211
212  double k2 =
213    (path->points[i+2].p.x-path->points[i+1].p.x)*(path->points[i+3].p.y-path->points[i+2].p.y) -
214    (path->points[i+2].p.y-path->points[i+1].p.y)*(path->points[i+3].p.x-path->points[i+2].p.x);
215
216
217  // coefficients for the quadratic equations.
218  double a = k0 - k1 + k2;
219
220  double b = k2 - k0;
221
222  double c = 0.25 * (k0 + k1 + k2);
223
224  *nx = quad_roots (a,b,c,t0,t1,tx);
225}
226
227static void add_quad (const hpgs_point *p1,
228                      const hpgs_point *d1,
229                      const hpgs_point *p2,
230                      const hpgs_point *d2,
231                      int *nx, hpgs_point *points)
232{
233  double det = d1->y * d2->x - d1->x * d2->y;
234
235  if (fabs(det) < 1.0e-8)
236    {
237      points[*nx].x = 0.5 * (p1->x + p2->x);
238      points[*nx].y = 0.5 * (p1->y + p2->y);
239      ++*nx;
240    }
241  else
242    {
243      double s = (d2->x*(p2->y-p1->y)-
244                  d2->y*(p2->x-p1->x)  )/det;
245
246      points[*nx].x = p1->x + s*d1->x;
247      points[*nx].y = p1->y + s*d1->y;
248      ++*nx;
249    }
250
251  points[*nx].x = p2->x;
252  points[*nx].y = p2->y;
253  ++*nx;
254}
255
256/*!
257    Approximates the bezier curve segment in the
258    parameter interval from \c t0 to \c t1 with quadratic bezier splines.
259
260    The array \c points is used to return the quadratic bezier segments and
261    must have a dimension of at least 16 points.
262
263    \c nx returns the number of generated points, which is always an even
264    number. The first point is the control point of the first quadratic
265    spline, the second point the endpoint of the first quadratic spline
266    and so on.
267 */
268void hpgs_bezier_path_to_quadratic(const hpgs_paint_path *path, int i,
269                                   double t0, double t1,
270                                   int *nx, hpgs_point *points)
271{
272  // The second order spline p'(t) cross p''(t)
273  double k0 =
274    (path->points[i+1].p.x-path->points[i].p.x)*(path->points[i+2].p.y-path->points[i+1].p.y) -
275    (path->points[i+1].p.y-path->points[i].p.y)*(path->points[i+2].p.x-path->points[i+1].p.x);
276
277  double k1 =
278    (path->points[i+1].p.x-path->points[i].p.x)*(path->points[i+3].p.y-path->points[i+2].p.y) -
279    (path->points[i+1].p.y-path->points[i].p.y)*(path->points[i+3].p.x-path->points[i+2].p.x);
280
281  double k2 =
282    (path->points[i+2].p.x-path->points[i+1].p.x)*(path->points[i+3].p.y-path->points[i+2].p.y) -
283    (path->points[i+2].p.y-path->points[i+1].p.y)*(path->points[i+3].p.x-path->points[i+2].p.x);
284
285
286  double xmax,xmin,ymin,ymax;
287
288  // coefficients for the quadratic equations.
289  double a = k0 - k1 + k2;
290
291  double b = k2 - k0;
292
293  double c = 0.25 * (k0 + k1 + k2);
294
295  // at most three partition points:
296  // extrem, roots of the curvature spline
297  double tpart[5];
298
299  // get the roots.
300  int ipart;
301  int npart = quad_roots (a,b,c,t0,t1,tpart+1) + 1;
302
303  // add the extreme point.
304  if (fabs(a) > fabs(b))
305    {
306      double text = -0.5 * b/a;
307
308      if (npart == 3)
309        {
310          tpart[3] = tpart[2];
311          tpart[2] = text;
312          npart = 3;
313        }
314      else if (text > t0 && text < t1)
315        {
316          if (npart == 2 && text < tpart[1])
317            {
318              tpart[2] = tpart[1];
319              tpart[1] = text;
320            }
321          else
322            tpart[npart] = text;
323
324          ++npart;
325        }
326    }
327
328  tpart[0] = t0;
329  tpart[npart] = t1;
330
331  xmax = xmin = path->points[i].p.x;
332  if (path->points[i+1].p.x > xmax) xmax = path->points[i+1].p.x;
333  if (path->points[i+2].p.x > xmax) xmax = path->points[i+2].p.x;
334  if (path->points[i+3].p.x > xmax) xmax = path->points[i+3].p.x;
335  if (path->points[i+1].p.x < xmin) xmin = path->points[i+1].p.x;
336  if (path->points[i+2].p.x < xmin) xmin = path->points[i+2].p.x;
337  if (path->points[i+3].p.x < xmin) xmin = path->points[i+3].p.x;
338
339  ymax = ymin = path->points[i].p.y;
340  if (path->points[i+1].p.y > ymax) ymax = path->points[i+1].p.y;
341  if (path->points[i+2].p.y > ymax) ymax = path->points[i+2].p.y;
342  if (path->points[i+3].p.y > ymax) ymax = path->points[i+3].p.y;
343  if (path->points[i+1].p.y < ymin) ymin = path->points[i+1].p.y;
344  if (path->points[i+2].p.y < ymin) ymin = path->points[i+2].p.y;
345  if (path->points[i+3].p.y < ymin) ymin = path->points[i+3].p.y;
346
347  double curv0 =
348    k0*(0.5-tpart[0])*(0.5-tpart[0]) +
349    k1*(0.5-tpart[0])*(0.5+tpart[0]) +
350    k2*(0.5+tpart[0])*(0.5+tpart[0]);
351
352  *nx = 0;
353
354  for (ipart=0; ipart<npart; ++ipart)
355    {
356      double curv1 =
357        k0*(0.5-tpart[ipart+1])*(0.5-tpart[ipart+1]) +
358        k1*(0.5-tpart[ipart+1])*(0.5+tpart[ipart+1]) +
359        k2*(0.5+tpart[ipart+1])*(0.5+tpart[ipart+1]);
360
361      int n_splines = 1;
362
363      if (fabs(xmax-xmin)<1.0e-8 && fabs(ymax-ymin) < 1.0e-8)
364        n_splines = (int)(16.0 * fabs(curv0-curv1)/((xmax-xmin)*(ymax-ymin)));
365
366      int j;
367
368      if (n_splines < 1)            n_splines = 1;
369      else if (n_splines > 8/npart) n_splines = 8/npart;
370
371      hpgs_point d1;
372      hpgs_point p1;
373
374      hpgs_bezier_path_point(path,i,tpart[ipart],&p1);
375      hpgs_bezier_path_tangent(path,i,tpart[ipart],1,1.0e-3,&d1);
376
377      for (j=1;j<=n_splines;++j)
378        {
379          double t = tpart[ipart] + (tpart[ipart+1]-tpart[ipart]) *  j / (double)n_splines;
380
381          hpgs_point d2;
382          hpgs_point p2;
383
384          hpgs_bezier_path_point(path,i,t,&p2);
385          hpgs_bezier_path_tangent(path,i,t,j<n_splines ? 0 : -1,1.0e-3,&d2);
386
387          add_quad(&p1,&d1,&p2,&d2,nx,points);
388
389          p1 = p2;
390          d1 = d2;
391       }
392
393      curv0 = curv1;
394    }
395}
396
397/* Internal: Calculates the derivation of the curve length
398     of the bezier curve \c b at curve parameter \c t.
399*/
400static double bezier_dl(const hpgs_paint_path *path, int i, double t)
401{
402  double tp,tm;
403
404  tp = t + 0.5;
405  tm = 0.5 - t;
406
407  return 3.0*hypot((path->points[i+1].p.x-path->points[i].p.x)*tm*tm+
408		   2.0 * (path->points[i+2].p.x-path->points[i+1].p.x) * tm*tp+
409		   (path->points[i+3].p.x-path->points[i+2].p.x) * tp*tp,
410
411		   (path->points[i+1].p.y-path->points[i].p.y)*tm*tm+
412		   2.0 * (path->points[i+2].p.y-path->points[i+1].p.y) * tm*tp+
413		   (path->points[i+3].p.y-path->points[i+2].p.y) * tp*tp
414		   );
415}
416
417/*! Calculates the curve lengths at an equidistant grid
418    of curve parmeters by na numerical quadrature.
419
420    You must call this subroutine before using
421    \c hpgs_bezier_length_param.
422*/
423void hpgs_bezier_length_init(hpgs_bezier_length *b,
424                             const hpgs_paint_path *path, int i)
425{
426  int ip;
427
428  b->dl[0] = bezier_dl(path,i,-0.5);
429
430  b->l[0] = 0.0;
431
432  for (ip=1;ip<=16;++ip)
433    {
434      double dlm0,dlm1;
435
436      // Lobatto quadrature with order 4.
437      b->dl[ip] = bezier_dl(path,i,(ip-8.0)*0.0625);
438
439      dlm0 = bezier_dl(path,i,(ip-(8.5+0.22360679774997896964))*0.0625);
440      dlm1 = bezier_dl(path,i,(ip-(8.5-0.22360679774997896964))*0.0625);
441
442      b->l[ip] = b->l[ip-1] +
443	(b->dl[ip-1] + 5.0 * (dlm0 + dlm1) + b->dl[ip])/(16.0*12.0);
444    }
445}
446
447/*! Returns the curve parameter at the bezier curve \c b
448    corresponding to a curve length \c l.
449
450    You must have called \c hpgs_bezier_length_init before using
451    this function.
452*/
453double hpgs_bezier_length_param(const hpgs_bezier_length *b, double l)
454{
455  // binary search.
456
457  int i = 0;
458
459  if (l >= b->l[i+8]) i+=8;
460  if (l >= b->l[i+4]) i+=4;
461  if (l >= b->l[i+2]) i+=2;
462  if (l >= b->l[i+1])  ++i;
463
464  double dl = b->l[i+1]-b->l[i];
465  double ll = (l-b->l[i])/dl;
466
467  return
468    (i - 8.0 + ll*ll*(3.0-2.0*ll)) * 0.0625 +
469    (ll*(ll*(ll-2.0)+1.0) / b->dl[i] + ll*ll*(ll-1.0) / b->dl[i+1] ) * dl;
470}
471
472/*! Calculates the minimal x component of all four control points of
473    the bezier spline. */
474double hpgs_bezier_path_xmin (const hpgs_paint_path *path, int i)
475{
476  double x0 =
477    path->points[i].p.x < path->points[i+1].p.x ?
478    path->points[i].p.x : path->points[i+1].p.x;
479
480  double x1 =
481    path->points[i+2].p.x < path->points[i+3].p.x ?
482    path->points[i+2].p.x : path->points[i+3].p.x;
483
484  return x0 < x1 ? x0 : x1;
485}
486
487/*! Calculates the maximal x component of all four control points of
488    the bezier spline. */
489double hpgs_bezier_path_xmax (const hpgs_paint_path *path, int i)
490{
491  double x0 =
492    path->points[i].p.x > path->points[i+1].p.x ?
493    path->points[i].p.x : path->points[i+1].p.x;
494
495  double x1 =
496    path->points[i+2].p.x > path->points[i+3].p.x ?
497    path->points[i+2].p.x : path->points[i+3].p.x;
498
499  return x0 > x1 ? x0 : x1;
500}
501
502/*! Calculates the minimal y component of all four control points of
503    the bezier spline. */
504double hpgs_bezier_path_ymin (const hpgs_paint_path *path, int i)
505{
506  double y0 =
507    path->points[i].p.y < path->points[i+1].p.y ?
508    path->points[i].p.y : path->points[i+1].p.y;
509
510  double y1 =
511    path->points[i+2].p.y < path->points[i+3].p.y ?
512    path->points[i+2].p.y : path->points[i+3].p.y;
513
514  return y0 < y1 ? y0 : y1;
515}
516
517/*! Calculates the maximal y component of all four control points of
518    the bezier spline. */
519double hpgs_bezier_path_ymax (const hpgs_paint_path *path, int i)
520{
521  double y0 =
522    path->points[i].p.y > path->points[i+1].p.y ?
523    path->points[i].p.y : path->points[i+1].p.y;
524
525  double y1 =
526    path->points[i+2].p.y > path->points[i+3].p.y ?
527    path->points[i+2].p.y : path->points[i+3].p.y;
528
529  return y0 > y1 ? y0 : y1;
530}
531
532/*! Calculates the x component of a point on the bezier spline. */
533double hpgs_bezier_path_x (const hpgs_paint_path *path, int i, double t)
534{
535  double tp,tm;
536
537  tp = t + 0.5;
538  tm = 0.5 - t;
539
540  return
541    path->points[i].p.x *tm*tm*tm +
542    3.0 * tm*tp * (path->points[i+1].p.x * tm + path->points[i+2].p.x * tp) +
543    path->points[i+3].p.x * tp*tp*tp;
544}
545
546/*! Calculates the y component of a point on the bezier spline. */
547double hpgs_bezier_path_y (const hpgs_paint_path *path, int i, double t)
548{
549  double tp,tm;
550
551  tp = t + 0.5;
552  tm = 0.5 - t;
553
554  return
555    path->points[i].p.y *tm*tm*tm +
556    3.0 * tm*tp * (path->points[i+1].p.y * tm + path->points[i+2].p.y * tp) +
557    path->points[i+3].p.y * tp*tp*tp;
558}
559
560/*! Calculates the point on the bezier curve starting
561    at point \c i of \c path at curve parameter \c t.
562*/
563void hpgs_bezier_path_point(const hpgs_paint_path *path, int i,
564			    double t, hpgs_point *p)
565{
566  double tp,tm;
567
568  tp = t + 0.5;
569  tm = 0.5 - t;
570
571  p->x =
572    path->points[i].p.x *tm*tm*tm +
573    3.0 * tm*tp * (path->points[i+1].p.x * tm + path->points[i+2].p.x * tp) +
574    path->points[i+3].p.x * tp*tp*tp;
575
576  p->y =
577    path->points[i].p.y *tm*tm*tm +
578    3.0 * tm*tp * (path->points[i+1].p.y * tm + path->points[i+2].p.y * tp) +
579    path->points[i+3].p.y * tp*tp*tp;
580}
581