Enable projective transforms

This is a rebase and cleanup of a patch proposed by Maarten Bosmans in
2010 to add projective transformations to Cairo.  In dusting it off,
I've removed some whitespace changes and reordered a few things to make
the patch clearer; also, the changes to cairo_rel_line_to() needed to be
moved to _cairo_default_context_rel_line_to() when cairo_backend_t was
introduced (commit 83bfd85a).

His original email follows:

"""
This is a first attempt to add projective transformations to Cairo. It
is far from complete, but mostly meant to get the discussion going
about how such a feature should be implemented.

The patch does the following:
 - add px and py, the first two elements of the bottom row to
cairo_matrix_t. These where until now always assumed zero.
 - modify the function in cairo-matrix.c to respect these elements in
all matrix operations and checks.
 - change cairo_rel_line_to such that it does the correct thing.

So is adding px and py the right approach? (btw. may be there are some
more appropriate names for these elements)
As now eight elements of the 3x3 matrix are explicitly defined and
only the last one assumed 1, it may be better to switch all the matrix
code to use the pixman floating point matrix.

The main problem with adding projective transformations is that
translational invariance isn't guaranteed anymore. So that means that
probably all uses of cairo_matrix_transform_distance must be changed.
As an example I changed cairo_rel_line_to to add the distances to the
current point itself.

If it is decided that this patch uses the right approach, at least
these things have to be addressed for the patch to be complete:
 - Check the test suite, right now there will undoubtedly be failures.
 - Add some tests to check the projective transforms.
 - Adapt other rel_ function not to use transform_distance.
 - Check how the performance is affected, mainly for the
non-projective case. Although there are some divisions added by the
patch, I suspect that the matrix operations aren't that performance
critical.

Maarten
"""

http://lists.cairographics.org/archives/cairo/2010-August/020507.html

Signed-off-by: Bryce Harrington <b.harrington@samsung.com>
Signed-off-by: Bryce Harrington <bryce@osg.samsung.com>

4 files changed, 102 insertions, 57 deletions

diff --git a/src/cairo-default-context.c b/src/cairo-default-context.c
index 1e5067bf1..a75151990 100644
--- a/src/cairo-default-context.c
+++ b/src/cairo-default-context.c
@@ -805,14 +805,18 @@ static cairo_status_t
 _cairo_default_context_rel_line_to (void *abstract_cr, double dx, double dy)
 {
     cairo_default_context_t *cr = abstract_cr;
-    cairo_fixed_t dx_fixed, dy_fixed;
+    cairo_fixed_t x_fixed, y_fixed;
+    double x, y;
 
-    _cairo_gstate_user_to_backend_distance (cr->gstate, &dx, &dy);
+    cairo_get_current_point (cr, &x, &y);
+    x += dx;
+    y += dy;
+    _cairo_gstate_user_to_backend (cr->gstate, &x, &y);
 
-    dx_fixed = _cairo_fixed_from_double (dx);
-    dy_fixed = _cairo_fixed_from_double (dy);
+    x_fixed = _cairo_fixed_from_double (x);
+    y_fixed = _cairo_fixed_from_double (y);
 
-    return _cairo_path_fixed_rel_line_to (cr->path, dx_fixed, dy_fixed);
+    return _cairo_path_fixed_line_to (cr->path, x_fixed, y_fixed);
 }
 
 
diff --git a/src/cairo-matrix.c b/src/cairo-matrix.c
index ae498f515..2f0e439eb 100644
--- a/src/cairo-matrix.c
+++ b/src/cairo-matrix.c
@@ -72,9 +72,6 @@
 static void
 _cairo_matrix_scalar_multiply (cairo_matrix_t *matrix, double scalar);
 
