shape ops work in progress
first 100,000 random cubic/cubic intersections working
git-svn-id: http://skia.googlecode.com/svn/trunk@7380 2bbb7eff-a529-9590-31e7-b0007b416f81
diff --git a/experimental/Intersection/QuadraticUtilities.cpp b/experimental/Intersection/QuadraticUtilities.cpp
index 8b02b4c..475ec13 100644
--- a/experimental/Intersection/QuadraticUtilities.cpp
+++ b/experimental/Intersection/QuadraticUtilities.cpp
@@ -5,6 +5,7 @@
* found in the LICENSE file.
*/
#include "QuadraticUtilities.h"
+#include "SkTypes.h"
#include <math.h>
/*
@@ -20,12 +21,37 @@
*/
+
+int add_valid_ts(double s[], int realRoots, double* t) {
+ int foundRoots = 0;
+ for (int index = 0; index < realRoots; ++index) {
+ double tValue = s[index];
+ if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) {
+ if (approximately_less_than_zero(tValue)) {
+ tValue = 0;
+ } else if (approximately_greater_than_one(tValue)) {
+ tValue = 1;
+ }
+ for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
+ if (approximately_equal(t[idx2], tValue)) {
+ goto nextRoot;
+ }
+ }
+ t[foundRoots++] = tValue;
+ }
+nextRoot:
+ ;
+ }
+ return foundRoots;
+}
+
// note: caller expects multiple results to be sorted smaller first
// note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting
// analysis of the quadratic equation, suggesting why the following looks at
// the sign of B -- and further suggesting that the greatest loss of precision
// is in b squared less two a c
-int quadraticRoots(double A, double B, double C, double t[2]) {
+int quadraticRootsValidT(double A, double B, double C, double t[2]) {
+#if 0
B *= 2;
double square = B * B - 4 * A * C;
if (approximately_negative(square)) {
@@ -61,9 +87,64 @@
t[0] = ratio;
}
}
+#else
+ double s[2];
+ int realRoots = quadraticRootsReal(A, B, C, s);
+ int foundRoots = add_valid_ts(s, realRoots, t);
+#endif
return foundRoots;
}
+// unlike quadratic roots, this does not discard real roots <= 0 or >= 1
+int quadraticRootsReal(const double A, const double B, const double C, double s[2]) {
+ if (approximately_zero(A)) {
+ if (approximately_zero(B)) {
+ s[0] = 0;
+ return C == 0;
+ }
+ s[0] = -C / B;
+ return 1;
+ }
+ /* normal form: x^2 + px + q = 0 */
+ const double p = B / (2 * A);
+ const double q = C / A;
+ const double p2 = p * p;
+#if 0
+ double D = AlmostEqualUlps(p2, q) ? 0 : p2 - q;
+ if (D <= 0) {
+ if (D < 0) {
+ return 0;
+ }
+ s[0] = -p;
+ SkDebugf("[%d] %1.9g\n", 1, s[0]);
+ return 1;
+ }
+ double sqrt_D = sqrt(D);
+ s[0] = sqrt_D - p;
+ s[1] = -sqrt_D - p;
+ SkDebugf("[%d] %1.9g %1.9g\n", 2, s[0], s[1]);
+ return 2;
+#else
+ if (!AlmostEqualUlps(p2, q) && p2 < q) {
+ return 0;
+ }
+ double sqrt_D = 0;
+ if (p2 > q) {
+ sqrt_D = sqrt(p2 - q);
+ }
+ s[0] = sqrt_D - p;
+ s[1] = -sqrt_D - p;
+#if 0
+ if (AlmostEqualUlps(s[0], s[1])) {
+ SkDebugf("[%d] %1.9g\n", 1, s[0]);
+ } else {
+ SkDebugf("[%d] %1.9g %1.9g\n", 2, s[0], s[1]);
+ }
+#endif
+ return 1 + !AlmostEqualUlps(s[0], s[1]);
+#endif
+}
+
static double derivativeAtT(const double* quad, double t) {
double a = t - 1;
double b = 1 - 2 * t;