[cairo] [PATCH] Improve computation of filter scale factors

[otaylor at redhat.com]

From: "Owen W. Taylor" <otaylor at fishsoup.net>

The filter scale factors were being computed based on the lengths
of the unit vectors under the pattern transform. This works for
pure scales, and pure rotations, but produces a completely wrong
result for, say, scaling X by 3 and then rotating by 90 degrees.

Instead compute the filter scale factors by computing the bounding
box of a unit circle transformed into pattern coordinates, and
then rescaling that so that scale_x * scale_y = det(pattern matrix)
---
 src/cairo-image-source.c | 66 ++++++++++++++++++++++++++++++++++++++++++------
 1 file changed, 58 insertions(+), 8 deletions(-)

diff --git a/src/cairo-image-source.c b/src/cairo-image-source.c
index 661bc10..a9c7991 100644
--- a/src/cairo-image-source.c
+++ b/src/cairo-image-source.c
@@ -555,18 +555,68 @@ _pixman_image_set_properties (pixman_image_t *pixman_image,
     else
     {
 	double scale_x, scale_y;
+	double alpha_x, alpha_y, det;
 	int shrink_x, shrink_y;
 	pixman_filter_t pixman_filter;
 	pixman_kernel_t pixman_kernel_sample, pixman_kernel_reconstruct;
 
-	/* Compute scale factors as the length of basis vectors transformed by
-	 * the pattern matrix. These scale factors are from user to pattern space,
-	 * and as such they are greater than 1.0 for downscaling and less than 1.0
-	 * for upscaling.
-	 * TODO: this approach may not be completely correct if the matrix
-	 * contains a skew component. */
-	scale_x = hypot (pattern->matrix.xx, pattern->matrix.yx);
-	scale_y = hypot (pattern->matrix.yx, pattern->matrix.yy);
+        /* Determining scale factors: a square pixel in the destination comes
+	 * from an an area in the source image which is generally a
+	 * parallelogram.
+	 *
+	 *  . . . .               . . . . . + . .
+	 *  . + + .         -\    . . . . / .\. .
+	 *  . + + .         -/    . . .  /. . + .
+	 *  . . . .               . . . + . ./. .
+	 *                        . . .  \. / . .
+	 *                        . . . . + . . .
+	 *
+	 * with filtering we'll sample a larger area in the source which will be
+	 * roughly ellipsoidal. But the separable convolution filter support in
+	 * Pixman requirest the filter be axis aligned - so we need to come up
+	 * with an axis aligned area that has a good correspondance. Consider
+	 * the image of the unit circle under the pattern matrix - it will be an
+	 * ellipse. We find the axis-aligned ellipse with the same bounding box,
+	 * then scale it down until it has the same area as the original
+	 * transformed ellipse. (This area is π D, where D is the determinant of
+	 * the pattern matrix.) The two axes of the axis-aligned ellipse are our
+	 * scale factors.
+	 *
+	 * Pattern matrix: [X] = [Axx Axy] * [x]
+	 *                 [Y]   [Ayx Ayy]   [y]
+	 * Inverse:        [x] = 1/D * [ Ayy -Axy] [X]
+	 *                 [y] =       [-Ayx  Axx] [Y]
+	 *
+	 *
+	 * The unit circle in the destination: x² + y² = 1
+	 * Substituting source coordinates:
+	 *
+	 *              (Ayy•X - Axy•Y)² + (-Ayx•X + Axx•Y)² = D²
+	 *          (Ayy² + Ayx²)X² - 2 (Ayy•Axy + Ayx•Axx) XY + (Axx² + Axy²) Y² = D²
+	 * Define:      αy²•X²      -          2 β•XY          +     αx²•Y²       = D²
+	 *
+	 * To find the bounding box, consider this as a quadratic equation in
+	 * Y. The X extremes are where the quadratic equation has one solution -
+	 * i.e. the discriminant b²-4ac is 0.
+	 *
+	 *   4β²X0² - 4αx² (αy²•X0² - D²) = 0
+	 *
+	 * Solving for X0: X0² = D²αx²/(αx²αy² - β²)
+	 * Similarly:      Y0² = D²αy²/(αx²αy² - β²)
+	 *
+	 * And X0/Y0 = αx/αy. Want to find Sx/Sy in the same ratio such that Sx•Sy=D
+	 *
+	 * (αx/αy)Sy • Sy = D => Sy² = D•αy/αx; Sx² = D•αx/αy.
+	 */
+	alpha_x = hypot (pattern->matrix.xx, pattern->matrix.xy);
+	alpha_y = hypot (pattern->matrix.yx, pattern->matrix.yy);
+	det = pattern->matrix.xx * pattern->matrix.yy - pattern->matrix.xy * pattern->matrix.yx;
+
+	/* αx, αy are non-zero, since they can only be zero if Axx=Axy=0 or Ayx=Ayy=0,
+	 * assuming quality implementation of hypot(), and then D would be zero,
+	 * which is not allowed for a pattern matrix. */
+	scale_x = sqrt (fabs (det) * alpha_x / alpha_y);
+	scale_y = sqrt (fabs (det) * alpha_y / alpha_x);
 
 	/* Use convolution filtering if the transformation shrinks the image
 	 * by more than half a pixel */
-- 
1.9.0