src/cairo-matrix.c

   1 /* cairo - a vector graphics library with display and print output
   2  *
   3  * Copyright © 2002 University of Southern California
   4  *
   5  * This library is free software; you can redistribute it and/or
   6  * modify it either under the terms of the GNU Lesser General Public
   7  * License version 2.1 as published by the Free Software Foundation
   8  * (the "LGPL") or, at your option, under the terms of the Mozilla
   9  * Public License Version 1.1 (the "MPL"). If you do not alter this
  10  * notice, a recipient may use your version of this file under either
  11  * the MPL or the LGPL.
  12  *
  13  * You should have received a copy of the LGPL along with this library
  14  * in the file COPYING-LGPL-2.1; if not, write to the Free Software
  15  * Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA
  16  * You should have received a copy of the MPL along with this library
  17  * in the file COPYING-MPL-1.1
  18  *
  19  * The contents of this file are subject to the Mozilla Public License
  20  * Version 1.1 (the "License"); you may not use this file except in
  21  * compliance with the License. You may obtain a copy of the License at
  22  * http://www.mozilla.org/MPL/
  23  *
  24  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
  25  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
  26  * the specific language governing rights and limitations.
  27  *
  28  * The Original Code is the cairo graphics library.
  29  *
  30  * The Initial Developer of the Original Code is University of Southern
  31  * California.
  32  *
  33  * Contributor(s):
  34  *      Carl D. Worth <cworth@cworth.org>
  35  */
  36
  37 #include "cairoint.h"
  38 #include "cairo-error-private.h"
  39 #include <float.h>
  40
  41 #define PIXMAN_MAX_INT ((pixman_fixed_1 >> 1) - pixman_fixed_e) /* need to ensure deltas also fit */
  42
  43 #if _XOPEN_SOURCE >= 600 || defined (_ISOC99_SOURCE)
  44 #define ISFINITE(x) isfinite (x)
  45 #else
  46 #define ISFINITE(x) ((x) * (x) >= 0.) /* check for NaNs */
  47 #endif
  48
  49 /**
  50  * SECTION:cairo-matrix
  51  * @Title: cairo_matrix_t
  52  * @Short_Description: Generic matrix operations
  53  * @See_Also: #cairo_t
  54  *
  55  * #cairo_matrix_t is used throughout cairo to convert between different
  56  * coordinate spaces.  A #cairo_matrix_t holds an affine transformation,
  57  * such as a scale, rotation, shear, or a combination of these.
  58  * The transformation of a point (<literal>x</literal>,<literal>y</literal>)
  59  * is given by:
  60  *
  61  * <programlisting>
  62  * x_new = xx * x + xy * y + x0;
  63  * y_new = yx * x + yy * y + y0;
  64  * </programlisting>
  65  *
  66  * The current transformation matrix of a #cairo_t, represented as a
  67  * #cairo_matrix_t, defines the transformation from user-space
  68  * coordinates to device-space coordinates. See cairo_get_matrix() and
  69  * cairo_set_matrix().
  70  **/
  71
  72 static void
  73 _cairo_matrix_scalar_multiply (cairo_matrix_t *matrix, double scalar);
  74
  75 static void
  76 _cairo_matrix_compute_adjoint (cairo_matrix_t *matrix);
  77
  78 /**
  79  * cairo_matrix_init_identity:
  80  * @matrix: a #cairo_matrix_t
  81  *
  82  * Modifies @matrix to be an identity transformation.
  83  *
  84  * Since: 1.0
  85  **/
  86 void
  87 cairo_matrix_init_identity (cairo_matrix_t *matrix)
  88 {
  89     cairo_matrix_init (matrix,
  90                        1, 0,
  91                        0, 1,
  92                        0, 0);
  93 }
  94 slim_hidden_def(cairo_matrix_init_identity);
  95
  96 /**
  97  * cairo_matrix_init:
  98  * @matrix: a #cairo_matrix_t
  99  * @xx: xx component of the affine transformation
 100  * @yx: yx component of the affine transformation
 101  * @xy: xy component of the affine transformation
 102  * @yy: yy component of the affine transformation
 103  * @x0: X translation component of the affine transformation
 104  * @y0: Y translation component of the affine transformation
 105  *
 106  * Sets @matrix to be the affine transformation given by
 107  * @xx, @yx, @xy, @yy, @x0, @y0. The transformation is given
 108  * by:
 109  * <programlisting>
 110  *  x_new = xx * x + xy * y + x0;
 111  *  y_new = yx * x + yy * y + y0;
 112  * </programlisting>
 113  *
 114  * Since: 1.0
 115  **/
 116 void
 117 cairo_matrix_init (cairo_matrix_t *matrix,
 118                    double xx, double yx,
 119
 120                    double xy, double yy,
 121                    double x0, double y0)
 122 {
 123     matrix->xx = xx; matrix->yx = yx;
 124     matrix->xy = xy; matrix->yy = yy;
 125     matrix->x0 = x0; matrix->y0 = y0;
 126 }
 127 slim_hidden_def(cairo_matrix_init);
 128
 129 /**
 130  * _cairo_matrix_get_affine:
 131  * @matrix: a #cairo_matrix_t
 132  * @xx: location to store xx component of matrix
 133  * @yx: location to store yx component of matrix
 134  * @xy: location to store xy component of matrix
 135  * @yy: location to store yy component of matrix
 136  * @x0: location to store x0 (X-translation component) of matrix, or %NULL
 137  * @y0: location to store y0 (Y-translation component) of matrix, or %NULL
 138  *
 139  * Gets the matrix values for the affine transformation that @matrix represents.
 140  * See cairo_matrix_init().
 141  *
 142  *
 143  * This function is a leftover from the old public API, but is still
 144  * mildly useful as an internal means for getting at the matrix
 145  * members in a positional way. For example, when reassigning to some
 146  * external matrix type, or when renaming members to more meaningful
 147  * names (such as a,b,c,d,e,f) for particular manipulations.
 148  **/
 149 void
 150 _cairo_matrix_get_affine (const cairo_matrix_t *matrix,
 151                           double *xx, double *yx,
 152                           double *xy, double *yy,
 153                           double *x0, double *y0)
 154 {
 155     *xx  = matrix->xx;
 156     *yx  = matrix->yx;
 157
 158     *xy  = matrix->xy;
 159     *yy  = matrix->yy;
 160
 161     if (x0)
 162         *x0 = matrix->x0;
 163     if (y0)
 164         *y0 = matrix->y0;
 165 }
 166
 167 /**
 168  * cairo_matrix_init_translate:
 169  * @matrix: a #cairo_matrix_t
 170  * @tx: amount to translate in the X direction
 171  * @ty: amount to translate in the Y direction
 172  *
 173  * Initializes @matrix to a transformation that translates by @tx and
 174  * @ty in the X and Y dimensions, respectively.
 175  *
 176  * Since: 1.0
 177  **/
 178 void
 179 cairo_matrix_init_translate (cairo_matrix_t *matrix,
 180                              double tx, double ty)
 181 {
 182     cairo_matrix_init (matrix,
 183                        1, 0,
 184                        0, 1,
 185                        tx, ty);
 186 }
 187 slim_hidden_def(cairo_matrix_init_translate);
 188
 189 /**
 190  * cairo_matrix_translate:
 191  * @matrix: a #cairo_matrix_t
 192  * @tx: amount to translate in the X direction
 193  * @ty: amount to translate in the Y direction
 194  *
 195  * Applies a translation by @tx, @ty to the transformation in
 196  * @matrix. The effect of the new transformation is to first translate
 197  * the coordinates by @tx and @ty, then apply the original transformation
 198  * to the coordinates.
 199  *
 200  * Since: 1.0
 201  **/
 202 void
 203 cairo_matrix_translate (cairo_matrix_t *matrix, double tx, double ty)
 204 {
 205     cairo_matrix_t tmp;
 206
 207     cairo_matrix_init_translate (&tmp, tx, ty);
 208
 209     cairo_matrix_multiply (matrix, &tmp, matrix);
 210 }
 211 slim_hidden_def (cairo_matrix_translate);
 212
 213 /**
 214  * cairo_matrix_init_scale:
 215  * @matrix: a #cairo_matrix_t
 216  * @sx: scale factor in the X direction
 217  * @sy: scale factor in the Y direction
 218  *
 219  * Initializes @matrix to a transformation that scales by @sx and @sy
 220  * in the X and Y dimensions, respectively.
 221  *
 222  * Since: 1.0
 223  **/
 224 void
 225 cairo_matrix_init_scale (cairo_matrix_t *matrix,
 226                          double sx, double sy)
 227 {
 228     cairo_matrix_init (matrix,
 229                        sx,  0,
 230                        0, sy,
 231                        0, 0);
 232 }
 233 slim_hidden_def(cairo_matrix_init_scale);
 234
 235 /**
 236  * cairo_matrix_scale:
 237  * @matrix: a #cairo_matrix_t
 238  * @sx: scale factor in the X direction
 239  * @sy: scale factor in the Y direction
 240  *
 241  * Applies scaling by @sx, @sy to the transformation in @matrix. The
 242  * effect of the new transformation is to first scale the coordinates
 243  * by @sx and @sy, then apply the original transformation to the coordinates.
 244  *
 245  * Since: 1.0
 246  **/
 247 void
 248 cairo_matrix_scale (cairo_matrix_t *matrix, double sx, double sy)
 249 {
 250     cairo_matrix_t tmp;
 251
 252     cairo_matrix_init_scale (&tmp, sx, sy);
 253
 254     cairo_matrix_multiply (matrix, &tmp, matrix);
 255 }
 256 slim_hidden_def(cairo_matrix_scale);
 257
 258 /**
 259  * cairo_matrix_init_rotate:
 260  * @matrix: a #cairo_matrix_t
 261  * @radians: angle of rotation, in radians. The direction of rotation
 262  * is defined such that positive angles rotate in the direction from
 263  * the positive X axis toward the positive Y axis. With the default
 264  * axis orientation of cairo, positive angles rotate in a clockwise
 265  * direction.
 266  *
 267  * Initialized @matrix to a transformation that rotates by @radians.
 268  *
 269  * Since: 1.0
 270  **/
 271 void
 272 cairo_matrix_init_rotate (cairo_matrix_t *matrix,
 273                           double radians)
 274 {
 275     double  s;
 276     double  c;
 277
 278     s = sin (radians);
 279     c = cos (radians);
 280
 281     cairo_matrix_init (matrix,
 282                        c, s,
 283                        -s, c,
 284                        0, 0);
 285 }
 286 slim_hidden_def(cairo_matrix_init_rotate);
 287
 288 /**
 289  * cairo_matrix_rotate:
 290  * @matrix: a #cairo_matrix_t
 291  * @radians: angle of rotation, in radians. The direction of rotation
 292  * is defined such that positive angles rotate in the direction from
 293  * the positive X axis toward the positive Y axis. With the default
 294  * axis orientation of cairo, positive angles rotate in a clockwise
 295  * direction.
 296  *
 297  * Applies rotation by @radians to the transformation in
 298  * @matrix. The effect of the new transformation is to first rotate the
 299  * coordinates by @radians, then apply the original transformation
 300  * to the coordinates.
 301  *
 302  * Since: 1.0
 303  **/
 304 void
 305 cairo_matrix_rotate (cairo_matrix_t *matrix, double radians)
 306 {
 307     cairo_matrix_t tmp;
 308
 309     cairo_matrix_init_rotate (&tmp, radians);
 310
 311     cairo_matrix_multiply (matrix, &tmp, matrix);
 312 }
 313
 314 /**
 315  * cairo_matrix_multiply:
 316  * @result: a #cairo_matrix_t in which to store the result
 317  * @a: a #cairo_matrix_t
 318  * @b: a #cairo_matrix_t
 319  *
 320  * Multiplies the affine transformations in @a and @b together
 321  * and stores the result in @result. The effect of the resulting
 322  * transformation is to first apply the transformation in @a to the
 323  * coordinates and then apply the transformation in @b to the
 324  * coordinates.
 325  *
 326  * It is allowable for @result to be identical to either @a or @b.
 327  *
 328  * Since: 1.0
 329  **/
 330 /*
 331  * XXX: The ordering of the arguments to this function corresponds
 332  *      to [row_vector]*A*B. If we want to use column vectors instead,
 333  *      then we need to switch the two arguments and fix up all
 334  *      uses.
 335  */
 336 void
 337 cairo_matrix_multiply (cairo_matrix_t *result, const cairo_matrix_t *a, const cairo_matrix_t *b)
 338 {
 339     cairo_matrix_t r;
 340
 341     r.xx = a->xx * b->xx + a->yx * b->xy;
 342     r.yx = a->xx * b->yx + a->yx * b->yy;
 343
 344     r.xy = a->xy * b->xx + a->yy * b->xy;
 345     r.yy = a->xy * b->yx + a->yy * b->yy;
 346
 347     r.x0 = a->x0 * b->xx + a->y0 * b->xy + b->x0;
 348     r.y0 = a->x0 * b->yx + a->y0 * b->yy + b->y0;
 349
 350     *result = r;
 351 }
 352 slim_hidden_def(cairo_matrix_multiply);
 353
 354 void
 355 _cairo_matrix_multiply (cairo_matrix_t *r,
 356                         const cairo_matrix_t *a,
 357                         const cairo_matrix_t *b)
 358 {
 359     r->xx = a->xx * b->xx + a->yx * b->xy;
 360     r->yx = a->xx * b->yx + a->yx * b->yy;
 361
 362     r->xy = a->xy * b->xx + a->yy * b->xy;
 363     r->yy = a->xy * b->yx + a->yy * b->yy;
 364
 365     r->x0 = a->x0 * b->xx + a->y0 * b->xy + b->x0;
 366     r->y0 = a->x0 * b->yx + a->y0 * b->yy + b->y0;
 367 }
 368
 369 /**
 370  * cairo_matrix_transform_distance:
 371  * @matrix: a #cairo_matrix_t
 372  * @dx: X component of a distance vector. An in/out parameter
 373  * @dy: Y component of a distance vector. An in/out parameter
 374  *
 375  * Transforms the distance vector (@dx,@dy) by @matrix. This is
 376  * similar to cairo_matrix_transform_point() except that the translation
 377  * components of the transformation are ignored. The calculation of
 378  * the returned vector is as follows:
 379  *
 380  * <programlisting>
 381  * dx2 = dx1 * a + dy1 * c;
 382  * dy2 = dx1 * b + dy1 * d;
 383  * </programlisting>
 384  *
 385  * Affine transformations are position invariant, so the same vector
 386  * always transforms to the same vector. If (@x1,@y1) transforms
 387  * to (@x2,@y2) then (@x1+@dx1,@y1+@dy1) will transform to
 388  * (@x1+@dx2,@y1+@dy2) for all values of @x1 and @x2.
 389  *
 390  * Since: 1.0
 391  **/
 392 void
 393 cairo_matrix_transform_distance (const cairo_matrix_t *matrix, double *dx, double *dy)
 394 {
 395     double new_x, new_y;
 396
 397     new_x = (matrix->xx * *dx + matrix->xy * *dy);
 398     new_y = (matrix->yx * *dx + matrix->yy * *dy);
 399
 400     *dx = new_x;
 401     *dy = new_y;
 402 }
 403 slim_hidden_def(cairo_matrix_transform_distance);
 404
 405 /**
 406  * cairo_matrix_transform_point:
 407  * @matrix: a #cairo_matrix_t
 408  * @x: X position. An in/out parameter
 409  * @y: Y position. An in/out parameter
 410  *
 411  * Transforms the point (@x, @y) by @matrix.
 412  *
 413  * Since: 1.0
 414  **/
 415 void
 416 cairo_matrix_transform_point (const cairo_matrix_t *matrix, double *x, double *y)
 417 {
 418     cairo_matrix_transform_distance (matrix, x, y);
 419
 420     *x += matrix->x0;
 421     *y += matrix->y0;
 422 }
 423 slim_hidden_def(cairo_matrix_transform_point);
 424
 425 void
 426 _cairo_matrix_transform_bounding_box (const cairo_matrix_t *matrix,
 427                                       double *x1, double *y1,
 428                                       double *x2, double *y2,
 429                                       cairo_bool_t *is_tight)
 430 {
 431     int i;
 432     double quad_x[4], quad_y[4];
 433     double min_x, max_x;
 434     double min_y, max_y;
 435
 436     if (matrix->xy == 0. && matrix->yx == 0.) {
 437         /* non-rotation/skew matrix, just map the two extreme points */
 438
 439         if (matrix->xx != 1.) {
 440             quad_x[0] = *x1 * matrix->xx;
 441             quad_x[1] = *x2 * matrix->xx;
 442             if (quad_x[0] < quad_x[1]) {
 443                 *x1 = quad_x[0];
 444                 *x2 = quad_x[1];
 445             } else {
 446                 *x1 = quad_x[1];
 447                 *x2 = quad_x[0];
 448             }
 449         }
 450         if (matrix->x0 != 0.) {
 451             *x1 += matrix->x0;
 452             *x2 += matrix->x0;
 453         }
 454
 455         if (matrix->yy != 1.) {
 456             quad_y[0] = *y1 * matrix->yy;
 457             quad_y[1] = *y2 * matrix->yy;
 458             if (quad_y[0] < quad_y[1]) {
 459                 *y1 = quad_y[0];
 460                 *y2 = quad_y[1];
 461             } else {
 462                 *y1 = quad_y[1];
 463                 *y2 = quad_y[0];
 464             }
 465         }
 466         if (matrix->y0 != 0.) {
 467             *y1 += matrix->y0;
 468             *y2 += matrix->y0;
 469         }
 470
 471         if (is_tight)
 472             *is_tight = TRUE;
 473
 474         return;
 475     }
 476
 477     /* general matrix */
 478     quad_x[0] = *x1;
 479     quad_y[0] = *y1;
 480     cairo_matrix_transform_point (matrix, &quad_x[0], &quad_y[0]);
 481
 482     quad_x[1] = *x2;
 483     quad_y[1] = *y1;
 484     cairo_matrix_transform_point (matrix, &quad_x[1], &quad_y[1]);
 485
 486     quad_x[2] = *x1;
 487     quad_y[2] = *y2;
 488     cairo_matrix_transform_point (matrix, &quad_x[2], &quad_y[2]);
 489
 490     quad_x[3] = *x2;
 491     quad_y[3] = *y2;
 492     cairo_matrix_transform_point (matrix, &quad_x[3], &quad_y[3]);
 493
 494     min_x = max_x = quad_x[0];
 495     min_y = max_y = quad_y[0];
 496
 497     for (i=1; i < 4; i++) {
 498         if (quad_x[i] < min_x)
 499             min_x = quad_x[i];
 500         if (quad_x[i] > max_x)
 501             max_x = quad_x[i];
 502
 503         if (quad_y[i] < min_y)
 504             min_y = quad_y[i];
 505         if (quad_y[i] > max_y)
 506             max_y = quad_y[i];
 507     }
 508
 509     *x1 = min_x;
 510     *y1 = min_y;
 511     *x2 = max_x;
 512     *y2 = max_y;
 513
 514     if (is_tight) {
 515         /* it's tight if and only if the four corner points form an axis-aligned
 516            rectangle.
 517            And that's true if and only if we can derive corners 0 and 3 from
 518            corners 1 and 2 in one of two straightforward ways...
 519            We could use a tolerance here but for now we'll fall back to FALSE in the case
 520            of floating point error.
 521         */
 522         *is_tight =
 523             (quad_x[1] == quad_x[0] && quad_y[1] == quad_y[3] &&
 524              quad_x[2] == quad_x[3] && quad_y[2] == quad_y[0]) ||
 525             (quad_x[1] == quad_x[3] && quad_y[1] == quad_y[0] &&
 526              quad_x[2] == quad_x[0] && quad_y[2] == quad_y[3]);
 527     }
 528 }
 529
 530 cairo_private void
 531 _cairo_matrix_transform_bounding_box_fixed (const cairo_matrix_t *matrix,
 532                                             cairo_box_t          *bbox,
 533                                             cairo_bool_t *is_tight)
 534 {
 535     double x1, y1, x2, y2;
 536
 537     _cairo_box_to_doubles (bbox, &x1, &y1, &x2, &y2);
 538     _cairo_matrix_transform_bounding_box (matrix, &x1, &y1, &x2, &y2, is_tight);
 539     _cairo_box_from_doubles (bbox, &x1, &y1, &x2, &y2);
 540 }
 541
 542 static void
 543 _cairo_matrix_scalar_multiply (cairo_matrix_t *matrix, double scalar)
 544 {
 545     matrix->xx *= scalar;
 546     matrix->yx *= scalar;
 547
 548     matrix->xy *= scalar;
 549     matrix->yy *= scalar;
 550
 551     matrix->x0 *= scalar;
 552     matrix->y0 *= scalar;
 553 }
 554
 555 /* This function isn't a correct adjoint in that the implicit 1 in the
 556    homogeneous result should actually be ad-bc instead. But, since this
 557    adjoint is only used in the computation of the inverse, which
 558    divides by det (A)=ad-bc anyway, everything works out in the end. */
 559 static void
 560 _cairo_matrix_compute_adjoint (cairo_matrix_t *matrix)
 561 {
 562     /* adj (A) = transpose (C:cofactor (A,i,j)) */
 563     double a, b, c, d, tx, ty;
 564
 565     _cairo_matrix_get_affine (matrix,
 566                               &a,  &b,
 567                               &c,  &d,
 568                               &tx, &ty);
 569
 570     cairo_matrix_init (matrix,
 571                        d, -b,
 572                        -c, a,
 573                        c*ty - d*tx, b*tx - a*ty);
 574 }
 575
 576 /**
 577  * cairo_matrix_invert:
 578  * @matrix: a #cairo_matrix_t
 579  *
 580  * Changes @matrix to be the inverse of its original value. Not
 581  * all transformation matrices have inverses; if the matrix
 582  * collapses points together (it is <firstterm>degenerate</firstterm>),
 583  * then it has no inverse and this function will fail.
 584  *
 585  * Returns: If @matrix has an inverse, modifies @matrix to
 586  *  be the inverse matrix and returns %CAIRO_STATUS_SUCCESS. Otherwise,
 587  *  returns %CAIRO_STATUS_INVALID_MATRIX.
 588  *
 589  * Since: 1.0
 590  **/
 591 cairo_status_t
 592 cairo_matrix_invert (cairo_matrix_t *matrix)
 593 {
 594     double det;
 595
 596     /* Simple scaling|translation matrices are quite common... */
 597     if (matrix->xy == 0. && matrix->yx == 0.) {
 598         matrix->x0 = -matrix->x0;
 599         matrix->y0 = -matrix->y0;
 600
 601         if (matrix->xx != 1.) {
 602             if (matrix->xx == 0.)
 603                 return _cairo_error (CAIRO_STATUS_INVALID_MATRIX);
 604
 605             matrix->xx = 1. / matrix->xx;
 606             matrix->x0 *= matrix->xx;
 607         }
 608
 609         if (matrix->yy != 1.) {
 610             if (matrix->yy == 0.)
 611                 return _cairo_error (CAIRO_STATUS_INVALID_MATRIX);
 612
 613             matrix->yy = 1. / matrix->yy;
 614             matrix->y0 *= matrix->yy;
 615         }
 616
 617         return CAIRO_STATUS_SUCCESS;
 618     }
 619
 620     /* inv (A) = 1/det (A) * adj (A) */
 621     det = _cairo_matrix_compute_determinant (matrix);
 622
 623     if (! ISFINITE (det))
 624         return _cairo_error (CAIRO_STATUS_INVALID_MATRIX);
 625
 626     if (det == 0)
 627         return _cairo_error (CAIRO_STATUS_INVALID_MATRIX);
 628
 629     _cairo_matrix_compute_adjoint (matrix);
 630     _cairo_matrix_scalar_multiply (matrix, 1 / det);
 631
 632     return CAIRO_STATUS_SUCCESS;
 633 }
 634 slim_hidden_def(cairo_matrix_invert);
 635
 636 cairo_bool_t
 637 _cairo_matrix_is_invertible (const cairo_matrix_t *matrix)
 638 {
 639     double det;
 640
 641     det = _cairo_matrix_compute_determinant (matrix);
 642
 643     return ISFINITE (det) && det != 0.;
 644 }
 645
 646 cairo_bool_t
 647 _cairo_matrix_is_scale_0 (const cairo_matrix_t *matrix)
 648 {
 649     return matrix->xx == 0. &&
 650            matrix->xy == 0. &&
 651            matrix->yx == 0. &&
 652            matrix->yy == 0.;
 653 }
 654
 655 double
 656 _cairo_matrix_compute_determinant (const cairo_matrix_t *matrix)
 657 {
 658     double a, b, c, d;
 659
 660     a = matrix->xx; b = matrix->yx;
 661     c = matrix->xy; d = matrix->yy;
 662
 663     return a*d - b*c;
 664 }
 665
 666 /**
 667  * _cairo_matrix_compute_basis_scale_factors:
 668  * @matrix: a matrix
 669  * @basis_scale: the scale factor in the direction of basis
 670  * @normal_scale: the scale factor in the direction normal to the basis
 671  * @x_basis: basis to use.  X basis if true, Y basis otherwise.
 672  *
 673  * Computes |Mv| and det(M)/|Mv| for v=[1,0] if x_basis is true, and v=[0,1]
 674  * otherwise, and M is @matrix.
 675  *
 676  * Return value: the scale factor of @matrix on the height of the font,
 677  * or 1.0 if @matrix is %NULL.
 678  **/
 679 cairo_status_t
 680 _cairo_matrix_compute_basis_scale_factors (const cairo_matrix_t *matrix,
 681                                            double *basis_scale, double *normal_scale,
 682                                            cairo_bool_t x_basis)
 683 {
 684     double det;
 685
 686     det = _cairo_matrix_compute_determinant (matrix);
 687
 688     if (! ISFINITE (det))
 689         return _cairo_error (CAIRO_STATUS_INVALID_MATRIX);
 690
 691     if (det == 0)
 692     {
 693         *basis_scale = *normal_scale = 0;
 694     }
 695     else
 696     {
 697         double x = x_basis != 0;
 698         double y = x == 0;
 699         double major, minor;
 700
 701         cairo_matrix_transform_distance (matrix, &x, &y);
 702         major = hypot (x, y);
 703         /*
 704          * ignore mirroring
 705          */
 706         if (det < 0)
 707             det = -det;
 708         if (major)
 709             minor = det / major;
 710         else
 711             minor = 0.0;
 712         if (x_basis)
 713         {
 714             *basis_scale = major;
 715             *normal_scale = minor;
 716         }
 717         else
 718         {
 719             *basis_scale = minor;
 720             *normal_scale = major;
 721         }
 722     }
 723
 724     return CAIRO_STATUS_SUCCESS;
 725 }
 726
 727 cairo_bool_t
 728 _cairo_matrix_is_integer_translation (const cairo_matrix_t *matrix,
 729                                       int *itx, int *ity)
 730 {
 731     if (_cairo_matrix_is_translation (matrix))
 732     {
 733         cairo_fixed_t x0_fixed = _cairo_fixed_from_double (matrix->x0);
 734         cairo_fixed_t y0_fixed = _cairo_fixed_from_double (matrix->y0);
 735
 736         if (_cairo_fixed_is_integer (x0_fixed) &&
 737             _cairo_fixed_is_integer (y0_fixed))
 738         {
 739             if (itx)
 740                 *itx = _cairo_fixed_integer_part (x0_fixed);
 741             if (ity)
 742                 *ity = _cairo_fixed_integer_part (y0_fixed);
 743
 744             return TRUE;
 745         }
 746     }
 747
 748     return FALSE;
 749 }
 750
 751 cairo_bool_t
 752 _cairo_matrix_has_unity_scale (const cairo_matrix_t *matrix)
 753 {
 754     if (matrix->xy == 0.0 && matrix->yx == 0.0) {
 755         if (! (matrix->xx == 1.0 || matrix->xx == -1.0))
 756             return FALSE;
 757         if (! (matrix->yy == 1.0 || matrix->yy == -1.0))
 758             return FALSE;
 759     } else if (matrix->xx == 0.0 && matrix->yy == 0.0) {
 760         if (! (matrix->xy == 1.0 || matrix->xy == -1.0))
 761             return FALSE;
 762         if (! (matrix->yx == 1.0 || matrix->yx == -1.0))
 763             return FALSE;
 764     } else
 765         return FALSE;
 766
 767     return TRUE;
 768 }
 769
 770 /* By pixel exact here, we mean a matrix that is composed only of
 771  * 90 degree rotations, flips, and integer translations and produces a 1:1
 772  * mapping between source and destination pixels. If we transform an image
 773  * with a pixel-exact matrix, filtering is not useful.
 774  */
 775 cairo_bool_t
 776 _cairo_matrix_is_pixel_exact (const cairo_matrix_t *matrix)
 777 {
 778     cairo_fixed_t x0_fixed, y0_fixed;
 779
 780     if (! _cairo_matrix_has_unity_scale (matrix))
 781         return FALSE;
 782
 783     x0_fixed = _cairo_fixed_from_double (matrix->x0);
 784     y0_fixed = _cairo_fixed_from_double (matrix->y0);
 785
 786     return _cairo_fixed_is_integer (x0_fixed) && _cairo_fixed_is_integer (y0_fixed);
 787 }
 788
 789 /*
 790   A circle in user space is transformed into an ellipse in device space.
 791
 792   The following is a derivation of a formula to calculate the length of the
 793   major axis for this ellipse; this is useful for error bounds calculations.
 794
 795   Thanks to Walter Brisken <wbrisken@aoc.nrao.edu> for this derivation:
 796
 797   1.  First some notation:
 798
 799   All capital letters represent vectors in two dimensions.  A prime '
 800   represents a transformed coordinate.  Matrices are written in underlined
 801   form, ie _R_.  Lowercase letters represent scalar real values.
 802
 803   2.  The question has been posed:  What is the maximum expansion factor
 804   achieved by the linear transformation
 805
 806   X' = X _R_
 807
 808   where _R_ is a real-valued 2x2 matrix with entries:
 809
 810   _R_ = [a b]
 811         [c d]  .
 812
 813   In other words, what is the maximum radius, MAX[ |X'| ], reached for any
 814   X on the unit circle ( |X| = 1 ) ?
 815
 816   3.  Some useful formulae
 817
 818   (A) through (C) below are standard double-angle formulae.  (D) is a lesser
 819   known result and is derived below:
 820
 821   (A)  sin²(θ) = (1 - cos(2*θ))/2
 822   (B)  cos²(θ) = (1 + cos(2*θ))/2
 823   (C)  sin(θ)*cos(θ) = sin(2*θ)/2
 824   (D)  MAX[a*cos(θ) + b*sin(θ)] = sqrt(a² + b²)
 825
 826   Proof of (D):
 827
 828   find the maximum of the function by setting the derivative to zero:
 829
 830        -a*sin(θ)+b*cos(θ) = 0
 831
 832   From this it follows that
 833
 834        tan(θ) = b/a
 835
 836   and hence
 837
 838        sin(θ) = b/sqrt(a² + b²)
 839
 840   and
 841
 842        cos(θ) = a/sqrt(a² + b²)
 843
 844   Thus the maximum value is
 845
 846        MAX[a*cos(θ) + b*sin(θ)] = (a² + b²)/sqrt(a² + b²)
 847                                    = sqrt(a² + b²)
 848
 849   4.  Derivation of maximum expansion
 850
 851   To find MAX[ |X'| ] we search brute force method using calculus.  The unit
 852   circle on which X is constrained is to be parameterized by t:
 853
 854        X(θ) = (cos(θ), sin(θ))
 855
 856   Thus
 857
 858        X'(θ) = X(θ) * _R_ = (cos(θ), sin(θ)) * [a b]
 859                                                [c d]
 860              = (a*cos(θ) + c*sin(θ), b*cos(θ) + d*sin(θ)).
 861
 862   Define
 863
 864        r(θ) = |X'(θ)|
 865
 866   Thus
 867
 868        r²(θ) = (a*cos(θ) + c*sin(θ))² + (b*cos(θ) + d*sin(θ))²
 869              = (a² + b²)*cos²(θ) + (c² + d²)*sin²(θ)
 870                  + 2*(a*c + b*d)*cos(θ)*sin(θ)
 871
 872   Now apply the double angle formulae (A) to (C) from above:
 873
 874        r²(θ) = (a² + b² + c² + d²)/2
 875              + (a² + b² - c² - d²)*cos(2*θ)/2
 876              + (a*c + b*d)*sin(2*θ)
 877              = f + g*cos(φ) + h*sin(φ)
 878
 879   Where
 880
 881        f = (a² + b² + c² + d²)/2
 882        g = (a² + b² - c² - d²)/2
 883        h = (a*c + d*d)
 884        φ = 2*θ
 885
 886   It is clear that MAX[ |X'| ] = sqrt(MAX[ r² ]).  Here we determine MAX[ r² ]
 887   using (D) from above:
 888
 889        MAX[ r² ] = f + sqrt(g² + h²)
 890
 891   And finally
 892
 893        MAX[ |X'| ] = sqrt( f + sqrt(g² + h²) )
 894
 895   Which is the solution to this problem.
 896
 897   Walter Brisken
 898   2004/10/08
 899
 900   (Note that the minor axis length is at the minimum of the above solution,
 901   which is just sqrt ( f - sqrt(g² + h²) ) given the symmetry of (D)).
 902
 903
 904   For another derivation of the same result, using Singular Value Decomposition,
 905   see doc/tutorial/src/singular.c.
 906 */
 907
 908 /* determine the length of the major axis of a circle of the given radius
 909    after applying the transformation matrix. */
 910 double
 911 _cairo_matrix_transformed_circle_major_axis (const cairo_matrix_t *matrix,
 912                                              double radius)
 913 {
 914     double  a, b, c, d, f, g, h, i, j;
 915
 916     if (_cairo_matrix_has_unity_scale (matrix))
 917         return radius;
 918
 919     _cairo_matrix_get_affine (matrix,
 920                               &a, &b,
 921                               &c, &d,
 922                               NULL, NULL);
 923
 924     i = a*a + b*b;
 925     j = c*c + d*d;
 926
 927     f = 0.5 * (i + j);
 928     g = 0.5 * (i - j);
 929     h = a*c + b*d;
 930
 931     return radius * sqrt (f + hypot (g, h));
 932
 933     /*
 934      * we don't need the minor axis length, which is
 935      * double min = radius * sqrt (f - sqrt (g*g+h*h));
 936      */
 937 }
 938
 939 static const pixman_transform_t pixman_identity_transform = {{
 940         {1 << 16,        0,       0},
 941         {       0, 1 << 16,       0},
 942         {       0,       0, 1 << 16}
 943     }};
 944
 945 static cairo_status_t
 946 _cairo_matrix_to_pixman_matrix (const cairo_matrix_t    *matrix,
 947                                 pixman_transform_t      *pixman_transform,
 948                                 double xc,
 949                                 double yc)
 950 {
 951     cairo_matrix_t inv;
 952     unsigned max_iterations;
 953
 954     pixman_transform->matrix[0][0] = _cairo_fixed_16_16_from_double (matrix->xx);
 955     pixman_transform->matrix[0][1] = _cairo_fixed_16_16_from_double (matrix->xy);
 956     pixman_transform->matrix[0][2] = _cairo_fixed_16_16_from_double (matrix->x0);
 957
 958     pixman_transform->matrix[1][0] = _cairo_fixed_16_16_from_double (matrix->yx);
 959     pixman_transform->matrix[1][1] = _cairo_fixed_16_16_from_double (matrix->yy);
 960     pixman_transform->matrix[1][2] = _cairo_fixed_16_16_from_double (matrix->y0);
 961
 962     pixman_transform->matrix[2][0] = 0;
 963     pixman_transform->matrix[2][1] = 0;
 964     pixman_transform->matrix[2][2] = 1 << 16;
 965
 966     /* The conversion above breaks cairo's translation invariance:
 967      * a translation of (a, b) in device space translates to
 968      * a translation of (xx * a + xy * b, yx * a + yy * b)
 969      * for cairo, while pixman uses rounded versions of xx ... yy.
 970      * This error increases as a and b get larger.
 971      *
 972      * To compensate for this, we fix the point (xc, yc) in pattern
 973      * space and adjust pixman's transform to agree with cairo's at
 974      * that point.
 975      */
 976
 977     if (_cairo_matrix_has_unity_scale (matrix))
 978         return CAIRO_STATUS_SUCCESS;
 979
 980     if (unlikely (fabs (matrix->xx) > PIXMAN_MAX_INT ||
 981                   fabs (matrix->xy) > PIXMAN_MAX_INT ||
 982                   fabs (matrix->x0) > PIXMAN_MAX_INT ||
 983                   fabs (matrix->yx) > PIXMAN_MAX_INT ||
 984                   fabs (matrix->yy) > PIXMAN_MAX_INT ||
 985                   fabs (matrix->y0) > PIXMAN_MAX_INT))
 986     {
 987         return _cairo_error (CAIRO_STATUS_INVALID_MATRIX);
 988     }
 989
 990     /* Note: If we can't invert the transformation, skip the adjustment. */
 991     inv = *matrix;
 992     if (cairo_matrix_invert (&inv) != CAIRO_STATUS_SUCCESS)
 993         return CAIRO_STATUS_SUCCESS;
 994
 995     /* find the pattern space coordinate that maps to (xc, yc) */
 996     max_iterations = 5;
 997     do {
 998         double x,y;
 999         pixman_vector_t vector;
1000         cairo_fixed_16_16_t dx, dy;
1001
1002         vector.vector[0] = _cairo_fixed_16_16_from_double (xc);
1003         vector.vector[1] = _cairo_fixed_16_16_from_double (yc);
1004         vector.vector[2] = 1 << 16;
1005
1006         /* If we can't transform the reference point, skip the adjustment. */
1007         if (! pixman_transform_point_3d (pixman_transform, &vector))
1008             return CAIRO_STATUS_SUCCESS;
1009
1010         x = pixman_fixed_to_double (vector.vector[0]);
1011         y = pixman_fixed_to_double (vector.vector[1]);
1012         cairo_matrix_transform_point (&inv, &x, &y);
1013
1014         /* Ideally, the vector should now be (xc, yc).
1015          * We can now compensate for the resulting error.
1016          */
1017         x -= xc;
1018         y -= yc;
1019         cairo_matrix_transform_distance (matrix, &x, &y);
1020         dx = _cairo_fixed_16_16_from_double (x);
1021         dy = _cairo_fixed_16_16_from_double (y);
1022         pixman_transform->matrix[0][2] -= dx;
1023         pixman_transform->matrix[1][2] -= dy;
1024
1025         if (dx == 0 && dy == 0)
1026             return CAIRO_STATUS_SUCCESS;
1027     } while (--max_iterations);
1028
1029     /* We didn't find an exact match between cairo and pixman, but
1030      * the matrix should be mostly correct */
1031     return CAIRO_STATUS_SUCCESS;
1032 }
1033
1034 static inline double
1035 _pixman_nearest_sample (double d)
1036 {
1037     return ceil (d - .5);
1038 }
1039
1040 /**
1041  * _cairo_matrix_is_pixman_translation:
1042  * @matrix: a matrix
1043  * @filter: the filter to be used on the pattern transformed by @matrix
1044  * @x_offset: the translation in the X direction
1045  * @y_offset: the translation in the Y direction
1046  *
1047  * Checks if @matrix translated by (x_offset, y_offset) can be
1048  * represented using just an offset (within the range pixman can
1049  * accept) and an identity matrix.
1050  *
1051  * Passing a non-zero value in x_offset/y_offset has the same effect
1052  * as applying cairo_matrix_translate(matrix, x_offset, y_offset) and
1053  * setting x_offset and y_offset to 0.
1054  *
1055  * Upon return x_offset and y_offset contain the translation vector if
1056  * the return value is %TRUE. If the return value is %FALSE, they will
1057  * not be modified.
1058  *
1059  * Return value: %TRUE if @matrix can be represented as a pixman
1060  * translation, %FALSE otherwise.
1061  **/
1062 cairo_bool_t
1063 _cairo_matrix_is_pixman_translation (const cairo_matrix_t     *matrix,
1064                                      cairo_filter_t            filter,
1065                                      int                      *x_offset,
1066                                      int                      *y_offset)
1067 {
1068     double tx, ty;
1069
1070     if (!_cairo_matrix_is_translation (matrix))
1071         return FALSE;
1072
1073     if (matrix->x0 == 0. && matrix->y0 == 0.)
1074         return TRUE;
1075
1076     tx = matrix->x0 + *x_offset;
1077     ty = matrix->y0 + *y_offset;
1078
1079     if (filter == CAIRO_FILTER_FAST || filter == CAIRO_FILTER_NEAREST) {
1080         tx = _pixman_nearest_sample (tx);
1081         ty = _pixman_nearest_sample (ty);
1082     } else if (tx != floor (tx) || ty != floor (ty)) {
1083         return FALSE;
1084     }
1085
1086     if (fabs (tx) > PIXMAN_MAX_INT || fabs (ty) > PIXMAN_MAX_INT)
1087         return FALSE;
1088
1089     *x_offset = _cairo_lround (tx);
1090     *y_offset = _cairo_lround (ty);
1091     return TRUE;
1092 }
1093
1094 /**
1095  * _cairo_matrix_to_pixman_matrix_offset:
1096  * @matrix: a matrix
1097  * @filter: the filter to be used on the pattern transformed by @matrix
1098  * @xc: the X coordinate of the point to fix in pattern space
1099  * @yc: the Y coordinate of the point to fix in pattern space
1100  * @out_transform: the transformation which best approximates @matrix
1101  * @x_offset: the translation in the X direction
1102  * @y_offset: the translation in the Y direction
1103  *
1104  * This function tries to represent @matrix translated by (x_offset,
1105  * y_offset) as a %pixman_transform_t and an translation.
1106  *
1107  * Passing a non-zero value in x_offset/y_offset has the same effect
1108  * as applying cairo_matrix_translate(matrix, x_offset, y_offset) and
1109  * setting x_offset and y_offset to 0.
1110  *
1111  * If it is possible to represent the matrix with an identity
1112  * %pixman_transform_t and a translation within the valid range for
1113  * pixman, this function will set @out_transform to be the identity,
1114  * @x_offset and @y_offset to be the translation vector and will
1115  * return %CAIRO_INT_STATUS_NOTHING_TO_DO. Otherwise it will try to
1116  * evenly divide the translational component of @matrix between
1117  * @out_transform and (@x_offset, @y_offset).
1118  *
1119  * Upon return x_offset and y_offset contain the translation vector.
1120  *
1121  * Return value: %CAIRO_INT_STATUS_NOTHING_TO_DO if the out_transform
1122  * is the identity, %CAIRO_STATUS_INVALID_MATRIX if it was not
1123  * possible to represent @matrix as a pixman_transform_t without
1124  * overflows, %CAIRO_STATUS_SUCCESS otherwise.
1125  **/
1126 cairo_status_t
1127 _cairo_matrix_to_pixman_matrix_offset (const cairo_matrix_t     *matrix,
1128                                        cairo_filter_t            filter,
1129                                        double                    xc,
1130                                        double                    yc,
1131                                        pixman_transform_t       *out_transform,
1132                                        int                      *x_offset,
1133                                        int                      *y_offset)
1134 {
1135     cairo_bool_t is_pixman_translation;
1136
1137     is_pixman_translation = _cairo_matrix_is_pixman_translation (matrix,
1138                                                                  filter,
1139                                                                  x_offset,
1140                                                                  y_offset);
1141
1142     if (is_pixman_translation) {
1143         *out_transform = pixman_identity_transform;
1144         return CAIRO_INT_STATUS_NOTHING_TO_DO;
1145     } else {
1146         cairo_matrix_t m;
1147
1148         m = *matrix;
1149         cairo_matrix_translate (&m, *x_offset, *y_offset);
1150         if (m.x0 != 0.0 || m.y0 != 0.0) {
1151             double tx, ty, norm;
1152             int i, j;
1153
1154             /* pixman also limits the [xy]_offset to 16 bits so evenly
1155              * spread the bits between the two.
1156              *
1157              * To do this, find the solutions of:
1158              *   |x| = |x*m.xx + y*m.xy + m.x0|
1159              *   |y| = |x*m.yx + y*m.yy + m.y0|
1160              *
1161              * and select the one whose maximum norm is smallest.
1162              */
1163             tx = m.x0;
1164             ty = m.y0;
1165             norm = MAX (fabs (tx), fabs (ty));
1166
1167             for (i = -1; i < 2; i+=2) {
1168                 for (j = -1; j < 2; j+=2) {
1169                     double x, y, den, new_norm;
1170
1171                     den = (m.xx + i) * (m.yy + j) - m.xy * m.yx;
1172                     if (fabs (den) < DBL_EPSILON)
1173                         continue;
1174
1175                     x = m.y0 * m.xy - m.x0 * (m.yy + j);
1176                     y = m.x0 * m.yx - m.y0 * (m.xx + i);
1177
1178                     den = 1 / den;
1179                     x *= den;
1180                     y *= den;
1181
1182                     new_norm = MAX (fabs (x), fabs (y));
1183                     if (norm > new_norm) {
1184                         norm = new_norm;
1185                         tx = x;
1186                         ty = y;
1187                     }
1188                 }
1189             }
1190
1191             tx = floor (tx);
1192             ty = floor (ty);
1193             *x_offset = -tx;
1194             *y_offset = -ty;
1195             cairo_matrix_translate (&m, tx, ty);
1196         } else {
1197             *x_offset = 0;
1198             *y_offset = 0;
1199         }
1200
1201         return _cairo_matrix_to_pixman_matrix (&m, out_transform, xc, yc);
1202     }
1203 }