1 /* GIMP - The GNU Image Manipulation Program
2 * Copyright (C) 1995-2001 Spencer Kimball, Peter Mattis, and others
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17
18 #include "config.h"
19
20 #include <string.h>
21
22 #include <glib-object.h>
23
24 #include "libgimpmath/gimpmath.h"
25
26 #include "core-types.h"
27
28 #include "gimp-transform-utils.h"
29 #include "gimpcoords.h"
30 #include "gimpcoords-interpolate.h"
31
32
33 #define EPSILON 1e-6
34
35
36 void
gimp_transform_get_rotate_center(gint x,gint y,gint width,gint height,gboolean auto_center,gdouble * center_x,gdouble * center_y)37 gimp_transform_get_rotate_center (gint x,
38 gint y,
39 gint width,
40 gint height,
41 gboolean auto_center,
42 gdouble *center_x,
43 gdouble *center_y)
44 {
45 g_return_if_fail (center_x != NULL);
46 g_return_if_fail (center_y != NULL);
47
48 if (auto_center)
49 {
50 *center_x = (gdouble) x + (gdouble) width / 2.0;
51 *center_y = (gdouble) y + (gdouble) height / 2.0;
52 }
53 }
54
55 void
gimp_transform_get_flip_axis(gint x,gint y,gint width,gint height,GimpOrientationType flip_type,gboolean auto_center,gdouble * axis)56 gimp_transform_get_flip_axis (gint x,
57 gint y,
58 gint width,
59 gint height,
60 GimpOrientationType flip_type,
61 gboolean auto_center,
62 gdouble *axis)
63 {
64 g_return_if_fail (axis != NULL);
65
66 if (auto_center)
67 {
68 switch (flip_type)
69 {
70 case GIMP_ORIENTATION_HORIZONTAL:
71 *axis = ((gdouble) x + (gdouble) width / 2.0);
72 break;
73
74 case GIMP_ORIENTATION_VERTICAL:
75 *axis = ((gdouble) y + (gdouble) height / 2.0);
76 break;
77
78 default:
79 g_return_if_reached ();
80 break;
81 }
82 }
83 }
84
85 void
gimp_transform_matrix_flip(GimpMatrix3 * matrix,GimpOrientationType flip_type,gdouble axis)86 gimp_transform_matrix_flip (GimpMatrix3 *matrix,
87 GimpOrientationType flip_type,
88 gdouble axis)
89 {
90 g_return_if_fail (matrix != NULL);
91
92 switch (flip_type)
93 {
94 case GIMP_ORIENTATION_HORIZONTAL:
95 gimp_matrix3_translate (matrix, - axis, 0.0);
96 gimp_matrix3_scale (matrix, -1.0, 1.0);
97 gimp_matrix3_translate (matrix, axis, 0.0);
98 break;
99
100 case GIMP_ORIENTATION_VERTICAL:
101 gimp_matrix3_translate (matrix, 0.0, - axis);
102 gimp_matrix3_scale (matrix, 1.0, -1.0);
103 gimp_matrix3_translate (matrix, 0.0, axis);
104 break;
105
106 case GIMP_ORIENTATION_UNKNOWN:
107 break;
108 }
109 }
110
111 void
gimp_transform_matrix_flip_free(GimpMatrix3 * matrix,gdouble x1,gdouble y1,gdouble x2,gdouble y2)112 gimp_transform_matrix_flip_free (GimpMatrix3 *matrix,
113 gdouble x1,
114 gdouble y1,
115 gdouble x2,
116 gdouble y2)
117 {
118 gdouble angle;
119
120 g_return_if_fail (matrix != NULL);
121
122 angle = atan2 (y2 - y1, x2 - x1);
123
124 gimp_matrix3_identity (matrix);
125 gimp_matrix3_translate (matrix, -x1, -y1);
126 gimp_matrix3_rotate (matrix, -angle);
127 gimp_matrix3_scale (matrix, 1.0, -1.0);
128 gimp_matrix3_rotate (matrix, angle);
129 gimp_matrix3_translate (matrix, x1, y1);
130 }
131
132 void
gimp_transform_matrix_rotate(GimpMatrix3 * matrix,GimpRotationType rotate_type,gdouble center_x,gdouble center_y)133 gimp_transform_matrix_rotate (GimpMatrix3 *matrix,
134 GimpRotationType rotate_type,
135 gdouble center_x,
136 gdouble center_y)
137 {
138 gdouble angle = 0;
139
140 switch (rotate_type)
141 {
142 case GIMP_ROTATE_90:
143 angle = G_PI_2;
144 break;
145 case GIMP_ROTATE_180:
146 angle = G_PI;
147 break;
148 case GIMP_ROTATE_270:
149 angle = - G_PI_2;
150 break;
151 }
152
153 gimp_transform_matrix_rotate_center (matrix, center_x, center_y, angle);
154 }
155
156 void
gimp_transform_matrix_rotate_rect(GimpMatrix3 * matrix,gint x,gint y,gint width,gint height,gdouble angle)157 gimp_transform_matrix_rotate_rect (GimpMatrix3 *matrix,
158 gint x,
159 gint y,
160 gint width,
161 gint height,
162 gdouble angle)
163 {
164 gdouble center_x;
165 gdouble center_y;
166
167 g_return_if_fail (matrix != NULL);
168
169 center_x = (gdouble) x + (gdouble) width / 2.0;
170 center_y = (gdouble) y + (gdouble) height / 2.0;
171
172 gimp_matrix3_translate (matrix, -center_x, -center_y);
173 gimp_matrix3_rotate (matrix, angle);
174 gimp_matrix3_translate (matrix, +center_x, +center_y);
175 }
176
177 void
gimp_transform_matrix_rotate_center(GimpMatrix3 * matrix,gdouble center_x,gdouble center_y,gdouble angle)178 gimp_transform_matrix_rotate_center (GimpMatrix3 *matrix,
179 gdouble center_x,
180 gdouble center_y,
181 gdouble angle)
182 {
183 g_return_if_fail (matrix != NULL);
184
185 gimp_matrix3_translate (matrix, -center_x, -center_y);
186 gimp_matrix3_rotate (matrix, angle);
187 gimp_matrix3_translate (matrix, +center_x, +center_y);
188 }
189
190 void
gimp_transform_matrix_scale(GimpMatrix3 * matrix,gint x,gint y,gint width,gint height,gdouble t_x,gdouble t_y,gdouble t_width,gdouble t_height)191 gimp_transform_matrix_scale (GimpMatrix3 *matrix,
192 gint x,
193 gint y,
194 gint width,
195 gint height,
196 gdouble t_x,
197 gdouble t_y,
198 gdouble t_width,
199 gdouble t_height)
200 {
201 gdouble scale_x = 1.0;
202 gdouble scale_y = 1.0;
203
204 g_return_if_fail (matrix != NULL);
205
206 if (width > 0)
207 scale_x = t_width / (gdouble) width;
208
209 if (height > 0)
210 scale_y = t_height / (gdouble) height;
211
212 gimp_matrix3_identity (matrix);
213 gimp_matrix3_translate (matrix, -x, -y);
214 gimp_matrix3_scale (matrix, scale_x, scale_y);
215 gimp_matrix3_translate (matrix, t_x, t_y);
216 }
217
218 void
gimp_transform_matrix_shear(GimpMatrix3 * matrix,gint x,gint y,gint width,gint height,GimpOrientationType orientation,gdouble amount)219 gimp_transform_matrix_shear (GimpMatrix3 *matrix,
220 gint x,
221 gint y,
222 gint width,
223 gint height,
224 GimpOrientationType orientation,
225 gdouble amount)
226 {
227 gdouble center_x;
228 gdouble center_y;
229
230 g_return_if_fail (matrix != NULL);
231
232 if (width == 0)
233 width = 1;
234
235 if (height == 0)
236 height = 1;
237
238 center_x = (gdouble) x + (gdouble) width / 2.0;
239 center_y = (gdouble) y + (gdouble) height / 2.0;
240
241 gimp_matrix3_identity (matrix);
242 gimp_matrix3_translate (matrix, -center_x, -center_y);
243
244 if (orientation == GIMP_ORIENTATION_HORIZONTAL)
245 gimp_matrix3_xshear (matrix, amount / height);
246 else
247 gimp_matrix3_yshear (matrix, amount / width);
248
249 gimp_matrix3_translate (matrix, +center_x, +center_y);
250 }
251
252 void
gimp_transform_matrix_perspective(GimpMatrix3 * matrix,gint x,gint y,gint width,gint height,gdouble t_x1,gdouble t_y1,gdouble t_x2,gdouble t_y2,gdouble t_x3,gdouble t_y3,gdouble t_x4,gdouble t_y4)253 gimp_transform_matrix_perspective (GimpMatrix3 *matrix,
254 gint x,
255 gint y,
256 gint width,
257 gint height,
258 gdouble t_x1,
259 gdouble t_y1,
260 gdouble t_x2,
261 gdouble t_y2,
262 gdouble t_x3,
263 gdouble t_y3,
264 gdouble t_x4,
265 gdouble t_y4)
266 {
267 GimpMatrix3 trafo;
268 gdouble scalex;
269 gdouble scaley;
270
271 g_return_if_fail (matrix != NULL);
272
273 scalex = scaley = 1.0;
274
275 if (width > 0)
276 scalex = 1.0 / (gdouble) width;
277
278 if (height > 0)
279 scaley = 1.0 / (gdouble) height;
280
281 gimp_matrix3_translate (matrix, -x, -y);
282 gimp_matrix3_scale (matrix, scalex, scaley);
283
284 /* Determine the perspective transform that maps from
285 * the unit cube to the transformed coordinates
286 */
287 {
288 gdouble dx1, dx2, dx3, dy1, dy2, dy3;
289
290 dx1 = t_x2 - t_x4;
291 dx2 = t_x3 - t_x4;
292 dx3 = t_x1 - t_x2 + t_x4 - t_x3;
293
294 dy1 = t_y2 - t_y4;
295 dy2 = t_y3 - t_y4;
296 dy3 = t_y1 - t_y2 + t_y4 - t_y3;
297
298 /* Is the mapping affine? */
299 if ((dx3 == 0.0) && (dy3 == 0.0))
300 {
301 trafo.coeff[0][0] = t_x2 - t_x1;
302 trafo.coeff[0][1] = t_x4 - t_x2;
303 trafo.coeff[0][2] = t_x1;
304 trafo.coeff[1][0] = t_y2 - t_y1;
305 trafo.coeff[1][1] = t_y4 - t_y2;
306 trafo.coeff[1][2] = t_y1;
307 trafo.coeff[2][0] = 0.0;
308 trafo.coeff[2][1] = 0.0;
309 }
310 else
311 {
312 gdouble det1, det2;
313
314 det1 = dx3 * dy2 - dy3 * dx2;
315 det2 = dx1 * dy2 - dy1 * dx2;
316
317 trafo.coeff[2][0] = (det2 == 0.0) ? 1.0 : det1 / det2;
318
319 det1 = dx1 * dy3 - dy1 * dx3;
320
321 trafo.coeff[2][1] = (det2 == 0.0) ? 1.0 : det1 / det2;
322
323 trafo.coeff[0][0] = t_x2 - t_x1 + trafo.coeff[2][0] * t_x2;
324 trafo.coeff[0][1] = t_x3 - t_x1 + trafo.coeff[2][1] * t_x3;
325 trafo.coeff[0][2] = t_x1;
326
327 trafo.coeff[1][0] = t_y2 - t_y1 + trafo.coeff[2][0] * t_y2;
328 trafo.coeff[1][1] = t_y3 - t_y1 + trafo.coeff[2][1] * t_y3;
329 trafo.coeff[1][2] = t_y1;
330 }
331
332 trafo.coeff[2][2] = 1.0;
333 }
334
335 gimp_matrix3_mult (&trafo, matrix);
336 }
337
338 /* modified gaussian algorithm
339 * solves a system of linear equations
340 *
341 * Example:
342 * 1x + 2y + 4z = 25
343 * 2x + 1y = 4
344 * 3x + 5y + 2z = 23
345 * Solution: x=1, y=2, z=5
346 *
347 * Input:
348 * matrix = { 1,2,4,25,2,1,0,4,3,5,2,23 }
349 * s = 3 (Number of variables)
350 * Output:
351 * return value == TRUE (TRUE, if there is a single unique solution)
352 * solution == { 1,2,5 } (if the return value is FALSE, the content
353 * of solution is of no use)
354 */
355 static gboolean
mod_gauss(gdouble matrix[],gdouble solution[],gint s)356 mod_gauss (gdouble matrix[],
357 gdouble solution[],
358 gint s)
359 {
360 gint p[s]; /* row permutation */
361 gint i, j, r, temp;
362 gdouble q;
363 gint t = s + 1;
364
365 for (i = 0; i < s; i++)
366 {
367 p[i] = i;
368 }
369
370 for (r = 0; r < s; r++)
371 {
372 /* make sure that (r,r) is not 0 */
373 if (fabs (matrix[p[r] * t + r]) <= EPSILON)
374 {
375 /* we need to permutate rows */
376 for (i = r + 1; i <= s; i++)
377 {
378 if (i == s)
379 {
380 /* if this happens, the linear system has zero or
381 * more than one solutions.
382 */
383 return FALSE;
384 }
385
386 if (fabs (matrix[p[i] * t + r]) > EPSILON)
387 break;
388 }
389
390 temp = p[r];
391 p[r] = p[i];
392 p[i] = temp;
393 }
394
395 /* make (r,r) == 1 */
396 q = 1.0 / matrix[p[r] * t + r];
397 matrix[p[r] * t + r] = 1.0;
398
399 for (j = r + 1; j < t; j++)
400 {
401 matrix[p[r] * t + j] *= q;
402 }
403
404 /* make that all entries in column r are 0 (except (r,r)) */
405 for (i = 0; i < s; i++)
406 {
407 if (i == r)
408 continue;
409
410 for (j = r + 1; j < t ; j++)
411 {
412 matrix[p[i] * t + j] -= matrix[p[r] * t + j] * matrix[p[i] * t + r];
413 }
414
415 /* we don't need to execute the following line
416 * since we won't access this element again:
417 *
418 * matrix[p[i] * t + r] = 0.0;
419 */
420 }
421 }
422
423 for (i = 0; i < s; i++)
424 {
425 solution[i] = matrix[p[i] * t + s];
426 }
427
428 return TRUE;
429 }
430
431 /* multiplies 'matrix' by the matrix that transforms a set of 4 'input_points'
432 * to corresponding 'output_points', if such matrix exists, and is valid (i.e.,
433 * keeps the output points in front of the camera).
434 *
435 * returns TRUE if successful.
436 */
437 gboolean
gimp_transform_matrix_generic(GimpMatrix3 * matrix,const GimpVector2 input_points[4],const GimpVector2 output_points[4])438 gimp_transform_matrix_generic (GimpMatrix3 *matrix,
439 const GimpVector2 input_points[4],
440 const GimpVector2 output_points[4])
441 {
442 GimpMatrix3 trafo;
443 gdouble coeff[8 * 9];
444 gboolean negative = -1;
445 gint i;
446 gboolean result = TRUE;
447
448 g_return_val_if_fail (matrix != NULL, FALSE);
449 g_return_val_if_fail (input_points != NULL, FALSE);
450 g_return_val_if_fail (output_points != NULL, FALSE);
451
452 /* find the matrix that transforms 'input_points' to 'output_points', whose
453 * (3, 3) coeffcient is 1, by solving a system of linear equations whose
454 * solution is the remaining 8 coefficients.
455 */
456 for (i = 0; i < 4; i++)
457 {
458 coeff[i * 9 + 0] = input_points[i].x;
459 coeff[i * 9 + 1] = input_points[i].y;
460 coeff[i * 9 + 2] = 1.0;
461 coeff[i * 9 + 3] = 0.0;
462 coeff[i * 9 + 4] = 0.0;
463 coeff[i * 9 + 5] = 0.0;
464 coeff[i * 9 + 6] = -input_points[i].x * output_points[i].x;
465 coeff[i * 9 + 7] = -input_points[i].y * output_points[i].x;
466 coeff[i * 9 + 8] = output_points[i].x;
467
468 coeff[(i + 4) * 9 + 0] = 0.0;
469 coeff[(i + 4) * 9 + 1] = 0.0;
470 coeff[(i + 4) * 9 + 2] = 0.0;
471 coeff[(i + 4) * 9 + 3] = input_points[i].x;
472 coeff[(i + 4) * 9 + 4] = input_points[i].y;
473 coeff[(i + 4) * 9 + 5] = 1.0;
474 coeff[(i + 4) * 9 + 6] = -input_points[i].x * output_points[i].y;
475 coeff[(i + 4) * 9 + 7] = -input_points[i].y * output_points[i].y;
476 coeff[(i + 4) * 9 + 8] = output_points[i].y;
477 }
478
479 /* if there is no solution, bail */
480 if (! mod_gauss (coeff, (gdouble *) trafo.coeff, 8))
481 return FALSE;
482
483 trafo.coeff[2][2] = 1.0;
484
485 /* make sure that none of the input points maps to a point at infinity, and
486 * that all output points are on the same side of the camera.
487 */
488 for (i = 0; i < 4; i++)
489 {
490 gdouble w;
491 gboolean neg;
492
493 w = trafo.coeff[2][0] * input_points[i].x +
494 trafo.coeff[2][1] * input_points[i].y +
495 trafo.coeff[2][2];
496
497 if (fabs (w) <= EPSILON)
498 result = FALSE;
499
500 neg = (w < 0.0);
501
502 if (negative < 0)
503 {
504 negative = neg;
505 }
506 else if (neg != negative)
507 {
508 result = FALSE;
509 break;
510 }
511 }
512
513 /* if the output points are all behind the camera, negate the matrix, which
514 * would map the input points to the corresponding points in front of the
515 * camera.
516 */
517 if (negative > 0)
518 {
519 gint r;
520 gint c;
521
522 for (r = 0; r < 3; r++)
523 {
524 for (c = 0; c < 3; c++)
525 {
526 trafo.coeff[r][c] = -trafo.coeff[r][c];
527 }
528 }
529 }
530
531 /* append the transformation to 'matrix' */
532 gimp_matrix3_mult (&trafo, matrix);
533
534 return result;
535 }
536
537 gboolean
gimp_transform_polygon_is_convex(gdouble x1,gdouble y1,gdouble x2,gdouble y2,gdouble x3,gdouble y3,gdouble x4,gdouble y4)538 gimp_transform_polygon_is_convex (gdouble x1,
539 gdouble y1,
540 gdouble x2,
541 gdouble y2,
542 gdouble x3,
543 gdouble y3,
544 gdouble x4,
545 gdouble y4)
546 {
547 gdouble z1, z2, z3, z4;
548
549 /* We test if the transformed polygon is convex. if z1 and z2 have
550 * the same sign as well as z3 and z4 the polygon is convex.
551 */
552 z1 = ((x2 - x1) * (y4 - y1) -
553 (x4 - x1) * (y2 - y1));
554 z2 = ((x4 - x1) * (y3 - y1) -
555 (x3 - x1) * (y4 - y1));
556 z3 = ((x4 - x2) * (y3 - y2) -
557 (x3 - x2) * (y4 - y2));
558 z4 = ((x3 - x2) * (y1 - y2) -
559 (x1 - x2) * (y3 - y2));
560
561 return (z1 * z2 > 0) && (z3 * z4 > 0);
562 }
563
564 /* transforms the polygon or polyline, whose vertices are given by 'vertices',
565 * by 'matrix', performing clipping by the near plane. 'closed' indicates
566 * whether the vertices represent a polygon ('closed == TRUE') or a polyline
567 * ('closed == FALSE').
568 *
569 * returns the transformed vertices in 't_vertices', and their count in
570 * 'n_t_vertices'. the minimal possible number of transformed vertices is 0,
571 * which happens when the entire input is clipped. in general, the maximal
572 * possible number of transformed vertices is '3 * n_vertices / 2' (rounded
573 * down), however, for convex polygons the number is 'n_vertices + 1', and for
574 * a single line segment ('n_vertices == 2' and 'closed == FALSE') the number
575 * is 2.
576 *
577 * 't_vertices' may not alias 'vertices', except when transforming a single
578 * line segment.
579 */
580 void
gimp_transform_polygon(const GimpMatrix3 * matrix,const GimpVector2 * vertices,gint n_vertices,gboolean closed,GimpVector2 * t_vertices,gint * n_t_vertices)581 gimp_transform_polygon (const GimpMatrix3 *matrix,
582 const GimpVector2 *vertices,
583 gint n_vertices,
584 gboolean closed,
585 GimpVector2 *t_vertices,
586 gint *n_t_vertices)
587 {
588 GimpVector3 curr;
589 gboolean curr_visible;
590 gint i;
591
592 g_return_if_fail (matrix != NULL);
593 g_return_if_fail (vertices != NULL);
594 g_return_if_fail (n_vertices >= 0);
595 g_return_if_fail (t_vertices != NULL);
596 g_return_if_fail (n_t_vertices != NULL);
597
598 *n_t_vertices = 0;
599
600 if (n_vertices == 0)
601 return;
602
603 curr.x = matrix->coeff[0][0] * vertices[0].x +
604 matrix->coeff[0][1] * vertices[0].y +
605 matrix->coeff[0][2];
606 curr.y = matrix->coeff[1][0] * vertices[0].x +
607 matrix->coeff[1][1] * vertices[0].y +
608 matrix->coeff[1][2];
609 curr.z = matrix->coeff[2][0] * vertices[0].x +
610 matrix->coeff[2][1] * vertices[0].y +
611 matrix->coeff[2][2];
612
613 curr_visible = (curr.z >= GIMP_TRANSFORM_NEAR_Z);
614
615 for (i = 0; i < n_vertices; i++)
616 {
617 if (curr_visible)
618 {
619 t_vertices[(*n_t_vertices)++] = (GimpVector2) { curr.x / curr.z,
620 curr.y / curr.z };
621 }
622
623 if (i < n_vertices - 1 || closed)
624 {
625 GimpVector3 next;
626 gboolean next_visible;
627 gint j = (i + 1) % n_vertices;
628
629 next.x = matrix->coeff[0][0] * vertices[j].x +
630 matrix->coeff[0][1] * vertices[j].y +
631 matrix->coeff[0][2];
632 next.y = matrix->coeff[1][0] * vertices[j].x +
633 matrix->coeff[1][1] * vertices[j].y +
634 matrix->coeff[1][2];
635 next.z = matrix->coeff[2][0] * vertices[j].x +
636 matrix->coeff[2][1] * vertices[j].y +
637 matrix->coeff[2][2];
638
639 next_visible = (next.z >= GIMP_TRANSFORM_NEAR_Z);
640
641 if (next_visible != curr_visible)
642 {
643 gdouble ratio = (curr.z - GIMP_TRANSFORM_NEAR_Z) / (curr.z - next.z);
644
645 t_vertices[(*n_t_vertices)++] =
646 (GimpVector2) { (curr.x + (next.x - curr.x) * ratio) / GIMP_TRANSFORM_NEAR_Z,
647 (curr.y + (next.y - curr.y) * ratio) / GIMP_TRANSFORM_NEAR_Z };
648 }
649
650 curr = next;
651 curr_visible = next_visible;
652 }
653 }
654 }
655
656 /* same as gimp_transform_polygon(), but using GimpCoords as the vertex type,
657 * instead of GimpVector2.
658 */
659 void
gimp_transform_polygon_coords(const GimpMatrix3 * matrix,const GimpCoords * vertices,gint n_vertices,gboolean closed,GimpCoords * t_vertices,gint * n_t_vertices)660 gimp_transform_polygon_coords (const GimpMatrix3 *matrix,
661 const GimpCoords *vertices,
662 gint n_vertices,
663 gboolean closed,
664 GimpCoords *t_vertices,
665 gint *n_t_vertices)
666 {
667 GimpVector3 curr;
668 gboolean curr_visible;
669 gint i;
670
671 g_return_if_fail (matrix != NULL);
672 g_return_if_fail (vertices != NULL);
673 g_return_if_fail (n_vertices >= 0);
674 g_return_if_fail (t_vertices != NULL);
675 g_return_if_fail (n_t_vertices != NULL);
676
677 *n_t_vertices = 0;
678
679 if (n_vertices == 0)
680 return;
681
682 curr.x = matrix->coeff[0][0] * vertices[0].x +
683 matrix->coeff[0][1] * vertices[0].y +
684 matrix->coeff[0][2];
685 curr.y = matrix->coeff[1][0] * vertices[0].x +
686 matrix->coeff[1][1] * vertices[0].y +
687 matrix->coeff[1][2];
688 curr.z = matrix->coeff[2][0] * vertices[0].x +
689 matrix->coeff[2][1] * vertices[0].y +
690 matrix->coeff[2][2];
691
692 curr_visible = (curr.z >= GIMP_TRANSFORM_NEAR_Z);
693
694 for (i = 0; i < n_vertices; i++)
695 {
696 if (curr_visible)
697 {
698 t_vertices[*n_t_vertices] = vertices[i];
699 t_vertices[*n_t_vertices].x = curr.x / curr.z;
700 t_vertices[*n_t_vertices].y = curr.y / curr.z;
701
702 (*n_t_vertices)++;
703 }
704
705 if (i < n_vertices - 1 || closed)
706 {
707 GimpVector3 next;
708 gboolean next_visible;
709 gint j = (i + 1) % n_vertices;
710
711 next.x = matrix->coeff[0][0] * vertices[j].x +
712 matrix->coeff[0][1] * vertices[j].y +
713 matrix->coeff[0][2];
714 next.y = matrix->coeff[1][0] * vertices[j].x +
715 matrix->coeff[1][1] * vertices[j].y +
716 matrix->coeff[1][2];
717 next.z = matrix->coeff[2][0] * vertices[j].x +
718 matrix->coeff[2][1] * vertices[j].y +
719 matrix->coeff[2][2];
720
721 next_visible = (next.z >= GIMP_TRANSFORM_NEAR_Z);
722
723 if (next_visible != curr_visible)
724 {
725 gdouble ratio = (curr.z - GIMP_TRANSFORM_NEAR_Z) / (curr.z - next.z);
726
727 gimp_coords_mix (1.0 - ratio, &vertices[i],
728 ratio, &vertices[j],
729 &t_vertices[*n_t_vertices]);
730
731 t_vertices[*n_t_vertices].x = (curr.x + (next.x - curr.x) * ratio) /
732 GIMP_TRANSFORM_NEAR_Z;
733 t_vertices[*n_t_vertices].y = (curr.y + (next.y - curr.y) * ratio) /
734 GIMP_TRANSFORM_NEAR_Z;
735
736 (*n_t_vertices)++;
737 }
738
739 curr = next;
740 curr_visible = next_visible;
741 }
742 }
743 }
744
745 /* returns the value of the polynomial 'poly', of degree 'degree', at 'x'. the
746 * coefficients of 'poly' should be specified in descending-degree order.
747 */
748 static gdouble
polynomial_eval(const gdouble * poly,gint degree,gdouble x)749 polynomial_eval (const gdouble *poly,
750 gint degree,
751 gdouble x)
752 {
753 gdouble y = poly[0];
754 gint i;
755
756 for (i = 1; i <= degree; i++)
757 y = y * x + poly[i];
758
759 return y;
760 }
761
762 /* derives the polynomial 'poly', of degree 'degree'.
763 *
764 * returns the derivative in 'result'.
765 */
766 static void
polynomial_derive(const gdouble * poly,gint degree,gdouble * result)767 polynomial_derive (const gdouble *poly,
768 gint degree,
769 gdouble *result)
770 {
771 while (degree)
772 *result++ = *poly++ * degree--;
773 }
774
775 /* finds the real odd-multiplicity root of the polynomial 'poly', of degree
776 * 'degree', inside the range '(x1, x2)'.
777 *
778 * returns TRUE if such a root exists, and stores its value in '*root'.
779 *
780 * 'poly' shall be monotonic in the range '(x1, x2)'.
781 */
782 static gboolean
polynomial_odd_root(const gdouble * poly,gint degree,gdouble x1,gdouble x2,gdouble * root)783 polynomial_odd_root (const gdouble *poly,
784 gint degree,
785 gdouble x1,
786 gdouble x2,
787 gdouble *root)
788 {
789 gdouble y1;
790 gdouble y2;
791 gint i;
792
793 y1 = polynomial_eval (poly, degree, x1);
794 y2 = polynomial_eval (poly, degree, x2);
795
796 if (y1 * y2 > -EPSILON)
797 {
798 /* the two endpoints have the same sign, or one of them is zero. there's
799 * no root inside the range.
800 */
801 return FALSE;
802 }
803 else if (y1 > 0.0)
804 {
805 gdouble t;
806
807 /* if the first endpoint is positive, swap the endpoints, so that the
808 * first endpoint is always negative, and the second endpoint is always
809 * positive.
810 */
811
812 t = x1;
813 x1 = x2;
814 x2 = t;
815 }
816
817 /* approximate the root using binary search */
818 for (i = 0; i < 53; i++)
819 {
820 gdouble x = (x1 + x2) / 2.0;
821 gdouble y = polynomial_eval (poly, degree, x);
822
823 if (y > 0.0)
824 x2 = x;
825 else
826 x1 = x;
827 }
828
829 *root = (x1 + x2) / 2.0;
830
831 return TRUE;
832 }
833
834 /* finds the real odd-multiplicity roots of the polynomial 'poly', of degree
835 * 'degree', inside the range '(x1, x2)'.
836 *
837 * returns the roots in 'roots', in ascending order, and their count in
838 * 'n_roots'.
839 */
840 static void
polynomial_odd_roots(const gdouble * poly,gint degree,gdouble x1,gdouble x2,gdouble * roots,gint * n_roots)841 polynomial_odd_roots (const gdouble *poly,
842 gint degree,
843 gdouble x1,
844 gdouble x2,
845 gdouble *roots,
846 gint *n_roots)
847 {
848 *n_roots = 0;
849
850 /* find the real degree of the polynomial (skip any leading coefficients that
851 * are 0)
852 */
853 for (; degree && fabs (*poly) < EPSILON; poly++, degree--);
854
855 #define ADD_ROOT(root) \
856 do \
857 { \
858 gdouble r = (root); \
859 \
860 if (r > x1 && r < x2) \
861 roots[(*n_roots)++] = r; \
862 } \
863 while (FALSE)
864
865 switch (degree)
866 {
867 /* constant case */
868 case 0:
869 break;
870
871 /* linear case */
872 case 1:
873 ADD_ROOT (-poly[1] / poly[0]);
874 break;
875
876 /* quadratic case */
877 case 2:
878 {
879 gdouble s = SQR (poly[1]) - 4 * poly[0] * poly[2];
880
881 if (s > EPSILON)
882 {
883 s = sqrt (s);
884
885 if (poly[0] < 0.0)
886 s = -s;
887
888 ADD_ROOT ((-poly[1] - s) / (2.0 * poly[0]));
889 ADD_ROOT ((-poly[1] + s) / (2.0 * poly[0]));
890 }
891
892 break;
893 }
894
895 /* general case */
896 default:
897 {
898 gdouble deriv[degree];
899 gdouble deriv_roots[degree - 1];
900 gint n_deriv_roots;
901 gdouble a;
902 gdouble b;
903 gint i;
904
905 /* find the odd roots of the derivative, i.e., the local extrema of the
906 * polynomial
907 */
908 polynomial_derive (poly, degree, deriv);
909 polynomial_odd_roots (deriv, degree - 1, x1, x2,
910 deriv_roots, &n_deriv_roots);
911
912 /* search for roots between each consecutive pair of extrema, including
913 * the endpoints
914 */
915 a = x1;
916
917 for (i = 0; i <= n_deriv_roots; i++)
918 {
919 if (i < n_deriv_roots)
920 b = deriv_roots[i];
921 else
922 b = x2;
923
924 *n_roots += polynomial_odd_root (poly, degree, a, b,
925 &roots[*n_roots]);
926
927 a = b;
928 }
929
930 break;
931 }
932 }
933
934 #undef ADD_ROOT
935 }
936
937 /* clips the cubic bezier segment, defined by the four control points 'bezier',
938 * to the halfplane 'ax + by + c >= 0'.
939 *
940 * returns the clipped set of bezier segments in 'c_bezier', and their count in
941 * 'n_c_bezier'. the minimal possible number of clipped segments is 0, which
942 * happens when the entire segment is clipped. the maximal possible number of
943 * clipped segments is 2.
944 *
945 * if the first clipped segment is an initial segment of 'bezier', sets
946 * '*start_in' to TRUE, otherwise to FALSE. if the last clipped segment is a
947 * final segment of 'bezier', sets '*end_in' to TRUE, otherwise to FALSE.
948 *
949 * 'c_bezier' may not alias 'bezier'.
950 */
951 static void
clip_bezier(const GimpCoords bezier[4],gdouble a,gdouble b,gdouble c,GimpCoords c_bezier[2][4],gint * n_c_bezier,gboolean * start_in,gboolean * end_in)952 clip_bezier (const GimpCoords bezier[4],
953 gdouble a,
954 gdouble b,
955 gdouble c,
956 GimpCoords c_bezier[2][4],
957 gint *n_c_bezier,
958 gboolean *start_in,
959 gboolean *end_in)
960 {
961 gdouble dot[4];
962 gdouble poly[4];
963 gdouble roots[5];
964 gint n_roots;
965 gint n_positive;
966 gint i;
967
968 n_positive = 0;
969
970 for (i = 0; i < 4; i++)
971 {
972 dot[i] = a * bezier[i].x + b * bezier[i].y + c;
973
974 n_positive += (dot[i] >= 0.0);
975 }
976
977 if (n_positive == 0)
978 {
979 /* all points are out -- the entire segment is out */
980
981 *n_c_bezier = 0;
982 *start_in = FALSE;
983 *end_in = FALSE;
984
985 return;
986 }
987 else if (n_positive == 4)
988 {
989 /* all points are in -- the entire segment is in */
990
991 memcpy (c_bezier[0], bezier, sizeof (GimpCoords[4]));
992
993 *n_c_bezier = 1;
994 *start_in = TRUE;
995 *end_in = TRUE;
996
997 return;
998 }
999
1000 /* find the points of intersection of the segment with the 'ax + by + c = 0'
1001 * line
1002 */
1003 poly[0] = dot[3] - 3.0 * dot[2] + 3.0 * dot[1] - dot[0];
1004 poly[1] = 3.0 * (dot[2] - 2.0 * dot[1] + dot[0]);
1005 poly[2] = 3.0 * (dot[1] - dot[0]);
1006 poly[3] = dot[0];
1007
1008 roots[0] = 0.0;
1009 polynomial_odd_roots (poly, 3, 0.0, 1.0, roots + 1, &n_roots);
1010 roots[++n_roots] = 1.0;
1011
1012 /* construct the list of segments that are inside the halfplane */
1013 *n_c_bezier = 0;
1014 *start_in = (polynomial_eval (poly, 3, roots[1] / 2.0) > 0.0);
1015 *end_in = (*start_in + n_roots + 1) % 2;
1016
1017 for (i = ! *start_in; i < n_roots; i += 2)
1018 {
1019 gdouble t0 = roots[i];
1020 gdouble t1 = roots[i + 1];
1021
1022 gimp_coords_interpolate_bezier_at (bezier, t0,
1023 &c_bezier[*n_c_bezier][0],
1024 &c_bezier[*n_c_bezier][1]);
1025 gimp_coords_interpolate_bezier_at (bezier, t1,
1026 &c_bezier[*n_c_bezier][3],
1027 &c_bezier[*n_c_bezier][2]);
1028
1029 gimp_coords_mix (1.0, &c_bezier[*n_c_bezier][0],
1030 (t1 - t0) / 3.0, &c_bezier[*n_c_bezier][1],
1031 &c_bezier[*n_c_bezier][1]);
1032 gimp_coords_mix (1.0, &c_bezier[*n_c_bezier][3],
1033 (t0 - t1) / 3.0, &c_bezier[*n_c_bezier][2],
1034 &c_bezier[*n_c_bezier][2]);
1035
1036 (*n_c_bezier)++;
1037 }
1038 }
1039
1040 /* transforms the cubic bezier segment, defined by the four control points
1041 * 'bezier', by 'matrix', subdividing it as necessary to avoid diverging too
1042 * much from the real transformed curve. at most 'depth' subdivisions are
1043 * performed.
1044 *
1045 * appends the transformed sequence of bezier segments to 't_beziers'.
1046 *
1047 * 'bezier' shall be fully clipped to the near plane.
1048 */
1049 static void
transform_bezier_coords(const GimpMatrix3 * matrix,const GimpCoords bezier[4],GQueue * t_beziers,gint depth)1050 transform_bezier_coords (const GimpMatrix3 *matrix,
1051 const GimpCoords bezier[4],
1052 GQueue *t_beziers,
1053 gint depth)
1054 {
1055 GimpCoords *t_bezier;
1056 gint n;
1057
1058 /* check if we need to split the segment */
1059 if (depth > 0)
1060 {
1061 GimpVector2 v[4];
1062 GimpVector2 c[2];
1063 GimpVector2 b;
1064 gint i;
1065
1066 for (i = 0; i < 4; i++)
1067 v[i] = (GimpVector2) { bezier[i].x, bezier[i].y };
1068
1069 gimp_vector2_sub (&c[0], &v[1], &v[0]);
1070 gimp_vector2_sub (&c[1], &v[2], &v[3]);
1071
1072 gimp_vector2_sub (&b, &v[3], &v[0]);
1073 gimp_vector2_mul (&b, 1.0 / gimp_vector2_inner_product (&b, &b));
1074
1075 for (i = 0; i < 2; i++)
1076 {
1077 /* split the segment if one of the control points is too far from the
1078 * line connecting the anchors
1079 */
1080 if (fabs (gimp_vector2_cross_product (&c[i], &b).x) > 0.5)
1081 {
1082 GimpCoords mid_position;
1083 GimpCoords mid_velocity;
1084 GimpCoords sub[4];
1085
1086 gimp_coords_interpolate_bezier_at (bezier, 0.5,
1087 &mid_position, &mid_velocity);
1088
1089 /* first half */
1090 sub[0] = bezier[0];
1091 sub[3] = mid_position;
1092
1093 gimp_coords_mix (0.5, &sub[0],
1094 0.5, &bezier[1],
1095 &sub[1]);
1096 gimp_coords_mix (1.0, &sub[3],
1097 -1.0 / 6.0, &mid_velocity,
1098 &sub[2]);
1099
1100 transform_bezier_coords (matrix, sub, t_beziers, depth - 1);
1101
1102 /* second half */
1103 sub[0] = mid_position;
1104 sub[3] = bezier[3];
1105
1106 gimp_coords_mix (1.0, &sub[0],
1107 +1.0 / 6.0, &mid_velocity,
1108 &sub[1]);
1109 gimp_coords_mix (0.5, &sub[3],
1110 0.5, &bezier[2],
1111 &sub[2]);
1112
1113 transform_bezier_coords (matrix, sub, t_beziers, depth - 1);
1114
1115 return;
1116 }
1117 }
1118 }
1119
1120 /* transform the segment by transforming each of the individual points. note
1121 * that, for non-affine transforms, this is only an approximation of the real
1122 * transformed curve, but due to subdivision it should be good enough.
1123 */
1124 t_bezier = g_new (GimpCoords, 4);
1125
1126 /* note that while the segments themselves are clipped to the near plane,
1127 * their control points may still get transformed behind the camera. we
1128 * therefore clip the control points to the near plane as well, which is not
1129 * too meaningful, but avoids erroneously transforming them behind the
1130 * camera.
1131 */
1132 gimp_transform_polygon_coords (matrix, bezier, 2, FALSE,
1133 t_bezier, &n);
1134 gimp_transform_polygon_coords (matrix, bezier + 2, 2, FALSE,
1135 t_bezier + 2, &n);
1136
1137 g_queue_push_tail (t_beziers, t_bezier);
1138 }
1139
1140 /* transforms the cubic bezier segment, defined by the four control points
1141 * 'bezier', by 'matrix', performing clipping by the near plane and subdividing
1142 * as necessary.
1143 *
1144 * returns the transformed set of bezier-segment sequences in 't_beziers', as
1145 * GQueues of GimpCoords[4] bezier-segments, and the number of sequences in
1146 * 'n_t_beziers'. the minimal possible number of transformed sequences is 0,
1147 * which happens when the entire segment is clipped. the maximal possible
1148 * number of transformed sequences is 2. each sequence has at least one
1149 * segment.
1150 *
1151 * if the first transformed segment is an initial segment of 'bezier', sets
1152 * '*start_in' to TRUE, otherwise to FALSE. if the last transformed segment is
1153 * a final segment of 'bezier', sets '*end_in' to TRUE, otherwise to FALSE.
1154 */
1155 void
gimp_transform_bezier_coords(const GimpMatrix3 * matrix,const GimpCoords bezier[4],GQueue * t_beziers[2],gint * n_t_beziers,gboolean * start_in,gboolean * end_in)1156 gimp_transform_bezier_coords (const GimpMatrix3 *matrix,
1157 const GimpCoords bezier[4],
1158 GQueue *t_beziers[2],
1159 gint *n_t_beziers,
1160 gboolean *start_in,
1161 gboolean *end_in)
1162 {
1163 GimpCoords c_bezier[2][4];
1164 gint i;
1165
1166 g_return_if_fail (matrix != NULL);
1167 g_return_if_fail (bezier != NULL);
1168 g_return_if_fail (t_beziers != NULL);
1169 g_return_if_fail (n_t_beziers != NULL);
1170 g_return_if_fail (start_in != NULL);
1171 g_return_if_fail (end_in != NULL);
1172
1173 /* if the matrix is affine, transform the easy way */
1174 if (gimp_matrix3_is_affine (matrix))
1175 {
1176 GimpCoords *t_bezier;
1177
1178 t_beziers[0] = g_queue_new ();
1179 *n_t_beziers = 1;
1180
1181 t_bezier = g_new (GimpCoords, 1);
1182 g_queue_push_tail (t_beziers[0], t_bezier);
1183
1184 for (i = 0; i < 4; i++)
1185 {
1186 t_bezier[i] = bezier[i];
1187
1188 gimp_matrix3_transform_point (matrix,
1189 bezier[i].x, bezier[i].y,
1190 &t_bezier[i].x, &t_bezier[i].y);
1191 }
1192
1193 return;
1194 }
1195
1196 /* clip the segment to the near plane */
1197 clip_bezier (bezier,
1198 matrix->coeff[2][0],
1199 matrix->coeff[2][1],
1200 matrix->coeff[2][2] - GIMP_TRANSFORM_NEAR_Z,
1201 c_bezier, n_t_beziers,
1202 start_in, end_in);
1203
1204 /* transform each of the resulting segments */
1205 for (i = 0; i < *n_t_beziers; i++)
1206 {
1207 t_beziers[i] = g_queue_new ();
1208
1209 transform_bezier_coords (matrix, c_bezier[i], t_beziers[i], 3);
1210 }
1211 }
1212