Gradient fixes

* Port fix for bug 7685 from pixman. Patch by Carl Worth

* Add projective version of radial gradient code.

* Make sure that all Pict*Gradient types have PictGradient as prefix,
  since code in various places relies on that.

1 files changed, 173 insertions, 30 deletions

diff --git a/fb/fbcompose.c b/fb/fbcompose.c
index 0faf7836d..24b552e0b 100644
--- a/fb/fbcompose.c
+++ b/fb/fbcompose.c
@@ -3148,13 +3148,128 @@ static void fbFetchSourcePict(PicturePtr pict, int x, int y, int width, CARD32 *
             }
         }
     } else {
+
+/*
+ * In the radial gradient problem we are given two circles (c₁,r₁) and
+ * (c₂,r₂) that define the gradient itself. Then, for any point p, we
+ * must compute the value(s) of t within [0.0, 1.0] representing the
+ * circle(s) that would color the point.
+ *
+ * There are potentially two values of t since the point p can be
+ * colored by both sides of the circle, (which happens whenever one
+ * circle is not entirely contained within the other).
+ *
+ * If we solve for a value of t that is outside of [0.0, 1.0] then we
+ * use the extend mode (NONE, REPEAT, REFLECT, or PAD) to map to a
+ * value within [0.0, 1.0].
+ *
+ * Here is an illustration of the problem:
+ *
+ *              p₂
+ *           p  •
+ *           •   ╲
+ *        ·       ╲r₂
+ *  p₁ ·           ╲
+ *  •              θ╲
+ *   ╲             ╌╌•
+ *    ╲r₁        ·   c₂
+ *    θ╲    ·
+ *    ╌╌•
+ *      c₁
+ *
+ * Given (c₁,r₁), (c₂,r₂) and p, we must find an angle θ such that two
+ * points p₁ and p₂ on the two circles are collinear with p. Then, the
+ * desired value of t is the ratio of the length of p₁p to the length
+ * of p₁p₂.
+ *
+ * So, we have six unknown values: (p₁x, p₁y), (p₂x, p₂y), θ and t.
+ * We can also write six equations that constrain the problem:
+ *
+ * Point p₁ is a distance r₁ from c₁ at an angle of θ:
+ *
+ *	1. p₁x = c₁x + r₁·cos θ
+ *	2. p₁y = c₁y + r₁·sin θ
+ *
+ * Point p₂ is a distance r₂ from c₂ at an angle of θ:
+ *
+ *	3. p₂x = c₂x + r2·cos θ
+ *	4. p₂y = c₂y + r2·sin θ
+ *
+ * Point p lies at a fraction t along the line segment p₁p₂:
+ *
+ *	5. px = t·p₂x + (1-t)·p₁x
+ *	6. py = t·p₂y + (1-t)·p₁y
+ *
+ * To solve, first subtitute 1-4 into 5 and 6:
+ *
+ * px = t·(c₂x + r₂·cos θ) + (1-t)·(c₁x + r₁·cos θ)
+ * py = t·(c₂y + r₂·sin θ) + (1-t)·(c₁y + r₁·sin θ)
+ *
+ * Then solve each for cos θ and sin θ expressed as a function of t:
+ *
+ * cos θ = (-(c₂x - c₁x)·t + (px - c₁x)) / ((r₂-r₁)·t + r₁)
+ * sin θ = (-(c₂y - c₁y)·t + (py - c₁y)) / ((r₂-r₁)·t + r₁)
+ *
+ * To simplify this a bit, we define new variables for several of the
+ * common terms as shown below:
+ *
+ *              p₂
+ *           p  •
+ *           •   ╲
+ *        ·  ┆    ╲r₂
+ *  p₁ ·     ┆     ╲
+ *  •     pdy┆      ╲
+ *   ╲       ┆       •c₂
+ *    ╲r₁    ┆   ·   ┆
+ *     ╲    ·┆       ┆cdy
+ *      •╌╌╌╌┴╌╌╌╌╌╌╌┘
+ *    c₁  pdx   cdx
+ *
+ * cdx = (c₂x - c₁x)
+ * cdy = (c₂y - c₁y)
+ *  dr =  r₂-r₁
+ * pdx =  px - c₁x
+ * pdy =  py - c₁y
+ *
+ * Note that cdx, cdy, and dr do not depend on point p at all, so can
+ * be pre-computed for the entire gradient. The simplifed equations
+ * are now:
+ *
+ * cos θ = (-cdx·t + pdx) / (dr·t + r₁)
+ * sin θ = (-cdy·t + pdy) / (dr·t + r₁)
+ *
+ * Finally, to get a single function of t and eliminate the last
+ * unknown θ, we use the identity sin²θ + cos²θ = 1. First, square
+ * each equation, (we knew a quadratic was coming since it must be
+ * possible to obtain two solutions in some cases):
+ *
+ * cos²θ = (cdx²t² - 2·cdx·pdx·t + pdx²) / (dr²·t² + 2·r₁·dr·t + r₁²)
+ * sin²θ = (cdy²t² - 2·cdy·pdy·t + pdy²) / (dr²·t² + 2·r₁·dr·t + r₁²)
+ *
+ * Then add both together, set the result equal to 1, and express as a
+ * standard quadratic equation in t of the form At² + Bt + C = 0
+ *
+ * (cdx² + cdy² - dr²)·t² - 2·(cdx·pdx + cdy·pdy + r₁·dr)·t + (pdx² + pdy² - r₁²) = 0
+ *
+ * In other words:
+ *
+ * A = cdx² + cdy² - dr²
+ * B = -2·(pdx·cdx + pdy·cdy + r₁·dr)
+ * C = pdx² + pdy² - r₁²
+ *
+ * And again, notice that A does not depend on p, so can be
+ * precomputed. From here we just use the quadratic formula to solve
+ * for t:
+ *
+ * t = (-2·B ± ⎷(B² - 4·A·C)) / 2·A
+ */
         /* radial or conical */
         Bool affine = TRUE;
         double cx = 1.;
         double cy = 0.;
         double cz = 0.;
-        double rx = x;
-        double ry = y;
+	double rx = x + 0.5;
+	double ry = y + 0.5;
         double rz = 1.;
 
         if (pict->transform) {
@@ -3176,25 +3291,38 @@ static void fbFetchSourcePict(PicturePtr pict, int x, int y, int width, CARD32 *
         }
 
         if (pGradient->type == SourcePictTypeRadial) {
+	    PictRadialGradient *radial;
+	    radial = &pGradient->radial;
             if (affine) {
-                rx -= pGradient->radial.fx;
-                ry -= pGradient->radial.fy;
-
                 while (buffer < end) {
-		    double b, c, det, s;
-		    
 		    if (!mask || *mask++ & maskBits)
 		    {
-			xFixed_48_16  t;
+			double pdx, pdy;
+			double B, C;
+			double det;
+			double c1x = radial->c1.x / 65536.0;
+			double c1y = radial->c1.y / 65536.0;
+			double r1  = radial->c1.radius / 65536.0;
+                        xFixed_48_16 t;
+
+			pdx = rx - c1x;
+			pdy = ry - c1y;
+
+			B = -2 * (  pdx * radial->cdx
+				  + pdy * radial->cdy
+				  + r1 * radial->dr);
+			C = (pdx * pdx + pdy * pdy - r1 * r1);
+
+                        det = (B * B) - (4 * radial->A * C);
+			if (det < 0.0)
+			    det = 0.0;
+
+			if (radial->A < 0)
+			    t = (xFixed_48_16) ((- B - sqrt(det)) / (2.0 * radial->A) * 65536);
+			else
+			    t = (xFixed_48_16) ((- B + sqrt(det)) / (2.0 * radial->A) * 65536);
 			
-			b = 2*(rx*pGradient->radial.dx + ry*pGradient->radial.dy);
-			c = -(rx*rx + ry*ry);
-			det = (b * b) - (4 * pGradient->radial.a * c);
-			s = (-b + sqrt(det))/(2. * pGradient->radial.a);
-			
-			t = (xFixed_48_16)((s*pGradient->radial.m + pGradient->radial.b)*65536);
-			
-			WRITE(buffer, _gradient_walker_pixel (&walker, t));
+ 			WRITE(buffer, _gradient_walker_pixel (&walker, t));
 		    }
 		    ++buffer;
 		    
@@ -3202,35 +3330,50 @@ static void fbFetchSourcePict(PicturePtr pict, int x, int y, int width, CARD32 *
                     ry += cy;
                 }
             } else {
+		/* projective */
                 while (buffer < end) {
-                    double x, y;
-                    double b, c, det, s;
-
 		    if (!mask || *mask++ & maskBits)
 		    {
-			xFixed_48_16  t;
-			
+			double pdx, pdy;
+			double B, C;
+			double det;
+			double c1x = radial->c1.x / 65536.0;
+			double c1y = radial->c1.y / 65536.0;
+			double r1  = radial->c1.radius / 65536.0;
+                        xFixed_48_16 t;
+			double x, y;
+
 			if (rz != 0) {
 			    x = rx/rz;
 			    y = ry/rz;
 			} else {
 			    x = y = 0.;
 			}
-			x -= pGradient->radial.fx;
-			y -= pGradient->radial.fy;
-			b = 2*(x*pGradient->radial.dx + y*pGradient->radial.dy);
-			c = -(x*x + y*y);
-			det = (b * b) - (4 * pGradient->radial.a * c);
-			s = (-b + sqrt(det))/(2. * pGradient->radial.a);
-			t = (xFixed_48_16)((s*pGradient->radial.m + pGradient->radial.b)*65536);
 			
-			WRITE(buffer, _gradient_walker_pixel (&walker, t));
+			pdx = x - c1x;
+			pdy = y - c1y;
+
+			B = -2 * (  pdx * radial->cdx
+				  + pdy * radial->cdy
+				  + r1 * radial->dr);
+			C = (pdx * pdx + pdy * pdy - r1 * r1);
+
+                        det = (B * B) - (4 * radial->A * C);
+			if (det < 0.0)
+			    det = 0.0;
+
+			if (radial->A < 0)
+			    t = (xFixed_48_16) ((- B - sqrt(det)) / (2.0 * radial->A) * 65536);
+			else
+			    t = (xFixed_48_16) ((- B + sqrt(det)) / (2.0 * radial->A) * 65536);
+			
+ 			WRITE(buffer, _gradient_walker_pixel (&walker, t));
 		    }
 		    ++buffer;
 		    
                     rx += cx;
                     ry += cy;
-                    rz += cz;
+		    rz += cz;
                 }
             }
         } else /* SourcePictTypeConical */ {