1 /*
2 * Source code from Xfig 3.2.4 modified to work with arrays of doubles
3 * instead linked lists of F_points and to remove some globals(!)
4 * See copyright etc below.
5 *
6 * #included from engine.c
7 * That manages the R_alloc stack.
8 */
9
10 /*
11 * From w_drawprim.h
12 */
13
14 #define MAXNUMPTS 25000
15
16 /*
17 * From u_draw.c
18 */
19
20 /*
21 * FIG : Facility for Interactive Generation of figures
22 * Copyright (c) 1985-1988 by Supoj Sutanthavibul
23 * Parts Copyright (c) 1989-2002 by Brian V. Smith
24 * Parts Copyright (c) 1991 by Paul King
25 * Parts Copyright (c) 1992 by James Tough
26 * Parts Copyright (c) 1998 by Georg Stemmer
27 * Parts Copyright (c) 1995 by C. Blanc and C. Schlick
28 *
29 * Any party obtaining a copy of these files is granted, free of charge, a
30 * full and unrestricted irrevocable, world-wide, paid up, royalty-free,
31 * nonexclusive right and license to deal in this software and
32 * documentation files (the "Software"), including without limitation the
33 * rights to use, copy, modify, merge, publish and/or distribute copies of
34 * the Software, and to permit persons who receive copies from any such
35 * party to do so, with the only requirement being that this copyright
36 * notice remain intact.
37 *
38 */
39
40 /************** POLYGON/CURVE DRAWING FACILITIES ****************/
41
42 static int npoints;
43 static int max_points;
44 static double *xpoints;
45 static double *ypoints;
46
47 /************* Code begins here *************/
48
49 /* R_allocs or mallocs global arrays */
50 static Rboolean
add_point(double x,double y,pGEDevDesc dd)51 add_point(double x, double y, pGEDevDesc dd)
52 {
53 if (npoints >= max_points) {
54 int tmp_n;
55 double *tmp_px;
56 double *tmp_py;
57 tmp_n = max_points + 200;
58 /* too many points, return false */
59 if (tmp_n > MAXNUMPTS) {
60 error(_("add_point - reached MAXNUMPTS (%d)"),tmp_n);
61 }
62 if (max_points == 0) {
63 tmp_px = (double *) R_alloc(tmp_n, sizeof(double));
64 tmp_py = (double *) R_alloc(tmp_n, sizeof(double));
65 } else {
66 tmp_px = (double *) S_realloc((char *) xpoints,
67 tmp_n, max_points,
68 sizeof(double));
69 tmp_py = (double *) S_realloc((char *) ypoints,
70 tmp_n, max_points,
71 sizeof(double));
72 }
73 if (tmp_px == NULL || tmp_py == NULL) {
74 error(_("insufficient memory to allocate point array"));
75 }
76 xpoints = tmp_px;
77 ypoints = tmp_py;
78 max_points = tmp_n;
79 }
80 /* ignore identical points */
81 if (npoints > 0 && xpoints[npoints-1] == x && ypoints[npoints-1] == y)
82 return TRUE;
83 /*
84 * Convert back from 1200ppi to DEVICE coordinates
85 */
86 xpoints[npoints] = toDeviceX(x / 1200, GE_INCHES, dd);
87 ypoints[npoints] = toDeviceY(y / 1200, GE_INCHES, dd);
88 npoints = npoints + 1;
89 return TRUE;
90 }
91
92 /*
93 * From u_draw_spline.c
94 */
95
96 /*
97 * FIG : Facility for Interactive Generation of figures
98 * Copyright (c) 1985-1988 by Supoj Sutanthavibul
99 * Parts Copyright (c) 1989-2002 by Brian V. Smith
100 * Parts Copyright (c) 1991 by Paul King
101 *
102 * Any party obtaining a copy of these files is granted, free of charge, a
103 * full and unrestricted irrevocable, world-wide, paid up, royalty-free,
104 * nonexclusive right and license to deal in this software and
105 * documentation files (the "Software"), including without limitation the
106 * rights to use, copy, modify, merge, publish and/or distribute copies of
107 * the Software, and to permit persons who receive copies from any such
108 * party to do so, with the only requirement being that this copyright
109 * notice remain intact.
110 *
111 */
112
113 /* THIS FILE IS #included FROM u_draw.c and u_geom.c */
114
115 /********************* CURVES FOR SPLINES *****************************
116
117 The following spline drawing routines are from
118
119 "X-splines : A Spline Model Designed for the End User"
120
121 by Carole BLANC and Christophe SCHLICK,
122 in Proceedings of SIGGRAPH ' 95
123
124 ***********************************************************************/
125
126 #define HIGH_PRECISION 0.5
127 #define LOW_PRECISION 1.0
128 #define ZOOM_PRECISION 5.0
129 #define ARROW_START 4
130 #define MAX_SPLINE_STEP 0.2
131
132 /***********************************************************************/
133
134 #define Q(s) (-(s))
135 #define EQN_NUMERATOR(dim) \
136 (A_blend[0]*dim[0]+A_blend[1]*dim[1]+A_blend[2]*dim[2]+A_blend[3]*dim[3])
137
138 static double
f_blend(double numerator,double denominator)139 f_blend(double numerator, double denominator)
140 {
141 double p = 2 * denominator * denominator;
142 double u = numerator / denominator;
143 double u2 = u * u;
144
145 return (u * u2 * (10 - p + (2*p - 15)*u + (6 - p)*u2));
146 }
147
148 static double
g_blend(double u,double q)149 g_blend(double u, double q) /* p equals 2 */
150 {
151 return(u*(q + u*(2*q + u*(8 - 12*q + u*(14*q - 11 + u*(4 - 5*q))))));
152 }
153
154 static double
h_blend(double u,double q)155 h_blend(double u, double q)
156 {
157 double u2=u*u;
158 return (u * (q + u * (2 * q + u2 * (-2*q - u*q))));
159 }
160
161 static void
negative_s1_influence(double t,double s1,double * A0,double * A2)162 negative_s1_influence(double t, double s1, double *A0, double *A2)
163 {
164 *A0 = h_blend(-t, Q(s1));
165 *A2 = g_blend(t, Q(s1));
166 }
167
168 static void
negative_s2_influence(double t,double s2,double * A1,double * A3)169 negative_s2_influence(double t, double s2, double *A1, double *A3)
170 {
171 *A1 = g_blend(1-t, Q(s2));
172 *A3 = h_blend(t-1, Q(s2));
173 }
174
175 static void
positive_s1_influence(double k,double t,double s1,double * A0,double * A2)176 positive_s1_influence(double k, double t, double s1, double *A0, double *A2)
177 {
178 double Tk;
179
180 Tk = k+1+s1;
181 *A0 = (t+k+1<Tk) ? f_blend(t+k+1-Tk, k-Tk) : 0.0;
182
183 Tk = k+1-s1;
184 *A2 = f_blend(t+k+1-Tk, k+2-Tk);
185 }
186
187 static void
positive_s2_influence(double k,double t,double s2,double * A1,double * A3)188 positive_s2_influence(double k, double t, double s2, double *A1, double *A3)
189 {
190 double Tk;
191
192 Tk = k+2+s2;
193 *A1 = f_blend(t+k+1-Tk, k+1-Tk);
194
195 Tk = k+2-s2;
196 *A3 = (t+k+1>Tk) ? f_blend(t+k+1-Tk, k+3-Tk) : 0.0;
197 }
198
199 static void
point_adding(double * A_blend,double * px,double * py,pGEDevDesc dd)200 point_adding(double *A_blend, double *px, double *py,
201 pGEDevDesc dd)
202 {
203 double weights_sum;
204
205 weights_sum = A_blend[0] + A_blend[1] + A_blend[2] + A_blend[3];
206 add_point(EQN_NUMERATOR(px) / (weights_sum),
207 EQN_NUMERATOR(py) / (weights_sum),
208 dd);
209 }
210
211 static void
point_computing(double * A_blend,double * px,double * py,double * x,double * y)212 point_computing(double *A_blend,
213 double *px, double *py,
214 double *x, double *y)
215 {
216 double weights_sum;
217
218 weights_sum = A_blend[0] + A_blend[1] + A_blend[2] + A_blend[3];
219
220 *x = EQN_NUMERATOR(px) / (weights_sum);
221 *y = EQN_NUMERATOR(py) / (weights_sum);
222 }
223
224 static double
step_computing(int k,double * px,double * py,double s1,double s2,double precision,pGEDevDesc dd)225 step_computing(int k,
226 double *px, double *py,
227 double s1, double s2,
228 double precision,
229 pGEDevDesc dd)
230 {
231 double A_blend[4];
232 double xstart, ystart, xend, yend, xmid, ymid, xlength, ylength;
233 double start_to_end_dist, devWidth, devHeight, devDiag, number_of_steps;
234 double step, angle_cos, scal_prod, xv1, xv2, yv1, yv2, sides_length_prod;
235
236 /* This function computes the step used to draw the segment (p1, p2)
237 (xv1, yv1) : coordinates of the vector from middle to origin
238 (xv2, yv2) : coordinates of the vector from middle to extremity */
239
240 if ((s1 == 0) && (s2 == 0))
241 return(1.0); /* only one step in case of linear segment */
242
243 /* compute coordinates of the origin */
244 if (s1>0) {
245 if (s2<0) {
246 positive_s1_influence(k, 0.0, s1, &A_blend[0], &A_blend[2]);
247 negative_s2_influence(0.0, s2, &A_blend[1], &A_blend[3]);
248 } else {
249 positive_s1_influence(k, 0.0, s1, &A_blend[0], &A_blend[2]);
250 positive_s2_influence(k, 0.0, s2, &A_blend[1], &A_blend[3]);
251 }
252 point_computing(A_blend, px, py, &xstart, &ystart);
253 } else {
254 xstart = px[1];
255 ystart = py[1];
256 }
257
258 /* compute coordinates of the extremity */
259 if (s2>0) {
260 if (s1<0) {
261 negative_s1_influence(1.0, s1, &A_blend[0], &A_blend[2]);
262 positive_s2_influence(k, 1.0, s2, &A_blend[1], &A_blend[3]);
263 } else {
264 positive_s1_influence(k, 1.0, s1, &A_blend[0], &A_blend[2]);
265 positive_s2_influence(k, 1.0, s2, &A_blend[1], &A_blend[3]);
266 }
267 point_computing(A_blend, px, py, &xend, ¥d);
268 } else {
269 xend = px[2];
270 yend = py[2];
271 }
272
273 /* compute coordinates of the middle */
274 if (s2>0) {
275 if (s1<0) {
276 negative_s1_influence(0.5, s1, &A_blend[0], &A_blend[2]);
277 positive_s2_influence(k, 0.5, s2, &A_blend[1], &A_blend[3]);
278 } else {
279 positive_s1_influence(k, 0.5, s1, &A_blend[0], &A_blend[2]);
280 positive_s2_influence(k, 0.5, s2, &A_blend[1], &A_blend[3]);
281 }
282 } else if (s1<0) {
283 negative_s1_influence(0.5, s1, &A_blend[0], &A_blend[2]);
284 negative_s2_influence(0.5, s2, &A_blend[1], &A_blend[3]);
285 } else {
286 positive_s1_influence(k, 0.5, s1, &A_blend[0], &A_blend[2]);
287 negative_s2_influence(0.5, s2, &A_blend[1], &A_blend[3]);
288 }
289 point_computing(A_blend, px, py, &xmid, &ymid);
290
291 xv1 = xstart - xmid;
292 yv1 = ystart - ymid;
293 xv2 = xend - xmid;
294 yv2 = yend - ymid;
295
296 scal_prod = xv1*xv2 + yv1*yv2;
297
298 sides_length_prod = sqrt((xv1*xv1 + yv1*yv1)*(xv2*xv2 + yv2*yv2));
299
300 /* compute cosinus of origin-middle-extremity angle, which approximates the
301 curve of the spline segment */
302 if (sides_length_prod == 0.0)
303 angle_cos = 0.0;
304 else
305 angle_cos = scal_prod/sides_length_prod;
306
307 xlength = xend - xstart;
308 ylength = yend - ystart;
309
310 start_to_end_dist = sqrt(xlength*xlength + ylength*ylength);
311
312 /* Paul 2009-01-25
313 * It is possible for origin and extremity to be very remote
314 * indeed (if the control points are located WAY off the device).
315 * In order to avoid having ridiculously many steps, limit
316 * the start_to_end_dist to being the length of the diagonal of the
317 * device.
318 */
319 devWidth = fromDeviceWidth(toDeviceWidth(1, GE_NDC, dd),
320 GE_INCHES, dd)*1200;
321 devHeight = fromDeviceHeight(toDeviceHeight(1, GE_NDC, dd),
322 GE_INCHES, dd)*1200;
323 devDiag = sqrt(devWidth* devWidth + devHeight*devHeight);
324 if (start_to_end_dist > devDiag)
325 start_to_end_dist = devDiag;
326
327 /* more steps if segment's origin and extremity are remote */
328 number_of_steps = sqrt(start_to_end_dist)/2;
329
330 /* more steps if the curve is high */
331 number_of_steps += (int)((1 + angle_cos)*10);
332
333 if (number_of_steps == 0)
334 step = 1;
335 else
336 step = precision/number_of_steps;
337
338 if ((step > MAX_SPLINE_STEP) || (step == 0))
339 step = MAX_SPLINE_STEP;
340
341 return (step);
342 }
343
344 static void
spline_segment_computing(double step,int k,double * px,double * py,double s1,double s2,pGEDevDesc dd)345 spline_segment_computing(double step, int k,
346 double *px, double *py,
347 double s1, double s2,
348 pGEDevDesc dd)
349 {
350 double A_blend[4];
351 double t;
352
353 if (s1<0) {
354 if (s2<0) {
355 for (t=0.0 ; t<1 ; t+=step) {
356 negative_s1_influence(t, s1, &A_blend[0], &A_blend[2]);
357 negative_s2_influence(t, s2, &A_blend[1], &A_blend[3]);
358
359 point_adding(A_blend, px, py, dd);
360 }
361 } else {
362 for (t=0.0 ; t<1 ; t+=step) {
363 negative_s1_influence(t, s1, &A_blend[0], &A_blend[2]);
364 positive_s2_influence(k, t, s2, &A_blend[1], &A_blend[3]);
365
366 point_adding(A_blend, px, py, dd);
367 }
368 }
369 } else if (s2<0) {
370 for (t=0.0 ; t<1 ; t+=step) {
371 positive_s1_influence(k, t, s1, &A_blend[0], &A_blend[2]);
372 negative_s2_influence(t, s2, &A_blend[1], &A_blend[3]);
373
374 point_adding(A_blend, px, py, dd);
375 }
376 } else {
377 for (t=0.0 ; t<1 ; t+=step) {
378 positive_s1_influence(k, t, s1, &A_blend[0], &A_blend[2]);
379 positive_s2_influence(k, t, s2, &A_blend[1], &A_blend[3]);
380
381 point_adding(A_blend, px, py, dd);
382 }
383 }
384 }
385
386 /*
387 * For adding last line segment when computing open spline
388 * WITHOUT end control points repeated
389 * (i.e., can't just connect to last control point)
390 */
391 static void
spline_last_segment_computing(double step,int k,double * px,double * py,double s1,double s2,pGEDevDesc dd)392 spline_last_segment_computing(double step, int k,
393 double *px, double *py,
394 double s1, double s2,
395 pGEDevDesc dd)
396 {
397 double A_blend[4];
398 double t = 1;
399
400 if (s1<0) {
401 if (s2<0) {
402 negative_s1_influence(t, s1, &A_blend[0], &A_blend[2]);
403 negative_s2_influence(t, s2, &A_blend[1], &A_blend[3]);
404
405 point_adding(A_blend, px, py, dd);
406 } else {
407 negative_s1_influence(t, s1, &A_blend[0], &A_blend[2]);
408 positive_s2_influence(k, t, s2, &A_blend[1], &A_blend[3]);
409
410 point_adding(A_blend, px, py, dd);
411 }
412 } else if (s2<0) {
413 positive_s1_influence(k, t, s1, &A_blend[0], &A_blend[2]);
414 negative_s2_influence(t, s2, &A_blend[1], &A_blend[3]);
415
416 point_adding(A_blend, px, py, dd);
417 } else {
418 positive_s1_influence(k, t, s1, &A_blend[0], &A_blend[2]);
419 positive_s2_influence(k, t, s2, &A_blend[1], &A_blend[3]);
420
421 point_adding(A_blend, px, py, dd);
422 }
423 }
424
425 /********************* MAIN METHODS *************************************/
426
427 /*
428 * x and y are in DEVICE coordinates
429 * xfig works in 1200ppi
430 * (http://www.csit.fsu.edu/~burkardt/data/fig/fig_format.html)
431 * so convert to 1200ppi so that step calculations are correct
432 */
433 #define COPY_CONTROL_POINT(PI, I, N) \
434 px[PI] = fromDeviceX(x[(I) % N], GE_INCHES, dd) * 1200; \
435 py[PI] = fromDeviceY(y[(I) % N], GE_INCHES, dd) * 1200; \
436 ps[PI] = s[(I) % N]
437
438 #define NEXT_CONTROL_POINTS(K, N) \
439 COPY_CONTROL_POINT(0, K, N); \
440 COPY_CONTROL_POINT(1, K + 1, N); \
441 COPY_CONTROL_POINT(2, K + 2, N); \
442 COPY_CONTROL_POINT(3, K + 3, N)
443
444 #define INIT_CONTROL_POINTS(N) \
445 COPY_CONTROL_POINT(0, N - 1, N); \
446 COPY_CONTROL_POINT(1, 0, N); \
447 COPY_CONTROL_POINT(2, 1, N); \
448 COPY_CONTROL_POINT(3, 2, N)
449
450 #define SPLINE_SEGMENT_LOOP(K, PX, PY, S1, S2, PREC) \
451 step = step_computing(K, PX, PY, S1, S2, PREC, dd); \
452 spline_segment_computing(step, K, PX, PY, S1, S2, dd)
453
454 static Rboolean
compute_open_spline(int n,double * x,double * y,double * s,Rboolean repEnds,double precision,pGEDevDesc dd)455 compute_open_spline(int n, double *x, double *y, double *s,
456 Rboolean repEnds,
457 double precision,
458 pGEDevDesc dd)
459 {
460 int k;
461 double step = 0.0 /* -Wall */;
462 double px[4];
463 double py[4];
464 double ps[4]={0.,0.,0.,0.};
465
466 max_points = 0;
467 npoints = 0;
468 xpoints = NULL;
469 ypoints = NULL;
470
471 if (repEnds && n < 2)
472 error(_("there must be at least two control points"));
473 if (!repEnds && n < 4)
474 error(_("there must be at least four control points"));
475
476 if (repEnds) {
477 /* first control point is needed twice for the first segment */
478 COPY_CONTROL_POINT(0, 0, n);
479 COPY_CONTROL_POINT(1, 0, n);
480 COPY_CONTROL_POINT(2, 1, n);
481
482 if (n == 2) {
483 COPY_CONTROL_POINT(3, 1, n);
484 } else {
485 COPY_CONTROL_POINT(3, 2, n);
486 }
487
488 for (k = 0 ; ; k++) {
489 SPLINE_SEGMENT_LOOP(k, px, py, ps[1], ps[2], precision);
490 /* Paul 2008-12-05
491 * >= rather than == to handle special case of n == 2
492 */
493 if (k + 3 >= n)
494 break;
495 NEXT_CONTROL_POINTS(k, n);
496 }
497
498 /* last control point is needed twice for the last segment */
499 COPY_CONTROL_POINT(0, n - 3, n);
500 COPY_CONTROL_POINT(1, n - 2, n);
501 COPY_CONTROL_POINT(2, n - 1, n);
502 COPY_CONTROL_POINT(3, n - 1, n);
503 SPLINE_SEGMENT_LOOP(k, px, py, ps[1], ps[2], precision);
504
505 add_point(px[3], py[3], dd);
506 } else {
507 for (k = 0 ; k + 3 < n ; k++) {
508 NEXT_CONTROL_POINTS(k, n);
509 SPLINE_SEGMENT_LOOP(k, px, py, ps[1], ps[2], precision);
510 }
511 spline_last_segment_computing(step, n - 4, px, py, ps[1], ps[2], dd);
512 }
513
514 return TRUE;
515 }
516
517 static Rboolean
compute_closed_spline(int n,double * x,double * y,double * s,double precision,pGEDevDesc dd)518 compute_closed_spline(int n, double *x, double *y, double *s,
519 double precision,
520 pGEDevDesc dd)
521 {
522 int k;
523 double step;
524 double px[4];
525 double py[4];
526 double ps[4];
527
528 max_points = 0;
529 npoints = 0;
530 xpoints = NULL;
531 ypoints = NULL;
532
533 if (n < 3)
534 error(_("There must be at least three control points"));
535
536 INIT_CONTROL_POINTS(n);
537
538 for (k = 0 ; k < n ; k++) {
539 SPLINE_SEGMENT_LOOP(k, px, py, ps[1], ps[2], precision);
540 NEXT_CONTROL_POINTS(k, n);
541 }
542
543 return TRUE;
544 }
545