diff options
Diffstat (limited to 'app/core/gimp-transform-utils.c')
| -rw-r--r-- | app/core/gimp-transform-utils.c | 1211 |
1 files changed, 1211 insertions, 0 deletions
diff --git a/app/core/gimp-transform-utils.c b/app/core/gimp-transform-utils.c new file mode 100644 index 0000000..555ff09 --- /dev/null +++ b/app/core/gimp-transform-utils.c @@ -0,0 +1,1211 @@ +/* GIMP - The GNU Image Manipulation Program + * Copyright (C) 1995-2001 Spencer Kimball, Peter Mattis, and others + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or + * (at your option) any later version. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <https://www.gnu.org/licenses/>. + */ + +#include "config.h" + +#include <string.h> + +#include <glib-object.h> + +#include "libgimpmath/gimpmath.h" + +#include "core-types.h" + +#include "gimp-transform-utils.h" +#include "gimpcoords.h" +#include "gimpcoords-interpolate.h" + + +#define EPSILON 1e-6 + + +void +gimp_transform_get_rotate_center (gint x, + gint y, + gint width, + gint height, + gboolean auto_center, + gdouble *center_x, + gdouble *center_y) +{ + g_return_if_fail (center_x != NULL); + g_return_if_fail (center_y != NULL); + + if (auto_center) + { + *center_x = (gdouble) x + (gdouble) width / 2.0; + *center_y = (gdouble) y + (gdouble) height / 2.0; + } +} + +void +gimp_transform_get_flip_axis (gint x, + gint y, + gint width, + gint height, + GimpOrientationType flip_type, + gboolean auto_center, + gdouble *axis) +{ + g_return_if_fail (axis != NULL); + + if (auto_center) + { + switch (flip_type) + { + case GIMP_ORIENTATION_HORIZONTAL: + *axis = ((gdouble) x + (gdouble) width / 2.0); + break; + + case GIMP_ORIENTATION_VERTICAL: + *axis = ((gdouble) y + (gdouble) height / 2.0); + break; + + default: + g_return_if_reached (); + break; + } + } +} + +void +gimp_transform_matrix_flip (GimpMatrix3 *matrix, + GimpOrientationType flip_type, + gdouble axis) +{ + g_return_if_fail (matrix != NULL); + + switch (flip_type) + { + case GIMP_ORIENTATION_HORIZONTAL: + gimp_matrix3_translate (matrix, - axis, 0.0); + gimp_matrix3_scale (matrix, -1.0, 1.0); + gimp_matrix3_translate (matrix, axis, 0.0); + break; + + case GIMP_ORIENTATION_VERTICAL: + gimp_matrix3_translate (matrix, 0.0, - axis); + gimp_matrix3_scale (matrix, 1.0, -1.0); + gimp_matrix3_translate (matrix, 0.0, axis); + break; + + case GIMP_ORIENTATION_UNKNOWN: + break; + } +} + +void +gimp_transform_matrix_flip_free (GimpMatrix3 *matrix, + gdouble x1, + gdouble y1, + gdouble x2, + gdouble y2) +{ + gdouble angle; + + g_return_if_fail (matrix != NULL); + + angle = atan2 (y2 - y1, x2 - x1); + + gimp_matrix3_identity (matrix); + gimp_matrix3_translate (matrix, -x1, -y1); + gimp_matrix3_rotate (matrix, -angle); + gimp_matrix3_scale (matrix, 1.0, -1.0); + gimp_matrix3_rotate (matrix, angle); + gimp_matrix3_translate (matrix, x1, y1); +} + +void +gimp_transform_matrix_rotate (GimpMatrix3 *matrix, + GimpRotationType rotate_type, + gdouble center_x, + gdouble center_y) +{ + gdouble angle = 0; + + switch (rotate_type) + { + case GIMP_ROTATE_90: + angle = G_PI_2; + break; + case GIMP_ROTATE_180: + angle = G_PI; + break; + case GIMP_ROTATE_270: + angle = - G_PI_2; + break; + } + + gimp_transform_matrix_rotate_center (matrix, center_x, center_y, angle); +} + +void +gimp_transform_matrix_rotate_rect (GimpMatrix3 *matrix, + gint x, + gint y, + gint width, + gint height, + gdouble angle) +{ + gdouble center_x; + gdouble center_y; + + g_return_if_fail (matrix != NULL); + + center_x = (gdouble) x + (gdouble) width / 2.0; + center_y = (gdouble) y + (gdouble) height / 2.0; + + gimp_matrix3_translate (matrix, -center_x, -center_y); + gimp_matrix3_rotate (matrix, angle); + gimp_matrix3_translate (matrix, +center_x, +center_y); +} + +void +gimp_transform_matrix_rotate_center (GimpMatrix3 *matrix, + gdouble center_x, + gdouble center_y, + gdouble angle) +{ + g_return_if_fail (matrix != NULL); + + gimp_matrix3_translate (matrix, -center_x, -center_y); + gimp_matrix3_rotate (matrix, angle); + gimp_matrix3_translate (matrix, +center_x, +center_y); +} + +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) +{ + gdouble scale_x = 1.0; + gdouble scale_y = 1.0; + + g_return_if_fail (matrix != NULL); + + if (width > 0) + scale_x = t_width / (gdouble) width; + + if (height > 0) + scale_y = t_height / (gdouble) height; + + gimp_matrix3_identity (matrix); + gimp_matrix3_translate (matrix, -x, -y); + gimp_matrix3_scale (matrix, scale_x, scale_y); + gimp_matrix3_translate (matrix, t_x, t_y); +} + +void +gimp_transform_matrix_shear (GimpMatrix3 *matrix, + gint x, + gint y, + gint width, + gint height, + GimpOrientationType orientation, + gdouble amount) +{ + gdouble center_x; + gdouble center_y; + + g_return_if_fail (matrix != NULL); + + if (width == 0) + width = 1; + + if (height == 0) + height = 1; + + center_x = (gdouble) x + (gdouble) width / 2.0; + center_y = (gdouble) y + (gdouble) height / 2.0; + + gimp_matrix3_identity (matrix); + gimp_matrix3_translate (matrix, -center_x, -center_y); + + if (orientation == GIMP_ORIENTATION_HORIZONTAL) + gimp_matrix3_xshear (matrix, amount / height); + else + gimp_matrix3_yshear (matrix, amount / width); + + gimp_matrix3_translate (matrix, +center_x, +center_y); +} + +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) +{ + GimpMatrix3 trafo; + gdouble scalex; + gdouble scaley; + + g_return_if_fail (matrix != NULL); + + scalex = scaley = 1.0; + + if (width > 0) + scalex = 1.0 / (gdouble) width; + + if (height > 0) + scaley = 1.0 / (gdouble) height; + + gimp_matrix3_translate (matrix, -x, -y); + gimp_matrix3_scale (matrix, scalex, scaley); + + /* Determine the perspective transform that maps from + * the unit cube to the transformed coordinates + */ + { + gdouble dx1, dx2, dx3, dy1, dy2, dy3; + + dx1 = t_x2 - t_x4; + dx2 = t_x3 - t_x4; + dx3 = t_x1 - t_x2 + t_x4 - t_x3; + + dy1 = t_y2 - t_y4; + dy2 = t_y3 - t_y4; + dy3 = t_y1 - t_y2 + t_y4 - t_y3; + + /* Is the mapping affine? */ + if ((dx3 == 0.0) && (dy3 == 0.0)) + { + trafo.coeff[0][0] = t_x2 - t_x1; + trafo.coeff[0][1] = t_x4 - t_x2; + trafo.coeff[0][2] = t_x1; + trafo.coeff[1][0] = t_y2 - t_y1; + trafo.coeff[1][1] = t_y4 - t_y2; + trafo.coeff[1][2] = t_y1; + trafo.coeff[2][0] = 0.0; + trafo.coeff[2][1] = 0.0; + } + else + { + gdouble det1, det2; + + det1 = dx3 * dy2 - dy3 * dx2; + det2 = dx1 * dy2 - dy1 * dx2; + + trafo.coeff[2][0] = (det2 == 0.0) ? 1.0 : det1 / det2; + + det1 = dx1 * dy3 - dy1 * dx3; + + trafo.coeff[2][1] = (det2 == 0.0) ? 1.0 : det1 / det2; + + trafo.coeff[0][0] = t_x2 - t_x1 + trafo.coeff[2][0] * t_x2; + trafo.coeff[0][1] = t_x3 - t_x1 + trafo.coeff[2][1] * t_x3; + trafo.coeff[0][2] = t_x1; + + trafo.coeff[1][0] = t_y2 - t_y1 + trafo.coeff[2][0] * t_y2; + trafo.coeff[1][1] = t_y3 - t_y1 + trafo.coeff[2][1] * t_y3; + trafo.coeff[1][2] = t_y1; + } + + trafo.coeff[2][2] = 1.0; + } + + gimp_matrix3_mult (&trafo, matrix); +} + +/* modified gaussian algorithm + * solves a system of linear equations + * + * Example: + * 1x + 2y + 4z = 25 + * 2x + 1y = 4 + * 3x + 5y + 2z = 23 + * Solution: x=1, y=2, z=5 + * + * Input: + * matrix = { 1,2,4,25,2,1,0,4,3,5,2,23 } + * s = 3 (Number of variables) + * Output: + * return value == TRUE (TRUE, if there is a single unique solution) + * solution == { 1,2,5 } (if the return value is FALSE, the content + * of solution is of no use) + */ +static gboolean +mod_gauss (gdouble matrix[], + gdouble solution[], + gint s) +{ + gint p[s]; /* row permutation */ + gint i, j, r, temp; + gdouble q; + gint t = s + 1; + + for (i = 0; i < s; i++) + { + p[i] = i; + } + + for (r = 0; r < s; r++) + { + /* make sure that (r,r) is not 0 */ + if (fabs (matrix[p[r] * t + r]) <= EPSILON) + { + /* we need to permutate rows */ + for (i = r + 1; i <= s; i++) + { + if (i == s) + { + /* if this happens, the linear system has zero or + * more than one solutions. + */ + return FALSE; + } + + if (fabs (matrix[p[i] * t + r]) > EPSILON) + break; + } + + temp = p[r]; + p[r] = p[i]; + p[i] = temp; + } + + /* make (r,r) == 1 */ + q = 1.0 / matrix[p[r] * t + r]; + matrix[p[r] * t + r] = 1.0; + + for (j = r + 1; j < t; j++) + { + matrix[p[r] * t + j] *= q; + } + + /* make that all entries in column r are 0 (except (r,r)) */ + for (i = 0; i < s; i++) + { + if (i == r) + continue; + + for (j = r + 1; j < t ; j++) + { + matrix[p[i] * t + j] -= matrix[p[r] * t + j] * matrix[p[i] * t + r]; + } + + /* we don't need to execute the following line + * since we won't access this element again: + * + * matrix[p[i] * t + r] = 0.0; + */ + } + } + + for (i = 0; i < s; i++) + { + solution[i] = matrix[p[i] * t + s]; + } + + return TRUE; +} + +/* multiplies 'matrix' by the matrix that transforms a set of 4 'input_points' + * to corresponding 'output_points', if such matrix exists, and is valid (i.e., + * keeps the output points in front of the camera). + * + * returns TRUE if successful. + */ +gboolean +gimp_transform_matrix_generic (GimpMatrix3 *matrix, + const GimpVector2 input_points[4], + const GimpVector2 output_points[4]) +{ + GimpMatrix3 trafo; + gdouble coeff[8 * 9]; + gboolean negative = -1; + gint i; + gboolean result = TRUE; + + g_return_val_if_fail (matrix != NULL, FALSE); + g_return_val_if_fail (input_points != NULL, FALSE); + g_return_val_if_fail (output_points != NULL, FALSE); + + /* find the matrix that transforms 'input_points' to 'output_points', whose + * (3, 3) coeffcient is 1, by solving a system of linear equations whose + * solution is the remaining 8 coefficients. + */ + for (i = 0; i < 4; i++) + { + coeff[i * 9 + 0] = input_points[i].x; + coeff[i * 9 + 1] = input_points[i].y; + coeff[i * 9 + 2] = 1.0; + coeff[i * 9 + 3] = 0.0; + coeff[i * 9 + 4] = 0.0; + coeff[i * 9 + 5] = 0.0; + coeff[i * 9 + 6] = -input_points[i].x * output_points[i].x; + coeff[i * 9 + 7] = -input_points[i].y * output_points[i].x; + coeff[i * 9 + 8] = output_points[i].x; + + coeff[(i + 4) * 9 + 0] = 0.0; + coeff[(i + 4) * 9 + 1] = 0.0; + coeff[(i + 4) * 9 + 2] = 0.0; + coeff[(i + 4) * 9 + 3] = input_points[i].x; + coeff[(i + 4) * 9 + 4] = input_points[i].y; + coeff[(i + 4) * 9 + 5] = 1.0; + coeff[(i + 4) * 9 + 6] = -input_points[i].x * output_points[i].y; + coeff[(i + 4) * 9 + 7] = -input_points[i].y * output_points[i].y; + coeff[(i + 4) * 9 + 8] = output_points[i].y; + } + + /* if there is no solution, bail */ + if (! mod_gauss (coeff, (gdouble *) trafo.coeff, 8)) + return FALSE; + + trafo.coeff[2][2] = 1.0; + + /* make sure that none of the input points maps to a point at infinity, and + * that all output points are on the same side of the camera. + */ + for (i = 0; i < 4; i++) + { + gdouble w; + gboolean neg; + + w = trafo.coeff[2][0] * input_points[i].x + + trafo.coeff[2][1] * input_points[i].y + + trafo.coeff[2][2]; + + if (fabs (w) <= EPSILON) + result = FALSE; + + neg = (w < 0.0); + + if (negative < 0) + { + negative = neg; + } + else if (neg != negative) + { + result = FALSE; + break; + } + } + + /* if the output points are all behind the camera, negate the matrix, which + * would map the input points to the corresponding points in front of the + * camera. + */ + if (negative > 0) + { + gint r; + gint c; + + for (r = 0; r < 3; r++) + { + for (c = 0; c < 3; c++) + { + trafo.coeff[r][c] = -trafo.coeff[r][c]; + } + } + } + + /* append the transformation to 'matrix' */ + gimp_matrix3_mult (&trafo, matrix); + + return result; +} + +gboolean +gimp_transform_polygon_is_convex (gdouble x1, + gdouble y1, + gdouble x2, + gdouble y2, + gdouble x3, + gdouble y3, + gdouble x4, + gdouble y4) +{ + gdouble z1, z2, z3, z4; + + /* We test if the transformed polygon is convex. if z1 and z2 have + * the same sign as well as z3 and z4 the polygon is convex. + */ + z1 = ((x2 - x1) * (y4 - y1) - + (x4 - x1) * (y2 - y1)); + z2 = ((x4 - x1) * (y3 - y1) - + (x3 - x1) * (y4 - y1)); + z3 = ((x4 - x2) * (y3 - y2) - + (x3 - x2) * (y4 - y2)); + z4 = ((x3 - x2) * (y1 - y2) - + (x1 - x2) * (y3 - y2)); + + return (z1 * z2 > 0) && (z3 * z4 > 0); +} + +/* transforms the polygon or polyline, whose vertices are given by 'vertices', + * by 'matrix', performing clipping by the near plane. 'closed' indicates + * whether the vertices represent a polygon ('closed == TRUE') or a polyline + * ('closed == FALSE'). + * + * returns the transformed vertices in 't_vertices', and their count in + * 'n_t_vertices'. the minimal possible number of transformed vertices is 0, + * which happens when the entire input is clipped. in general, the maximal + * possible number of transformed vertices is '3 * n_vertices / 2' (rounded + * down), however, for convex polygons the number is 'n_vertices + 1', and for + * a single line segment ('n_vertices == 2' and 'closed == FALSE') the number + * is 2. + * + * 't_vertices' may not alias 'vertices', except when transforming a single + * line segment. + */ +void +gimp_transform_polygon (const GimpMatrix3 *matrix, + const GimpVector2 *vertices, + gint n_vertices, + gboolean closed, + GimpVector2 *t_vertices, + gint *n_t_vertices) +{ + GimpVector3 curr; + gboolean curr_visible; + gint i; + + g_return_if_fail (matrix != NULL); + g_return_if_fail (vertices != NULL); + g_return_if_fail (n_vertices >= 0); + g_return_if_fail (t_vertices != NULL); + g_return_if_fail (n_t_vertices != NULL); + + *n_t_vertices = 0; + + if (n_vertices == 0) + return; + + curr.x = matrix->coeff[0][0] * vertices[0].x + + matrix->coeff[0][1] * vertices[0].y + + matrix->coeff[0][2]; + curr.y = matrix->coeff[1][0] * vertices[0].x + + matrix->coeff[1][1] * vertices[0].y + + matrix->coeff[1][2]; + curr.z = matrix->coeff[2][0] * vertices[0].x + + matrix->coeff[2][1] * vertices[0].y + + matrix->coeff[2][2]; + + curr_visible = (curr.z >= GIMP_TRANSFORM_NEAR_Z); + + for (i = 0; i < n_vertices; i++) + { + if (curr_visible) + { + t_vertices[(*n_t_vertices)++] = (GimpVector2) { curr.x / curr.z, + curr.y / curr.z }; + } + + if (i < n_vertices - 1 || closed) + { + GimpVector3 next; + gboolean next_visible; + gint j = (i + 1) % n_vertices; + + next.x = matrix->coeff[0][0] * vertices[j].x + + matrix->coeff[0][1] * vertices[j].y + + matrix->coeff[0][2]; + next.y = matrix->coeff[1][0] * vertices[j].x + + matrix->coeff[1][1] * vertices[j].y + + matrix->coeff[1][2]; + next.z = matrix->coeff[2][0] * vertices[j].x + + matrix->coeff[2][1] * vertices[j].y + + matrix->coeff[2][2]; + + next_visible = (next.z >= GIMP_TRANSFORM_NEAR_Z); + + if (next_visible != curr_visible) + { + gdouble ratio = (curr.z - GIMP_TRANSFORM_NEAR_Z) / (curr.z - next.z); + + t_vertices[(*n_t_vertices)++] = + (GimpVector2) { (curr.x + (next.x - curr.x) * ratio) / GIMP_TRANSFORM_NEAR_Z, + (curr.y + (next.y - curr.y) * ratio) / GIMP_TRANSFORM_NEAR_Z }; + } + + curr = next; + curr_visible = next_visible; + } + } +} + +/* same as gimp_transform_polygon(), but using GimpCoords as the vertex type, + * instead of GimpVector2. + */ +void +gimp_transform_polygon_coords (const GimpMatrix3 *matrix, + const GimpCoords *vertices, + gint n_vertices, + gboolean closed, + GimpCoords *t_vertices, + gint *n_t_vertices) +{ + GimpVector3 curr; + gboolean curr_visible; + gint i; + + g_return_if_fail (matrix != NULL); + g_return_if_fail (vertices != NULL); + g_return_if_fail (n_vertices >= 0); + g_return_if_fail (t_vertices != NULL); + g_return_if_fail (n_t_vertices != NULL); + + *n_t_vertices = 0; + + if (n_vertices == 0) + return; + + curr.x = matrix->coeff[0][0] * vertices[0].x + + matrix->coeff[0][1] * vertices[0].y + + matrix->coeff[0][2]; + curr.y = matrix->coeff[1][0] * vertices[0].x + + matrix->coeff[1][1] * vertices[0].y + + matrix->coeff[1][2]; + curr.z = matrix->coeff[2][0] * vertices[0].x + + matrix->coeff[2][1] * vertices[0].y + + matrix->coeff[2][2]; + + curr_visible = (curr.z >= GIMP_TRANSFORM_NEAR_Z); + + for (i = 0; i < n_vertices; i++) + { + if (curr_visible) + { + t_vertices[*n_t_vertices] = vertices[i]; + t_vertices[*n_t_vertices].x = curr.x / curr.z; + t_vertices[*n_t_vertices].y = curr.y / curr.z; + + (*n_t_vertices)++; + } + + if (i < n_vertices - 1 || closed) + { + GimpVector3 next; + gboolean next_visible; + gint j = (i + 1) % n_vertices; + + next.x = matrix->coeff[0][0] * vertices[j].x + + matrix->coeff[0][1] * vertices[j].y + + matrix->coeff[0][2]; + next.y = matrix->coeff[1][0] * vertices[j].x + + matrix->coeff[1][1] * vertices[j].y + + matrix->coeff[1][2]; + next.z = matrix->coeff[2][0] * vertices[j].x + + matrix->coeff[2][1] * vertices[j].y + + matrix->coeff[2][2]; + + next_visible = (next.z >= GIMP_TRANSFORM_NEAR_Z); + + if (next_visible != curr_visible) + { + gdouble ratio = (curr.z - GIMP_TRANSFORM_NEAR_Z) / (curr.z - next.z); + + gimp_coords_mix (1.0 - ratio, &vertices[i], + ratio, &vertices[j], + &t_vertices[*n_t_vertices]); + + t_vertices[*n_t_vertices].x = (curr.x + (next.x - curr.x) * ratio) / + GIMP_TRANSFORM_NEAR_Z; + t_vertices[*n_t_vertices].y = (curr.y + (next.y - curr.y) * ratio) / + GIMP_TRANSFORM_NEAR_Z; + + (*n_t_vertices)++; + } + + curr = next; + curr_visible = next_visible; + } + } +} + +/* returns the value of the polynomial 'poly', of degree 'degree', at 'x'. the + * coefficients of 'poly' should be specified in descending-degree order. + */ +static gdouble +polynomial_eval (const gdouble *poly, + gint degree, + gdouble x) +{ + gdouble y = poly[0]; + gint i; + + for (i = 1; i <= degree; i++) + y = y * x + poly[i]; + + return y; +} + +/* derives the polynomial 'poly', of degree 'degree'. + * + * returns the derivative in 'result'. + */ +static void +polynomial_derive (const gdouble *poly, + gint degree, + gdouble *result) +{ + while (degree) + *result++ = *poly++ * degree--; +} + +/* finds the real odd-multiplicity root of the polynomial 'poly', of degree + * 'degree', inside the range '(x1, x2)'. + * + * returns TRUE if such a root exists, and stores its value in '*root'. + * + * 'poly' shall be monotonic in the range '(x1, x2)'. + */ +static gboolean +polynomial_odd_root (const gdouble *poly, + gint degree, + gdouble x1, + gdouble x2, + gdouble *root) +{ + gdouble y1; + gdouble y2; + gint i; + + y1 = polynomial_eval (poly, degree, x1); + y2 = polynomial_eval (poly, degree, x2); + + if (y1 * y2 > -EPSILON) + { + /* the two endpoints have the same sign, or one of them is zero. there's + * no root inside the range. + */ + return FALSE; + } + else if (y1 > 0.0) + { + gdouble t; + + /* if the first endpoint is positive, swap the endpoints, so that the + * first endpoint is always negative, and the second endpoint is always + * positive. + */ + + t = x1; + x1 = x2; + x2 = t; + } + + /* approximate the root using binary search */ + for (i = 0; i < 53; i++) + { + gdouble x = (x1 + x2) / 2.0; + gdouble y = polynomial_eval (poly, degree, x); + + if (y > 0.0) + x2 = x; + else + x1 = x; + } + + *root = (x1 + x2) / 2.0; + + return TRUE; +} + +/* finds the real odd-multiplicity roots of the polynomial 'poly', of degree + * 'degree', inside the range '(x1, x2)'. + * + * returns the roots in 'roots', in ascending order, and their count in + * 'n_roots'. + */ +static void +polynomial_odd_roots (const gdouble *poly, + gint degree, + gdouble x1, + gdouble x2, + gdouble *roots, + gint *n_roots) +{ + *n_roots = 0; + + /* find the real degree of the polynomial (skip any leading coefficients that + * are 0) + */ + for (; degree && fabs (*poly) < EPSILON; poly++, degree--); + + #define ADD_ROOT(root) \ + do \ + { \ + gdouble r = (root); \ + \ + if (r > x1 && r < x2) \ + roots[(*n_roots)++] = r; \ + } \ + while (FALSE) + + switch (degree) + { + /* constant case */ + case 0: + break; + + /* linear case */ + case 1: + ADD_ROOT (-poly[1] / poly[0]); + break; + + /* quadratic case */ + case 2: + { + gdouble s = SQR (poly[1]) - 4 * poly[0] * poly[2]; + + if (s > EPSILON) + { + s = sqrt (s); + + if (poly[0] < 0.0) + s = -s; + + ADD_ROOT ((-poly[1] - s) / (2.0 * poly[0])); + ADD_ROOT ((-poly[1] + s) / (2.0 * poly[0])); + } + + break; + } + + /* general case */ + default: + { + gdouble deriv[degree]; + gdouble deriv_roots[degree - 1]; + gint n_deriv_roots; + gdouble a; + gdouble b; + gint i; + + /* find the odd roots of the derivative, i.e., the local extrema of the + * polynomial + */ + polynomial_derive (poly, degree, deriv); + polynomial_odd_roots (deriv, degree - 1, x1, x2, + deriv_roots, &n_deriv_roots); + + /* search for roots between each consecutive pair of extrema, including + * the endpoints + */ + a = x1; + + for (i = 0; i <= n_deriv_roots; i++) + { + if (i < n_deriv_roots) + b = deriv_roots[i]; + else + b = x2; + + *n_roots += polynomial_odd_root (poly, degree, a, b, + &roots[*n_roots]); + + a = b; + } + + break; + } + } + + #undef ADD_ROOT +} + +/* clips the cubic bezier segment, defined by the four control points 'bezier', + * to the halfplane 'ax + by + c >= 0'. + * + * returns the clipped set of bezier segments in 'c_bezier', and their count in + * 'n_c_bezier'. the minimal possible number of clipped segments is 0, which + * happens when the entire segment is clipped. the maximal possible number of + * clipped segments is 2. + * + * if the first clipped segment is an initial segment of 'bezier', sets + * '*start_in' to TRUE, otherwise to FALSE. if the last clipped segment is a + * final segment of 'bezier', sets '*end_in' to TRUE, otherwise to FALSE. + * + * 'c_bezier' may not alias 'bezier'. + */ +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) +{ + gdouble dot[4]; + gdouble poly[4]; + gdouble roots[5]; + gint n_roots; + gint n_positive; + gint i; + + n_positive = 0; + + for (i = 0; i < 4; i++) + { + dot[i] = a * bezier[i].x + b * bezier[i].y + c; + + n_positive += (dot[i] >= 0.0); + } + + if (n_positive == 0) + { + /* all points are out -- the entire segment is out */ + + *n_c_bezier = 0; + *start_in = FALSE; + *end_in = FALSE; + + return; + } + else if (n_positive == 4) + { + /* all points are in -- the entire segment is in */ + + memcpy (c_bezier[0], bezier, sizeof (GimpCoords[4])); + + *n_c_bezier = 1; + *start_in = TRUE; + *end_in = TRUE; + + return; + } + + /* find the points of intersection of the segment with the 'ax + by + c = 0' + * line + */ + poly[0] = dot[3] - 3.0 * dot[2] + 3.0 * dot[1] - dot[0]; + poly[1] = 3.0 * (dot[2] - 2.0 * dot[1] + dot[0]); + poly[2] = 3.0 * (dot[1] - dot[0]); + poly[3] = dot[0]; + + roots[0] = 0.0; + polynomial_odd_roots (poly, 3, 0.0, 1.0, roots + 1, &n_roots); + roots[++n_roots] = 1.0; + + /* construct the list of segments that are inside the halfplane */ + *n_c_bezier = 0; + *start_in = (polynomial_eval (poly, 3, roots[1] / 2.0) > 0.0); + *end_in = (*start_in + n_roots + 1) % 2; + + for (i = ! *start_in; i < n_roots; i += 2) + { + gdouble t0 = roots[i]; + gdouble t1 = roots[i + 1]; + + gimp_coords_interpolate_bezier_at (bezier, t0, + &c_bezier[*n_c_bezier][0], + &c_bezier[*n_c_bezier][1]); + gimp_coords_interpolate_bezier_at (bezier, t1, + &c_bezier[*n_c_bezier][3], + &c_bezier[*n_c_bezier][2]); + + gimp_coords_mix (1.0, &c_bezier[*n_c_bezier][0], + (t1 - t0) / 3.0, &c_bezier[*n_c_bezier][1], + &c_bezier[*n_c_bezier][1]); + gimp_coords_mix (1.0, &c_bezier[*n_c_bezier][3], + (t0 - t1) / 3.0, &c_bezier[*n_c_bezier][2], + &c_bezier[*n_c_bezier][2]); + + (*n_c_bezier)++; + } +} + +/* transforms the cubic bezier segment, defined by the four control points + * 'bezier', by 'matrix', subdividing it as necessary to avoid diverging too + * much from the real transformed curve. at most 'depth' subdivisions are + * performed. + * + * appends the transformed sequence of bezier segments to 't_beziers'. + * + * 'bezier' shall be fully clipped to the near plane. + */ +static void +transform_bezier_coords (const GimpMatrix3 *matrix, + const GimpCoords bezier[4], + GQueue *t_beziers, + gint depth) +{ + GimpCoords *t_bezier; + gint n; + + /* check if we need to split the segment */ + if (depth > 0) + { + GimpVector2 v[4]; + GimpVector2 c[2]; + GimpVector2 b; + gint i; + + for (i = 0; i < 4; i++) + v[i] = (GimpVector2) { bezier[i].x, bezier[i].y }; + + gimp_vector2_sub (&c[0], &v[1], &v[0]); + gimp_vector2_sub (&c[1], &v[2], &v[3]); + + gimp_vector2_sub (&b, &v[3], &v[0]); + gimp_vector2_mul (&b, 1.0 / gimp_vector2_inner_product (&b, &b)); + + for (i = 0; i < 2; i++) + { + /* split the segment if one of the control points is too far from the + * line connecting the anchors + */ + if (fabs (gimp_vector2_cross_product (&c[i], &b).x) > 0.5) + { + GimpCoords mid_position; + GimpCoords mid_velocity; + GimpCoords sub[4]; + + gimp_coords_interpolate_bezier_at (bezier, 0.5, + &mid_position, &mid_velocity); + + /* first half */ + sub[0] = bezier[0]; + sub[3] = mid_position; + + gimp_coords_mix (0.5, &sub[0], + 0.5, &bezier[1], + &sub[1]); + gimp_coords_mix (1.0, &sub[3], + -1.0 / 6.0, &mid_velocity, + &sub[2]); + + transform_bezier_coords (matrix, sub, t_beziers, depth - 1); + + /* second half */ + sub[0] = mid_position; + sub[3] = bezier[3]; + + gimp_coords_mix (1.0, &sub[0], + +1.0 / 6.0, &mid_velocity, + &sub[1]); + gimp_coords_mix (0.5, &sub[3], + 0.5, &bezier[2], + &sub[2]); + + transform_bezier_coords (matrix, sub, t_beziers, depth - 1); + + return; + } + } + } + + /* transform the segment by transforming each of the individual points. note + * that, for non-affine transforms, this is only an approximation of the real + * transformed curve, but due to subdivision it should be good enough. + */ + t_bezier = g_new (GimpCoords, 4); + + /* note that while the segments themselves are clipped to the near plane, + * their control points may still get transformed behind the camera. we + * therefore clip the control points to the near plane as well, which is not + * too meaningful, but avoids erroneously transforming them behind the + * camera. + */ + gimp_transform_polygon_coords (matrix, bezier, 2, FALSE, + t_bezier, &n); + gimp_transform_polygon_coords (matrix, bezier + 2, 2, FALSE, + t_bezier + 2, &n); + + g_queue_push_tail (t_beziers, t_bezier); +} + +/* transforms the cubic bezier segment, defined by the four control points + * 'bezier', by 'matrix', performing clipping by the near plane and subdividing + * as necessary. + * + * returns the transformed set of bezier-segment sequences in 't_beziers', as + * GQueues of GimpCoords[4] bezier-segments, and the number of sequences in + * 'n_t_beziers'. the minimal possible number of transformed sequences is 0, + * which happens when the entire segment is clipped. the maximal possible + * number of transformed sequences is 2. each sequence has at least one + * segment. + * + * if the first transformed segment is an initial segment of 'bezier', sets + * '*start_in' to TRUE, otherwise to FALSE. if the last transformed segment is + * a final segment of 'bezier', sets '*end_in' to TRUE, otherwise to FALSE. + */ +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) +{ + GimpCoords c_bezier[2][4]; + gint i; + + g_return_if_fail (matrix != NULL); + g_return_if_fail (bezier != NULL); + g_return_if_fail (t_beziers != NULL); + g_return_if_fail (n_t_beziers != NULL); + g_return_if_fail (start_in != NULL); + g_return_if_fail (end_in != NULL); + + /* if the matrix is affine, transform the easy way */ + if (gimp_matrix3_is_affine (matrix)) + { + GimpCoords *t_bezier; + + t_beziers[0] = g_queue_new (); + *n_t_beziers = 1; + + t_bezier = g_new (GimpCoords, 1); + g_queue_push_tail (t_beziers[0], t_bezier); + + for (i = 0; i < 4; i++) + { + t_bezier[i] = bezier[i]; + + gimp_matrix3_transform_point (matrix, + bezier[i].x, bezier[i].y, + &t_bezier[i].x, &t_bezier[i].y); + } + + return; + } + + /* clip the segment to the near plane */ + clip_bezier (bezier, + matrix->coeff[2][0], + matrix->coeff[2][1], + matrix->coeff[2][2] - GIMP_TRANSFORM_NEAR_Z, + c_bezier, n_t_beziers, + start_in, end_in); + + /* transform each of the resulting segments */ + for (i = 0; i < *n_t_beziers; i++) + { + t_beziers[i] = g_queue_new (); + + transform_bezier_coords (matrix, c_bezier[i], t_beziers[i], 3); + } +} |