-static void
-_cairo_matrix_compute_adjoint (cairo_matrix_t *matrix);
-
 /**
  * cairo_matrix_init_identity:
  * @matrix: a #cairo_matrix_t
@@ -116,13 +113,13 @@ slim_hidden_def(cairo_matrix_init_identity);
 void
 cairo_matrix_init (cairo_matrix_t *matrix,
 		   double xx, double yx,
-
 		   double xy, double yy,
 		   double x0, double y0)
 {
     matrix->xx = xx; matrix->yx = yx;
     matrix->xy = xy; matrix->yy = yy;
     matrix->x0 = x0; matrix->y0 = y0;
+    matrix->px = 0; matrix->py = 0;
 }
 slim_hidden_def(cairo_matrix_init);
 
@@ -337,15 +334,24 @@ void
 cairo_matrix_multiply (cairo_matrix_t *result, const cairo_matrix_t *a, const cairo_matrix_t *b)
 {
     cairo_matrix_t r;
+    double p0;
+
+    r.xx = a->xx * b->xx + a->yx * b->xy + a->px * b->x0;
+    r.yx = a->xx * b->yx + a->yx * b->yy + a->px * b->y0;
 
-    r.xx = a->xx * b->xx + a->yx * b->xy;
-    r.yx = a->xx * b->yx + a->yx * b->yy;
+    r.xy = a->xy * b->xx + a->yy * b->xy + a->py * b->x0;
+    r.yy = a->xy * b->yx + a->yy * b->yy + a->py * b->y0;
 
-    r.xy = a->xy * b->xx + a->yy * b->xy;
-    r.yy = a->xy * b->yx + a->yy * b->yy;
+    r.x0 = a->x0 * b->xx + a->y0 * b->xy +     1 * b->x0;
+    r.y0 = a->x0 * b->yx + a->y0 * b->yy +     1 * b->y0;
 
-    r.x0 = a->x0 * b->xx + a->y0 * b->xy + b->x0;
-    r.y0 = a->x0 * b->yx + a->y0 * b->yy + b->y0;
+    r.px = a->xx * b->px + a->yx * b->py + a->px *     1;
+    r.py = a->xy * b->px + a->yy * b->py + a->py *     1;
+    p0   = a->x0 * b->px + a->y0 * b->py +     1 *     1;
+
+    printf("cairo_matrix_multiply: xx=%f yx=%f xy=%f yy=%f x0=%f y0=%f\n", r.xx, r.yx, r.xy, r.yy, r.x0, r.y0);
+
+    _cairo_matrix_scalar_multiply (&r, 1 / p0);
 
     *result = r;
 }
@@ -356,14 +362,22 @@ _cairo_matrix_multiply (cairo_matrix_t *r,
 			const cairo_matrix_t *a,
 			const cairo_matrix_t *b)
 {
-    r->xx = a->xx * b->xx + a->yx * b->xy;
-    r->yx = a->xx * b->yx + a->yx * b->yy;
+    double p0;
+
+    r->xx = a->xx * b->xx + a->yx * b->xy + a->px * b->x0;
+    r->yx = a->xx * b->yx + a->yx * b->yy + a->px * b->y0;
+
+    r->xy = a->xy * b->xx + a->yy * b->xy + a->py * b->x0;
+    r->yy = a->xy * b->yx + a->yy * b->yy + a->py * b->y0;
 
-    r->xy = a->xy * b->xx + a->yy * b->xy;
-    r->yy = a->xy * b->yx + a->yy * b->yy;
+    r->x0 = a->x0 * b->xx + a->y0 * b->xy +     1 * b->x0;
+    r->y0 = a->x0 * b->yx + a->y0 * b->yy +     1 * b->y0;
 
-    r->x0 = a->x0 * b->xx + a->y0 * b->xy + b->x0;
-    r->y0 = a->x0 * b->yx + a->y0 * b->yy + b->y0;
+    r->px = a->xx * b->px + a->yx * b->py + a->px *     1;
+    r->py = a->xy * b->px + a->yy * b->py + a->py *     1;
+    p0    = a->x0 * b->px + a->y0 * b->py +     1 *     1;
+
+    _cairo_matrix_scalar_multiply (r, 1 / p0);
 }
 
 /**
@@ -392,10 +406,13 @@ _cairo_matrix_multiply (cairo_matrix_t *r,
 void
 cairo_matrix_transform_distance (const cairo_matrix_t *matrix, double *dx, double *dy)
 {
-    double new_x, new_y;
+    double new_x, new_y, new_p;
 
     new_x = (matrix->xx * *dx + matrix->xy * *dy);
     new_y = (matrix->yx * *dx + matrix->yy * *dy);
+    new_p = matrix->px * *dx + matrix->py * *dy + 1;
+    new_x = new_x / new_p;
+    new_y = new_y / new_p;
 
     *dx = new_x;
     *dy = new_y;
@@ -415,10 +432,11 @@ slim_hidden_def(cairo_matrix_transform_distance);
 void
 cairo_matrix_transform_point (const cairo_matrix_t *matrix, double *x, double *y)
 {
-    cairo_matrix_transform_distance (matrix, x, y);
+    double new_p;
 
-    *x += matrix->x0;
-    *y += matrix->y0;
+    new_p = matrix->px * *x + matrix->py * *y + 1;
+    *x = (matrix->xx * *x + matrix->xy * *y + matrix->x0) / new_p;
+    *y = (matrix->yx * *x + matrix->yy * *y + matrix->y0) / new_p;
 }
 slim_hidden_def(cairo_matrix_transform_point);
 
@@ -550,27 +568,9 @@ _cairo_matrix_scalar_multiply (cairo_matrix_t *matrix, double scalar)
 
     matrix->x0 *= scalar;
     matrix->y0 *= scalar;
-}
-
-/* This function isn't a correct adjoint in that the implicit 1 in the
-   homogeneous result should actually be ad-bc instead. But, since this
-   adjoint is only used in the computation of the inverse, which
-   divides by det (A)=ad-bc anyway, everything works out in the end. */
-static void
-_cairo_matrix_compute_adjoint (cairo_matrix_t *matrix)
-{
-    /* adj (A) = transpose (C:cofactor (A,i,j)) */
-    double a, b, c, d, tx, ty;
-
-    _cairo_matrix_get_affine (matrix,
-			      &a,  &b,
-			      &c,  &d,
-			      &tx, &ty);
 
-    cairo_matrix_init (matrix,
-		       d, -b,
-		       -c, a,
-		       c*ty - d*tx, b*tx - a*ty);
+    matrix->px *= scalar;
+    matrix->py *= scalar;
 }
 
 /**
@@ -592,9 +592,12 @@ cairo_status_t
 cairo_matrix_invert (cairo_matrix_t *matrix)
 {
     double det;
+    double xx, xy, yx, yy;
+    double x0, y0, px, py, p0;
 
     /* Simple scaling|translation matrices are quite common... */
