2010-09-09 Kenneth Russell <kbr@google.com>
[WebKit.git] / WebCore / platform / graphics / gpu / LoopBlinnMathUtils.cpp
1 /*
2  * Copyright (C) 2010 Google Inc. All rights reserved.
3  *
4  * Redistribution and use in source and binary forms, with or without
5  * modification, are permitted provided that the following conditions
6  * are met:
7  *
8  * 1.  Redistributions of source code must retain the above copyright
9  *     notice, this list of conditions and the following disclaimer.
10  * 2.  Redistributions in binary form must reproduce the above copyright
11  *     notice, this list of conditions and the following disclaimer in the
12  *     documentation and/or other materials provided with the distribution.
13  *
14  * THIS SOFTWARE IS PROVIDED BY APPLE AND ITS CONTRIBUTORS "AS IS" AND ANY
15  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
16  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
17  * DISCLAIMED. IN NO EVENT SHALL APPLE OR ITS CONTRIBUTORS BE LIABLE FOR ANY
18  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
19  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
20  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
21  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
22  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
23  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24  */
25
26 #include "config.h"
27
28 #include "LoopBlinnMathUtils.h"
29
30 #include "FloatPoint.h"
31 #include "MathExtras.h"
32 #include <algorithm>
33 #include <string.h> // for memcpy
34
35 namespace WebCore {
36 namespace LoopBlinnMathUtils {
37
38 namespace {
39
40 // Utility functions local to this file.
41 int orientation(const FloatPoint& p1,
42                 const FloatPoint& p2,
43                 const FloatPoint& p3)
44 {
45     float crossProduct = (p2.y() - p1.y()) * (p3.x() - p2.x()) - (p3.y() - p2.y()) * (p2.x() - p1.x());
46     return (crossProduct < 0.0f) ? -1 : ((crossProduct > 0.0f) ? 1 : 0);
47 }
48
49 bool edgeEdgeTest(const FloatSize& v0Delta,
50                   const FloatPoint& v0,
51                   const FloatPoint& u0,
52                   const FloatPoint& u1)
53 {
54     // This edge to edge test is based on Franlin Antonio's gem: "Faster
55     // Line Segment Intersection", in Graphics Gems III, pp. 199-202.
56     float ax = v0Delta.width();
57     float ay = v0Delta.height();
58     float bx = u0.x() - u1.x();
59     float by = u0.y() - u1.y();
60     float cx = v0.x() - u0.x();
61     float cy = v0.y() - u0.y();
62     float f = ay * bx - ax * by;
63     float d = by * cx - bx * cy;
64     if ((f > 0 && d >= 0 && d <= f) || (f < 0 && d <= 0 && d >= f)) {
65         float e = ax * cy - ay * cx;
66
67         // This additional test avoids reporting coincident edges, which
68         // is the behavior we want.
69         if (approxEqual(e, 0) || approxEqual(f, 0) || approxEqual(e, f))
70             return false;
71
72         if (f > 0)
73             return e >= 0 && e <= f;
74
75         return e <= 0 && e >= f;
76     }
77     return false;
78 }
79
80 bool edgeAgainstTriangleEdges(const FloatPoint& v0,
81                               const FloatPoint& v1,
82                               const FloatPoint& u0,
83                               const FloatPoint& u1,
84                               const FloatPoint& u2)
85 {
86     FloatSize delta = v1 - v0;
87     // Test edge u0, u1 against v0, v1.
88     if (edgeEdgeTest(delta, v0, u0, u1))
89         return true;
90     // Test edge u1, u2 against v0, v1.
91     if (edgeEdgeTest(delta, v0, u1, u2))
92         return true;
93     // Test edge u2, u1 against v0, v1.
94     if (edgeEdgeTest(delta, v0, u2, u0))
95         return true;
96     return false;
97 }
98
99 // A roundoff factor in the cubic classification and texture coordinate
100 // generation algorithms. It primarily determines the handling of corner
101 // cases during the classification process. Be careful when adjusting it;
102 // it has been determined empirically to work well. When changing it, you
103 // should look in particular at shapes that contain quadratic curves and
104 // ensure they still look smooth. Once pixel tests are running against this
105 // algorithm, they should provide sufficient coverage to ensure that
106 // adjusting the constant won't break anything.
107 const float Epsilon = 5.0e-4f;
108
109 } // anonymous namespace
110
111 // Exported routines
112
113 float roundToZero(float val)
114 {
115     if (val < Epsilon && val > -Epsilon)
116         return 0;
117     return val;
118 }
119
120 bool approxEqual(const FloatPoint& v0, const FloatPoint& v1)
121 {
122     return (v0 - v1).diagonalLengthSquared() < Epsilon * Epsilon;
123 }
124
125 bool approxEqual(const FloatPoint3D& v0, const FloatPoint3D& v1)
126 {
127     return (v0 - v1).lengthSquared() < Epsilon * Epsilon;
128 }
129
130 bool approxEqual(float f0, float f1)
131 {
132     return fabsf(f0 - f1) < Epsilon;
133 }
134
135 bool linesIntersect(const FloatPoint& p1,
136                     const FloatPoint& q1,
137                     const FloatPoint& p2,
138                     const FloatPoint& q2)
139 {
140     return (orientation(p1, q1, p2) != orientation(p1, q1, q2)
141             && orientation(p2, q2, p1) != orientation(p2, q2, q1));
142 }
143
144 bool pointInTriangle(const FloatPoint& point,
145                      const FloatPoint& a,
146                      const FloatPoint& b,
147                      const FloatPoint& c)
148 {
149     // Algorithm from http://www.blackpawn.com/texts/pointinpoly/default.html
150     float x0 = c.x() - a.x();
151     float y0 = c.y() - a.y();
152     float x1 = b.x() - a.x();
153     float y1 = b.y() - a.y();
154     float x2 = point.x() - a.x();
155     float y2 = point.y() - a.y();
156
157     float dot00 = x0 * x0 + y0 * y0;
158     float dot01 = x0 * x1 + y0 * y1;
159     float dot02 = x0 * x2 + y0 * y2;
160     float dot11 = x1 * x1 + y1 * y1;
161     float dot12 = x1 * x2 + y1 * y2;
162     float denominator = dot00 * dot11 - dot01 * dot01;
163     if (!denominator)
164         // Triangle is zero-area. Treat query point as not being inside.
165         return false;
166     // Compute
167     float inverseDenominator = 1.0f / denominator;
168     float u = (dot11 * dot02 - dot01 * dot12) * inverseDenominator;
169     float v = (dot00 * dot12 - dot01 * dot02) * inverseDenominator;
170
171     return (u > 0.0f) && (v > 0.0f) && (u + v < 1.0f);
172 }
173
174 bool trianglesOverlap(const FloatPoint& a1,
175                       const FloatPoint& b1,
176                       const FloatPoint& c1,
177                       const FloatPoint& a2,
178                       const FloatPoint& b2,
179                       const FloatPoint& c2)
180 {
181     // Derived from coplanar_tri_tri() at
182     // http://jgt.akpeters.com/papers/ShenHengTang03/tri_tri.html ,
183     // simplified for the 2D case and modified so that overlapping edges
184     // do not report overlapping triangles.
185
186     // Test all edges of triangle 1 against the edges of triangle 2.
187     if (edgeAgainstTriangleEdges(a1, b1, a2, b2, c2)
188         || edgeAgainstTriangleEdges(b1, c1, a2, b2, c2)
189         || edgeAgainstTriangleEdges(c1, a1, a2, b2, c2))
190         return true;
191     // Finally, test if tri1 is totally contained in tri2 or vice versa.
192     // The paper above only performs the first two point-in-triangle tests.
193     // Because we define that triangles sharing a vertex or edge don't
194     // overlap, we must perform additional tests to see whether one
195     // triangle is contained in the other.
196     if (pointInTriangle(a1, a2, b2, c2)
197         || pointInTriangle(a2, a1, b1, c1)
198         || pointInTriangle(b1, a2, b2, c2)
199         || pointInTriangle(b2, a1, b1, c1)
200         || pointInTriangle(c1, a2, b2, c2)
201         || pointInTriangle(c2, a1, b1, c1))
202         return true;
203     return false;
204 }
205
206 namespace {
207
208 // Helper routines for public XRay queries below. All of this code
209 // originated in Skia; see include/core/ and src/core/, SkScalar.h and
210 // SkGeometry.{cpp,h}.
211
212 const float NearlyZeroConstant = (1.0f / (1 << 12));
213
214 bool nearlyZero(float x, float tolerance = NearlyZeroConstant)
215 {
216     ASSERT(tolerance > 0.0f);
217     return ::fabsf(x) < tolerance;
218 }
219
220 // Linearly interpolate between a and b, based on t.
221 // If t is 0, return a; if t is 1, return b; else interpolate.
222 // t must be [0..1].
223 float interpolate(float a, float b, float t)
224 {
225     ASSERT(t >= 0 && t <= 1);
226     return a + (b - a) * t;
227 }
228
229 float evaluateCubic(float controlPoint0, float controlPoint1, float controlPoint2, float controlPoint3, float t)
230 {
231     ASSERT(t >= 0 && t <= 1);
232
233     if (!t)
234         return controlPoint0;
235
236     float ab = interpolate(controlPoint0, controlPoint1, t);
237     float bc = interpolate(controlPoint1, controlPoint2, t);
238     float cd = interpolate(controlPoint2, controlPoint3, t);
239     float abc = interpolate(ab, bc, t);
240     float bcd = interpolate(bc, cd, t);
241     return interpolate(abc, bcd, t);
242 }
243
244 // Evaluates the point on the source cubic specified by t, 0 <= t <= 1.0.
245 FloatPoint evaluateCubicAt(const FloatPoint cubic[4], float t)
246 {
247     return FloatPoint(evaluateCubic(cubic[0].x(), cubic[1].x(), cubic[2].x(), cubic[3].x(), t),
248                       evaluateCubic(cubic[0].y(), cubic[1].y(), cubic[2].y(), cubic[3].y(), t));
249 }
250
251 bool xRayCrossesMonotonicCubic(const XRay& xRay, const FloatPoint cubic[4], bool& ambiguous)
252 {
253     ambiguous = false;
254
255     // Find the minimum and maximum y of the extrema, which are the
256     // first and last points since this cubic is monotonic
257     float minY = std::min(cubic[0].y(), cubic[3].y());
258     float maxY = std::max(cubic[0].y(), cubic[3].y());
259
260     if (xRay.y() == cubic[0].y()
261         || xRay.y() < minY
262         || xRay.y() > maxY) {
263         // The query line definitely does not cross the curve
264         ambiguous = (xRay.y() == cubic[0].y());
265         return false;
266     }
267
268     const bool pointAtExtremum = (xRay.y() == cubic[3].y());
269
270     float minX = std::min(std::min(std::min(cubic[0].x(), cubic[1].x()),
271                                    cubic[2].x()),
272                           cubic[3].x());
273     if (xRay.x() < minX) {
274         // The query line definitely crosses the curve
275         ambiguous = pointAtExtremum;
276         return true;
277     }
278
279     float maxX = std::max(std::max(std::max(cubic[0].x(), cubic[1].x()),
280                                    cubic[2].x()),
281                           cubic[3].x());
282     if (xRay.x() > maxX)
283         // The query line definitely does not cross the curve
284         return false;
285
286     // Do a binary search to find the parameter value which makes y as
287     // close as possible to the query point. See whether the query
288     // line's origin is to the left of the associated x coordinate.
289
290     // MaxIterations is chosen as the number of mantissa bits for a float,
291     // since there's no way we are going to get more precision by
292     // iterating more times than that.
293     const int MaxIterations = 23;
294     FloatPoint evaluatedPoint;
295     int iter = 0;
296     float upperT;
297     float lowerT;
298     // Need to invert direction of t parameter if cubic goes up
299     // instead of down
300     if (cubic[3].y() > cubic[0].y()) {
301         upperT = 1;
302         lowerT = 0;
303     } else {
304         upperT = 0;
305         lowerT = 1;
306     }
307     do {
308         float t = 0.5f * (upperT + lowerT);
309         evaluatedPoint = evaluateCubicAt(cubic, t);
310         if (xRay.y() > evaluatedPoint.y())
311             lowerT = t;
312         else
313             upperT = t;
314     } while (++iter < MaxIterations && !nearlyZero(evaluatedPoint.y() - xRay.y()));
315
316     // FIXME: once we have more regression tests for this code,
317     // determine whether this should be using a fuzzy test.
318     if (xRay.x() <= evaluatedPoint.x()) {
319         ambiguous = pointAtExtremum;
320         return true;
321     }
322     return false;
323 }
324
325 // Divides the numerator by the denominator safely for the case where
326 // the result must lie in the range (0..1). Result indicates whether
327 // the result is valid.
328 bool safeUnitDivide(float numerator, float denominator, float& ratio)
329 {
330     if (numerator < 0) {
331         // Make the "numerator >= denominator" check below work.
332         numerator = -numerator;
333         denominator = -denominator;
334     }
335     if (!numerator || !denominator || numerator >= denominator)
336         return false;
337     float r = numerator / denominator;
338     if (isnan(r))
339         return false;
340     ASSERT(r >= 0 && r < 1);
341     if (!r) // catch underflow if numerator <<<< denominator
342         return false;
343     ratio = r;
344     return true;
345 }
346
347 // From Numerical Recipes in C.
348 //
349 //   q = -1/2 (b + sign(b) sqrt[b*b - 4*a*c])
350 //   x1 = q / a
351 //   x2 = c / q
352 //
353 // Returns the number of real roots of the equation [0..2]. Roots are
354 // returned in sorted order, smaller root first.
355 int findUnitQuadRoots(float a, float b, float c, float roots[2])
356 {
357     if (!a)
358         return safeUnitDivide(-c, b, roots[0]) ? 1 : 0;
359
360     float discriminant = b*b - 4*a*c;
361     if (discriminant < 0 || isnan(discriminant)) // complex roots
362         return 0;
363     discriminant = sqrtf(discriminant);
364
365     float q = (b < 0) ? -(b - discriminant) / 2 : -(b + discriminant) / 2;
366     int numberOfRoots = 0;
367     if (safeUnitDivide(q, a, roots[numberOfRoots]))
368         ++numberOfRoots;
369     if (safeUnitDivide(c, q, roots[numberOfRoots]))
370         ++numberOfRoots;
371     if (numberOfRoots == 2) {
372         // Seemingly have two roots. Check for equality and sort.
373         if (roots[0] == roots[1])
374             return 1;
375         if (roots[0] > roots[1])
376             std::swap(roots[0], roots[1]);
377     }
378     return numberOfRoots;
379 }
380
381 // Cubic'(t) = pt^2 + qt + r, where
382 //   p = 3(-a + 3(b - c) + d)
383 //   q = 6(a - 2b + c)
384 //   r = 3(b - a)
385 // Solve for t, keeping only those that fit between 0 < t < 1.
386 int findCubicExtrema(float a, float b, float c, float d, float tValues[2])
387 {
388     // Divide p, q, and r by 3 to simplify the equations.
389     float p = d - a + 3*(b - c);
390     float q = 2*(a - b - b + c);
391     float r = b - a;
392
393     return findUnitQuadRoots(p, q, r, tValues);
394 }
395
396 void interpolateCubicCoords(float controlPoint0, float controlPoint1, float controlPoint2, float controlPoint3, float* dst, float t)
397 {
398     float ab = interpolate(controlPoint0, controlPoint1, t);
399     float bc = interpolate(controlPoint1, controlPoint2, t);
400     float cd = interpolate(controlPoint2, controlPoint3, t);
401     float abc = interpolate(ab, bc, t);
402     float bcd = interpolate(bc, cd, t);
403     float abcd = interpolate(abc, bcd, t);
404
405     dst[0] = controlPoint0;
406     dst[2] = ab;
407     dst[4] = abc;
408     dst[6] = abcd;
409     dst[8] = bcd;
410     dst[10] = cd;
411     dst[12] = controlPoint3;
412 }
413
414 #ifndef NDEBUG
415 bool isUnitInterval(float x)
416 {
417     return x > 0 && x < 1;
418 }
419 #endif
420
421 void chopCubicAtTValues(const FloatPoint src[4], FloatPoint dst[], const float tValues[], int roots)
422 {
423 #ifndef NDEBUG
424     for (int i = 0; i < roots - 1; ++i) {
425         ASSERT(isUnitInterval(tValues[i]));
426         ASSERT(isUnitInterval(tValues[i+1]));
427         ASSERT(tValues[i] < tValues[i+1]);
428     }
429 #endif
430
431     if (!roots) {
432         // nothing to chop
433         for (int j = 0; j < 4; ++j)
434             dst[j] = src[j];
435         return;
436     }
437
438     float t = tValues[0];
439     FloatPoint tmp[4];
440     for (int j = 0; j < 4; ++j)
441         tmp[j] = src[j];
442
443     for (int i = 0; i < roots; ++i) {
444         chopCubicAt(tmp, dst, t);
445         if (i == roots - 1)
446             break;
447
448         dst += 3;
449         // Make tmp contain the remaining cubic (after the first chop).
450         for (int j = 0; j < 4; ++j)
451             tmp[j] = dst[j];
452
453         // Watch out for the case that the renormalized t isn't in range.
454         if (!safeUnitDivide(tValues[i+1] - tValues[i], 1.0f - tValues[i], t)) {
455             // If it isn't, just create a degenerate cubic.
456             dst[4] = dst[5] = dst[6] = tmp[3];
457             break;
458         }
459     }
460 }
461
462 void flattenDoubleCubicYExtrema(FloatPoint coords[7])
463 {
464     coords[2].setY(coords[3].y());
465     coords[4].setY(coords[3].y());
466 }
467
468 int chopCubicAtYExtrema(const FloatPoint src[4], FloatPoint dst[10])
469 {
470     float tValues[2];
471     int roots = findCubicExtrema(src[0].y(), src[1].y(), src[2].y(), src[3].y(), tValues);
472
473     chopCubicAtTValues(src, dst, tValues, roots);
474     if (roots) {
475         // we do some cleanup to ensure our Y extrema are flat
476         flattenDoubleCubicYExtrema(&dst[0]);
477         if (roots == 2)
478             flattenDoubleCubicYExtrema(&dst[3]);
479     }
480     return roots;
481 }
482
483 } // anonymous namespace
484
485 // Public cubic operations.
486
487 void chopCubicAt(const FloatPoint src[4], FloatPoint dst[7], float t)
488 {
489     ASSERT(t >= 0 && t <= 1);
490
491     float output[14];
492     interpolateCubicCoords(src[0].x(), src[1].x(), src[2].x(), src[3].x(), &output[0], t);
493     interpolateCubicCoords(src[0].y(), src[1].y(), src[2].y(), src[3].y(), &output[1], t);
494     for (int i = 0; i < 7; i++)
495         dst[i].set(output[2 * i], output[2 * i + 1]);
496 }
497
498 // Public XRay queries.
499
500 bool xRayCrossesLine(const XRay& xRay, const FloatPoint pts[2], bool& ambiguous)
501 {
502     ambiguous = false;
503
504     // Determine quick discards.
505     // Consider query line going exactly through point 0 to not
506     // intersect, for symmetry with xRayCrossesMonotonicCubic.
507     if (xRay.y() == pts[0].y()) {
508         ambiguous = true;
509         return false;
510     }
511     if (xRay.y() < pts[0].y() && xRay.y() < pts[1].y())
512         return false;
513     if (xRay.y() > pts[0].y() && xRay.y() > pts[1].y())
514         return false;
515     if (xRay.x() > pts[0].x() && xRay.x() > pts[1].x())
516         return false;
517     // Determine degenerate cases
518     if (nearlyZero(pts[0].y() - pts[1].y()))
519         return false;
520     if (nearlyZero(pts[0].x() - pts[1].x())) {
521         // We've already determined the query point lies within the
522         // vertical range of the line segment.
523         if (xRay.x() <= pts[0].x()) {
524             ambiguous = (xRay.y() == pts[1].y());
525             return true;
526         }
527         return false;
528     }
529     // Ambiguity check
530     if (xRay.y() == pts[1].y()) {
531         if (xRay.x() <= pts[1].x()) {
532             ambiguous = true;
533             return true;
534         }
535         return false;
536     }
537     // Full line segment evaluation
538     float deltaY = pts[1].y() - pts[0].y();
539     float deltaX = pts[1].x() - pts[0].x();
540     float slope = deltaY / deltaX;
541     float b = pts[0].y() - slope * pts[0].x();
542     // Solve for x coordinate at y = xRay.y()
543     float x = (xRay.y() - b) / slope;
544     return xRay.x() <= x;
545 }
546
547 int numXRayCrossingsForCubic(const XRay& xRay, const FloatPoint cubic[4], bool& ambiguous)
548 {
549     int numCrossings = 0;
550     FloatPoint monotonicCubics[10];
551     int numMonotonicCubics = 1 + chopCubicAtYExtrema(cubic, monotonicCubics);
552     ambiguous = false;
553     FloatPoint* monotonicCubicsPointer = &monotonicCubics[0];
554     for (int i = 0; i < numMonotonicCubics; ++i) {
555         if (xRayCrossesMonotonicCubic(xRay, monotonicCubicsPointer, ambiguous))
556             ++numCrossings;
557         if (ambiguous)
558             return 0;
559         monotonicCubicsPointer += 3;
560     }
561     return numCrossings;
562 }
563
564 } // namespace LoopBlinnMathUtils
565 } // namespace WebCore