-    if (matrix->xy == 0. && matrix->yx == 0.) {
+    if (matrix->xy == 0. && matrix->yx == 0. &&
+	matrix->px == 0. && matrix->py == 0. ) {
 	matrix->x0 = -matrix->x0;
 	matrix->y0 = -matrix->y0;
 
@@ -626,8 +629,27 @@ cairo_matrix_invert (cairo_matrix_t *matrix)
     if (det == 0)
 	return _cairo_error (CAIRO_STATUS_INVALID_MATRIX);
 
-    _cairo_matrix_compute_adjoint (matrix);
-    _cairo_matrix_scalar_multiply (matrix, 1 / det);
+    /* Invert the matrix */
+    _cairo_matrix_get_affine (matrix,
+			      &xx, &yx,
+			      &xy, &yy,
+			      &x0, &y0);
+    px = matrix->px;
+    py = matrix->py;
+
+    cairo_matrix_init (matrix,
+		       yy    - y0*py ,  y0*px - yx    ,
+		       x0*py - xy    ,  xx    - x0*px ,
+		       xy*y0 - x0*yy ,  x0*yx - xx*y0 );
+    matrix->px = yx*py - yy*px;
+    matrix->py = xy*px - xx*py;
+
+    /* The matrix does not have to be scaled with the determinant,
+     * because it is normalized such that the bottom right element
+     * is equal to 1.
+     */
+    p0 = xx*yy - xy*yx;
+    _cairo_matrix_scalar_multiply (matrix, 1 / p0);
 
     return CAIRO_STATUS_SUCCESS;
 }
@@ -656,11 +678,21 @@ double
 _cairo_matrix_compute_determinant (const cairo_matrix_t *matrix)
 {
     double a, b, c, d;
+    double det;
 
-    a = matrix->xx; b = matrix->yx;
-    c = matrix->xy; d = matrix->yy;
+    a = matrix->xx; b = matrix->xy;
+    c = matrix->yx; d = matrix->yy;
+    det = a*d - b*c;
 
-    return a*d - b*c;
+    a = matrix->xx; b = matrix->x0;
+    c = matrix->yx; d = matrix->y0;
+    det -= (a*d - b*c) * matrix->py;
+
+    a = matrix->xy; b = matrix->x0;
+    c = matrix->yy; d = matrix->y0;
+    det += (a*d - b*c) * matrix->px;
+
+    return det;
 }
 
 /**
@@ -759,8 +791,14 @@ _cairo_matrix_is_integer_translation (const cairo_matrix_t *matrix,
 cairo_bool_t
 _cairo_matrix_has_unity_scale (const cairo_matrix_t *matrix)
 {
+    double det;
+
+    if (! (matrix->px == 0.0 && matrix->py == 0.0))
+       return FALSE;
+
     /* check that the determinant is near +/-1 */
-    double det = _cairo_matrix_compute_determinant (matrix);
+    det = _cairo_matrix_compute_determinant (matrix);
+
     if (fabs (det * det - 1.0) < SCALING_EPSILON) {
 	/* check that one axis is close to zero */
 	if (fabs (matrix->xy) < SCALING_EPSILON  &&
@@ -968,8 +1006,8 @@ _cairo_matrix_to_pixman_matrix (const cairo_matrix_t	*matrix,
     pixman_transform->matrix[1][1] = _cairo_fixed_16_16_from_double (matrix->yy);
     pixman_transform->matrix[1][2] = _cairo_fixed_16_16_from_double (matrix->y0);
 
-    pixman_transform->matrix[2][0] = 0;
-    pixman_transform->matrix[2][1] = 0;
+    pixman_transform->matrix[2][0] = _cairo_fixed_16_16_from_double (matrix->px);
+    pixman_transform->matrix[2][1] = _cairo_fixed_16_16_from_double (matrix->py);
     pixman_transform->matrix[2][2] = 1 << 16;
 
     /* The conversion above breaks cairo's translation invariance:
@@ -1013,7 +1051,7 @@ _cairo_matrix_to_pixman_matrix (const cairo_matrix_t	*matrix,
 	vector.vector[2] = 1 << 16;
 
 	/* If we can't transform the reference point, skip the adjustment. */
-	if (! pixman_transform_point_3d (pixman_transform, &vector))
+	if (! pixman_transform_point (pixman_transform, &vector))
 	    return CAIRO_STATUS_SUCCESS;
 
 	x = pixman_fixed_to_double (vector.vector[0]);
diff --git a/src/cairo.h b/src/cairo.h
index 3104d47e4..7e73e1102 100644
--- a/src/cairo.h
+++ b/src/cairo.h
@@ -193,6 +193,7 @@ typedef struct _cairo_matrix {
     double xx; double yx;
     double xy; double yy;
     double x0; double y0;
+    double px; double py;
 } cairo_matrix_t;
 
 /**
diff --git a/src/cairoint.h b/src/cairoint.h
index 5bca00389..026278765 100644
--- a/src/cairoint.h
+++ b/src/cairoint.h
@@ -1733,14 +1733,16 @@ _cairo_matrix_is_identity (const cairo_matrix_t *matrix)
 {
     return (matrix->xx == 1.0 && matrix->yx == 0.0 &&
 	    matrix->xy == 0.0 && matrix->yy == 1.0 &&
-	    matrix->x0 == 0.0 && matrix->y0 == 0.0);
+	    matrix->x0 == 0.0 && matrix->y0 == 0.0 &&
+	    matrix->px == 0.0 && matrix->py == 0.0);
 }
 
 static inline cairo_bool_t
 _cairo_matrix_is_translation (const cairo_matrix_t *matrix)
 {
     return (matrix->xx == 1.0 && matrix->yx == 0.0 &&
-	    matrix->xy == 0.0 && matrix->yy == 1.0);
+	    matrix->xy == 0.0 && matrix->yy == 1.0 &&
+	    matrix->px == 0.0 && matrix->py == 0.0);
 }
 
 static inline cairo_bool_t