Hi Jim.
> I'm not sure how frequently we run into case 2 in practice, and the
> code
> being simpler and based on a mechanism that has survived testing
> pretty
> well is a big win. I'm willing to go with this. Plus, if we find ways
> to make the Curve methods go faster then it will help a lot of shape
> types.
>
> Have you used more than one case to test the performance differential,
> or did you find a single case that gave the intended result and then
> run
> just that one case N times? If so, then perhaps the performance
> differential is largely an idiosyncrasy of the particular test case?
>
> Either way, I think the existing patch is a good compromise and
> possibly
> close to being a fairly reliable win depending on what kind of
> performance testing you did.
For case (1) I call contains/intersects on 2x2 rectangles that
are outside of the bounding box of the curve. For (2) I created
N 2x2 rectangles with their top lefts on the curve. (3) was like
(1) but with the rectangles inside the curve. I always used the
same curve: (0, 0, 0, 1000, 1000, 1000, 1000, 0) so the rectangles
were extremely small compared to the curve. Changing the size of
the rectangles improved the performance for (2) by about 50%, since
fewer subdivisions had to be made to find that some point of the
curve was inside the rectangle.
I also wrote another test ({intersect,contain}sBigTest) that tests
both correctness and performance. It does the former by creating an
Area(CubicCurve2D) object and using the results of its
contains/intersects methods as reference results. Like the above tests
it uses a fixed curve for every call, but the rectangles have random dimensions
and positions. On 10000000 rectangles, the intersection version of this
test runs in about 8 seconds using the old intersects algorithm, and 1.4s
using the new one. Even better, this test finds disagreements between
the old CubicCurve2D.intersects and Area.intersects, but not between
the new CubicCurve2D.intersects and Area.intersects.
This is very good. Almost too good, which is why I've attached the file
with the tests.
I have yet to think about the optimizations, so I'll send an e-mail later
about those.
Regards,
Denis.
----- Original Message -----
> Hi Denis,
>
> On 1/12/2011 2:33 PM, Denis Lila wrote:
> > Hi Jim.
> >
> >> I replaced that test with
> >>
> >> if (!getBounds2D().contains(x, y, w, h)) {
> >> return false;
> >> }
> >>
> >> and that fixed both the bug and the performance of (1) without
> >> making
> >> (3) worse.
> >
> > It turns out that that test actually slows it down a bit. It was
> > helping
> > us when we were using the PathIterator, but after the introduction
> > of
> > the rectCrossings method you suggested, it's faster in all cases
> > (even (1))
> > to just do the general thing all the time. I've updated the webrev
> > with this.
>
> Great news! After reading the previous email I was going to ask that
> question, but you've already done the work.
>
> > Eliminating the PathIterator also had a large impact on the
> > intersect(rect)
> > method. Now it's about 20% faster in cases (1) and (3). (2) is still
> > slower
> > but only by a factor of 2-3. Here's the webrev for it:
> > http://icedtea.classpath.org/~dlila/webrevs/intersectsFix/webrev/
>
>
> The following is for further optimization experiments...
>
> > I think the reason (2) is slow is because rectCrossingsForCubic
> > recursively
> > searches through subdivisions of the input curve starting at t=0 and
> > in increasing t. Do you think it would be worth it to switch the
> > order
> > of the recursive calls depending on the distances of the two
> > subdivided curves relative to the rectangle (so that perhaps we
> > would get
> > to the subdivided curve that crosses one of the rectangle sides
> > sooner)?
>
> Not sure - how would you gauge "distance to the rectangle"? How about
> this quick test:
>
> if ((y0 <= y1) == (ymid <= ymin)) {
> // Either rightside up and first half is likely a fast rejection
> // or upside down and first half is possibly a reject
> do second half first
> } else {
> do first half first
> }
>
> Either way, it only saves a few tests for the branch that isn't taken.
> What if we optimized the fast rejection cases (which would make all
> test
> cases go faster) by trying to do some trivial min/max sharing for the
> Y
> case rejections. Minimally if the first Y rejection test finds that y0
> >= ymax then there is no need to test y0 <= ymin in the second set of
> rejection tests, so the following would cost no more than what we do
> now:
>
> // Assuming ymin strictly < ymax
> if (y0 >= ymax) {
> if (all others >= ymax) {
> return crossings;
> }
> // No need to do ymin testing since the first test would fail
> } else if (all <= ymin) {
> return crossings;
> }
>
> I'm not sure how many tests it might save in practice, though, but it
> would never cost any more tests (and the compiler can't optimize it
> away
> since it doesn't understand that we can require ymin<ymax as part of
> the
> method contract).
>
> Another solution:
>
> if (y0 <= y1) {
> // y0 is above if y1 is above
> // y1 is below if y0 is below
> test y1, yc1 and yc0 above
> test y0, yc0 and yc1 below
> } else {
> // reverse assumptions as above
> test y0, yc0 and yc1 above
> test y1, yc1, and yc0 below
> }
>
> (Note that it leads off every case above with a test of y0 or y1 since
> those tests are testing 2 rejection points against the rectangle, but
> the yc0 and yc1 tests only test a single point against the rectangle.)
> It only eliminates a total of 1 test, though since you still have to
> test y0 against y1. You can take it another step further by comparing
> yc0 against yc1:
>
> if (y0 <= y1) {
> // y0 is above if y1 is above
> // y1 is below if y0 is below
> if (yc0 <= yc1) {
> // similar assumptions about yc0,yc1 ordering
> test y1, yc1 above
> test y0, yc0 below
> } else {
> test y1, yc0 above
> test y0, yc1 below
> }
> } else {
> /* similar */
> }
>
> One downside with these "ordering the control points" approaches,
> though, is that the minimum number of tests in the rejection portion
> may
> go up, even if the max number goes down. It may simply be trading off
> average for consistency.
>
> Another idea: Once a curve is monotonic in Y then we can do very fast
> rejections. It might be worth testing for monotonicity (in Y mainly)
> along with the above/below rejections and switch to a faster monotonic
> method when that case occurs:
>
> if (y0 <= yc0 && yc0 <= yc1 && yc1 <= y1) {
> return rectCrossingsForMonotonicCubic(crossings, ...);
> } else if (reverse monotonic tests) {
> return 0 - rectCrossingsForMonotonicCubic(0 - crossings,
> reverse curve);
> }
> // Standard y rejection tests...etc
>
> ...jim
import java.awt.BasicStroke;
import java.awt.Color;
import java.awt.Graphics;
import java.awt.Graphics2D;
import java.awt.Panel;
import java.awt.Point;
import java.awt.Rectangle;
import java.awt.geom.Area;
import java.awt.geom.Line2D;
import java.awt.geom.Path2D;
import java.awt.geom.Point2D;
import java.awt.geom.QuadCurve2D;
import java.awt.geom.Rectangle2D;
import java.awt.image.BufferedImage;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.LinkedList;
import java.util.List;
import java.util.Random;
import javax.swing.JFrame;
import javax.swing.SwingUtilities;
import sun.awt.geom.Curve;
import static java.lang.System.out;
import static java.lang.Math.abs;
import static java.lang.Math.max;
import static java.lang.Math.ulp;
public class CubicSolver {
public static int solveCubicOld(double eqn[], double res[]) {
if (res == eqn) {
// Copy the eqn so that we don't clobber it with the
// roots.
eqn = new double[4];
System.arraycopy(res, 0, eqn, 0, 4);
}
// From Numerical Recipes, 5.6, Quadratic and Cubic Equations
double d = eqn[3];
if (d == 0.0) {
// The cubic has degenerated to quadratic (or line or ...).
return QuadCurve2D.solveQuadratic(eqn, res);
}
double a = eqn[2] / d;
double b = eqn[1] / d;
double c = eqn[0] / d;
int roots = 0;
double Q = (a * a - 3.0 * b) / 9.0;
double R = (2.0 * a * a * a - 9.0 * a * b + 27.0 * c) / 54.0;
double R2 = R * R;
double Q3 = Q * Q * Q;
a = a / 3.0;
if (R2 < Q3) {
double theta = Math.acos(R / Math.sqrt(Q3));
Q = -2.0 * Math.sqrt(Q);
res[roots++] = Q * Math.cos(theta / 3.0) - a;
res[roots++] = Q * Math.cos((theta + Math.PI * 2.0)/ 3.0) - a;
res[roots++] = Q * Math.cos((theta - Math.PI * 2.0)/ 3.0) - a;
} else {
boolean neg = (R < 0.0);
double S = Math.sqrt(R2 - Q3);
if (neg) {
R = -R;
}
double A = Math.pow(R + S, 1.0 / 3.0);
if (!neg) {
A = -A;
}
double B = (A == 0.0) ? 0.0 : (Q / A);
res[roots++] = (A + B) - a;
}
if (roots == 3) {
fixRoots(res, eqn);
}
return roots;
}
private static boolean iszero(double x, double err) {
return within(x, 0, err);
}
private static boolean within(final double x, final double y, final double err) {
final double d = y - x;
return (d <= err && d >= -err);
}
private static void fixRoots(double res[], double eqn[]) {
final double EPSILON = 1E-5;
for (int i = 0; i < 3; i++) {
double t = res[i];
if (Math.abs(t) < EPSILON) {
res[i] = findZero(t, 0, eqn);
} else if (Math.abs(t - 1) < EPSILON) {
res[i] = findZero(t, 1, eqn);
}
}
}
public static int solveCubicNew(double eqn[], double res[]) {
// From Graphics Gems:
// http://tog.acm.org/resources/GraphicsGems/gems/Roots3And4.c
final double d = eqn[3];
if (d == 0) {
return QuadCurve2D.solveQuadratic(eqn, res);
}
/* normal form: x^3 + Ax^2 + Bx + C = 0 */
final double A = eqn[2] / d;
final double B = eqn[1] / d;
final double C = eqn[0] / d;
// substitute x = y - A/3 to eliminate quadratic term:
// x^3 +px + q = 0
double sq_A = A * A;
double p = 1.0/3 * (-1.0/3 * sq_A + B);
double q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);
/* use Cardano's formula */
double cb_p = p * p * p;
double D = q * q + cb_p;
final double sub = 1.0/3 * A;
int num;
if (D < 0) { /* Casus irreducibilis: three real solutions */
final double phi = 1.0/3 * Math.acos(-q / Math.sqrt(-cb_p));
final double t = 2 * Math.sqrt(-p);
if (res == eqn) {
eqn = Arrays.copyOf(res, 4);
}
res[ 0 ] = ( t * Math.cos(phi));
res[ 1 ] = (-t * Math.cos(phi + Math.PI / 3));
res[ 2 ] = (-t * Math.cos(phi - Math.PI / 3));
num = 3;
for (int i = 0; i < num; ++i)
res[ i ] -= sub;
} else {
final double sqrt_D = Math.sqrt(D);
final double u = Math.cbrt(sqrt_D - q);
final double v = - Math.cbrt(sqrt_D + q);
final double uv = u+v;
res[ 0 ] = uv - sub;
num = 1;
final double err = 1200000000*ulp(abs(uv) + abs(sub));
if (iszero(D, err) || within(u, v, err)) {
if (res == eqn) {
eqn = Arrays.copyOf(res, 4);
}
res[1] = -(uv / 2) - sub;
num = 2;
}
}
if (num == 3 || num == 2) {
num = fixRoots(eqn, res, num);
}
if (num > 2 && (res[2] == res[1] || res[2] == res[0])) {
num--;
}
if (num > 1 && res[1] == res[0]) {
res[1] = res[--num]; // Copies res[2] to res[1] if needed
}
return num;
}
private static int fixRoots(double[] eqn, double[] res, final int num) {
double[] intervals = {eqn[1], 2*eqn[2], 3*eqn[3]};
int critCount = QuadCurve2D.solveQuadratic(intervals, intervals);
if (critCount == 2 && intervals[0] == intervals[1]) {
critCount--;
}
if (critCount == 2 && intervals[0] > intervals[1]) {
double tmp = intervals[0];
intervals[0] = intervals[1];
intervals[1] = tmp;
}
double xe = getRootUpperBound(eqn);
xe += max(ulp(xe), 1);
double x0 = -xe;
double fx0 = eqn[3] > 0 ? -1 : 1;
double fxe = -fx0;
if (num == 3) {
Arrays.sort(res, 0, num);
if (critCount == 2) {
res[0] = refineRootWithHint(eqn, x0, intervals[0], res[0]);
res[1] = refineRootWithHint(eqn, intervals[0], intervals[1], res[1]);
res[2] = refineRootWithHint(eqn, intervals[1], xe, res[2]);
return 3;
} else if (critCount == 1) {
double x1 = intervals[0];
double fx1 = solveEqn(eqn, 3, x1);
if (oppositeSigns(fx0, fx1)) {
res[0] = bisectRootWithHint(eqn, x0, x1, res[0]);
} else if (oppositeSigns(fx1, fxe)) {
res[0] = bisectRootWithHint(eqn, x1, xe, res[2]);
} else /* c1 must be 0 */ {
res[0] = x1;
}
} else if (critCount == 0) {
res[0] = bisectRootWithHint(eqn, x0, xe, res[1]);
}
} else if (num == 2 && critCount == 2) {
double badRoot = res[1];
double goodRoot = res[0];
double x1 = intervals[0];
double x2 = intervals[1];
double x = abs(x1 - goodRoot) > abs(x2 - goodRoot) ? x1 : x2;
double fx = solveEqn(eqn, 3, x);
if (iszero(fx, 10000000*ulp(x))) {
double badRootVal = solveEqn(eqn, 3, badRoot);
res[1] = abs(badRootVal) < abs(fx) ? badRoot : x;
return 2;
}
} // else there can only be one root - goodRoot, and it is already in res[0]
return 1;
}
// use newton's method.
private static double refineRootWithHint(double[] eqn, double min, double max, double t) {
if (!inInterval(t, min, max)) {
return t;
}
double[] deriv = {eqn[1], 2*eqn[2], 3*eqn[3]};
double slope;
double origt = t;
for (int i = 0; i < 3; i++) {
slope = solveEqn(deriv, 2, t);
double y = solveEqn(eqn, 3, t);
double delta = - (y / slope);
double newt = t + delta;
if (slope == 0 || y == 0 || t == newt) {
break;
}
t = newt;
}
if (within(t, origt, 1000*ulp(origt)) && inInterval(t, min, max)) {
return t;
}
return origt;
}
private static boolean inInterval(double t, double min, double max) {
return min <= t && t <= max;
}
public static int solveCubic(double eqn[], double res[]) {
// From Graphics Gems:
// http://tog.acm.org/resources/GraphicsGems/gems/Roots3And4.c
final double d = eqn[3];
if (d == 0) {
return QuadCurve2D.solveQuadratic(eqn, res);
}
if (res == eqn) {
// Copy the eqn so that we don't clobber it with the roots.
eqn = new double[4];
System.arraycopy(res, 0, eqn, 0, 4);
}
/* normal form: x^3 + Ax^2 + Bx + C = 0 */
final double A = eqn[2] / d;
final double B = eqn[1] / d;
final double C = eqn[0] / d;
// substitute x = y - A/3 to eliminate quadratic term:
// x^3 +px + q = 0
double sq_A = A * A;
double p = 1.0/3 * (-1.0/3 * sq_A + B);
double q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);
/* use Cardano's formula */
double cb_p = p * p * p;
double D = q * q + cb_p;
int num;
if (D < 0) { /* Casus irreducibilis: three real solutions */
final double phi = 1.0/3 * Math.acos(-q / Math.sqrt(-cb_p));
final double t = 2 * Math.sqrt(-p);
res[ 0 ] = ( t * Math.cos(phi));
res[ 1 ] = (-t * Math.cos(phi + Math.PI / 3));
res[ 2 ] = (-t * Math.cos(phi - Math.PI / 3));
num = 3;
} else { /* one real solution */
final double sqrt_D = Math.sqrt(D);
final double u = Math.cbrt(sqrt_D - q);
final double v = -Math.cbrt(sqrt_D + q);
res[ 0 ] = u + v;
num = 1;
if (within(u, v, 1e-5)) {
res[1] = -res[0] / 2;
num = 2;
}
}
/* resubstitute */
final double sub = 1.0/3 * A;
for (int i = 0; i < num; ++i)
res[ i ] -= sub;
return num;
}
private static double bisectRootWithHint(double[] eqn, double x0, double xe, double hint) {
double delta1 = Math.min(abs(hint - x0) / 64, 0.0625);
double delta2 = Math.min(abs(hint - xe) / 64, 0.0625);
double x02 = hint - delta1;
double xe2 = hint + delta2;
double fx02 = solveEqn(eqn, 3, x02);
double fxe2 = solveEqn(eqn, 3, xe2);
while(oppositeSigns(fx02, fxe2)) {
if (x02 >= xe2) {
return x02;
}
x0 = x02;
xe = xe2;
delta1 /= 64;
delta2 /= 64;
x02 = hint - delta1;
xe2 = hint + delta2;
fx02 = solveEqn(eqn, 3, x02);
fxe2 = solveEqn(eqn, 3, xe2);
}
if (fx02 == 0) {
return x02;
}
if (fxe2 == 0) {
return xe2;
}
return bisectRoot(eqn, x0, xe);
}
private static double bisectRoot(double[] eqn, double x0, double xe) {
double fx0 = solveEqn(eqn, 3, x0);
double m = x0 + (xe - x0) / 2;
double fm = solveEqn(eqn, 3, m);
while (m != x0 && m != xe) {
if (fm == 0) {
return m;
}
if (oppositeSigns(fx0, fm)) {
xe = m;
} else {
fx0 = fm;
x0 = m;
}
m = x0 + (xe-x0)/2;
fm = solveEqn(eqn, 3, m);
}
return m;
}
private static double getRootUpperBound(double[] eqn) {
final double d = eqn[3];
final double a = eqn[2];
final double b = eqn[1];
final double c = eqn[0];
double M = 1 + max(max(abs(a), abs(b)), abs(c)) / abs(d);
M += max(ulp(M), 1);
return M;
}
private static boolean oppositeSigns(double x1, double x2) {
return (x1 < 0 && x2 > 0) || (x1 > 0 && x2 < 0);
}
private static double solveEqn(double eqn[], int order, double t) {
double v = eqn[order];
while (--order >= 0) {
v = v * t + eqn[order];
}
return v;
}
private static double findZero(double t, double target, double eqn[]) {
double slopeqn[] = {eqn[1], 2*eqn[2], 3*eqn[3]};
double slope;
double origdelta = 0;
double origt = t;
while (true) {
slope = solveEqn(slopeqn, 2, t);
if (slope == 0) {
// At a local minima - must return
return t;
}
double y = solveEqn(eqn, 3, t);
if (y == 0) {
// Found it! - return it
return t;
}
// assert(slope != 0 && y != 0);
double delta = - (y / slope);
// assert(delta != 0);
if (origdelta == 0) {
origdelta = delta;
}
if (t < target) {
if (delta < 0) return t;
} else if (t > target) {
if (delta > 0) return t;
} else { /* t == target */
return (delta > 0
? (target + java.lang.Double.MIN_VALUE)
: (target - java.lang.Double.MIN_VALUE));
}
double newt = t + delta;
if (t == newt) {
// The deltas are so small that we aren't moving...
return t;
}
if (delta * origdelta < 0) {
// We have reversed our path.
int tag = (origt < t
? getTag(target, origt, t)
: getTag(target, t, origt));
if (tag != INSIDE) {
// Local minima found away from target - return the middle
return (origt + t) / 2;
}
// Local minima somewhere near target - move to target
// and let the slope determine the resulting t.
t = target;
} else {
t = newt;
}
}
}
private static final int BELOW = -2;
private static final int LOWEDGE = -1;
private static final int INSIDE = 0;
private static final int HIGHEDGE = 1;
private static final int ABOVE = 2;
private static int getTag(double coord, double low, double high) {
if (coord <= low) {
return (coord < low ? BELOW : LOWEDGE);
}
if (coord >= high) {
return (coord > high ? ABOVE : HIGHEDGE);
}
return INSIDE;
}
private static double[] makeEqnWithRoots(double r1, double r2, double r3) {
double eqn[] = new double[4];
eqn[3] = 1.0;
eqn[2] = - (r1 + r2 + r3);
eqn[1] = r1*r2 + r2*r3 + r3*r1;
eqn[0] = - (r1 * r2 * r3);
return eqn;
}
private static int oldTooMany = 0;
private static int oldTooFew = 0;
private static int newTooFew = 0;
private static int newTooMany = 0;
private static List<Double> oldDiffs = new ArrayList<Double>();
private static List<Double> newDiffs = new ArrayList<Double>();
private static int oldWins = 0;
private static int newWins = 0;
private static double rd(Random r, double min, double max) {
return r.nextDouble() * (max - min) + min;
}
private static void cubicSolvePerfTest(int n) {
double[][] eqns = new double[n][];
double[][] oldres = new double[n][4];
int[] oldNumRoots = new int[n];
double[][] newres = new double[n][4];
int[] newNumRoots = new int[n];
double min = -1000;
double max = 1000;
Random r = new Random();
for(int i = 0; i < n; i++) {
eqns[i] = makeEqnWithRoots(rd(r, min, max), rd(r, min, max), rd(r, min, max));
}
long t1 = System.nanoTime();
for(int i = 0; i < n; i++) {
newNumRoots[i] = solveCubicNew(eqns[i], newres[i]);
}
long t2 = System.nanoTime();
for(int i = 0; i < n; i++) {
oldNumRoots[i] = solveCubic(eqns[i], oldres[i]);
}
long t3 = System.nanoTime();
System.out.println("Test2, solve NewNewNew took: " + (t2 - t1));
System.out.println("Test2, solve NewNew took: " + (t3 - t2));
System.out.println(Arrays.toString(newres[abs(r.nextInt())%n]));
System.out.println(Arrays.toString(oldres[abs(r.nextInt())%n]));
compStats(eqns, oldres, oldNumRoots, newres, newNumRoots);
}
private static Point2D evalCurve(CubicCurve2D c, double t) {
double x1 = c.getX1();
double x2 = c.getX2();
double cx1 = c.getCtrlX1();
double cx2 = c.getCtrlX2();
double y1 = c.getY1();
double y2 = c.getY2();
double cy1 = c.getCtrlY1();
double cy2 = c.getCtrlY2();
final double t1 = 1-t;
final double x = x1*t1*t1*t1 + 3*cx1*t1*t1*t + 3*cx2*t1*t*t + x2*t*t*t;
final double y = y1*t1*t1*t1 + 3*cy1*t1*t1*t + 3*cy2*t1*t*t + y2*t*t*t;
return new Point2D.Double(x, y);
}
private static void compStats(double[][] eqns, double[][] oldres, int[] oldNumRoots, double[][] newres, int[] newNumRoots) {
int disagreements = 0;
int num = oldNumRoots.length;
List<Double> olddiffs = new LinkedList<Double>();
List<Double> newdiffs = new LinkedList<Double>();
for (int i = 0; i < num; i++) {
int oldNum = oldNumRoots[i];
int newNum = newNumRoots[i];
double[] eqn = eqns[i];
if (oldNum == newNum) {
for (int j = 0; j < oldNum; j++) {
double oldroot = oldres[i][j];
double newroot = newres[i][j];
olddiffs.add(abs(solveEqn(eqn, 3, oldroot)));
newdiffs.add(abs(solveEqn(eqn, 3, newroot)));
}
} else {
disagreements++;
}
}
int same = 0, newWin = 0, oldWin = 0;
double maxoldfiff = 0, maxnewdiff = 0;
double meanolddiff = 0, meannewdiff = 0;
for (int i = 0; i < olddiffs.size(); i++) {
double o = olddiffs.get(i);
double n = newdiffs.get(i);
meanolddiff += o;
meannewdiff += n;
maxoldfiff = max(maxoldfiff, o);
maxnewdiff = max(maxnewdiff, n);
if (o < n) {
oldWin++;
} else if (n < o) {
newWin++;
} else {
same++;
}
}
out.println("Out of " + num + " equations:");
out.printf("%d disagreements on the number of roots%n", disagreements);
out.printf("sameDiffs = %d, newWins = %d, oldwins = %d%n", same, newWin, oldWin);
out.printf("maxolddiff = %.22f, maxnewdiff = %.22f%n", maxoldfiff, maxnewdiff);
out.printf("totalolddiff = %.22f, totalnewdiff = %.22f%n", meanolddiff, meannewdiff);
meanolddiff /= olddiffs.size();
meannewdiff /= newdiffs.size();
out.printf("meanolddiff = %.22f, meannewdiff = %.22f%n", meanolddiff, meannewdiff);
}
public static void trySolve3(double r1, double r2, double r3, int count) {
// (x-r1)*(x-r2)*(x-r3)
// = (xx - r1x - r2x + r1r2)*(x-r3)
// = xxx - r1xx - r2xx + r1r2x - r3xx + r1r3x + r2r3x - r1r2r3
// = xxx - (r1 + r2 + r3)xx + (r1r2 + r2r3 + r3r1)x - r1r2r3
// System.out.println("solving: (x - "+r1+") * (x - "+r2+") * (x - "+r3+") = 0");
double eqn[] = makeEqnWithRoots(r1, r2, r3);
double resOld[] = new double[4];
int n = solveCubic(eqn, resOld);
Arrays.sort(resOld, 0, n);
boolean otf = false;
if (n < count) {
otf = true;
oldTooFew++;
} else if (n > count) {
oldTooMany++;
}
double resNew[] = new double[4];
int m = solveCubicNew(eqn, resNew);
Arrays.sort(resNew, 0, m);
boolean ntf = false;
if (m < count) {
ntf = true;
newTooFew++;
} else if (m > count) {
newTooMany++;
}
if (n == count && m == count) {
double[] roots = {r1, r2, r3};
Arrays.sort(roots, 0, 3);
if (r1 == r2) {
roots[1] = roots[2];
}
double dold = 0;
double dnew = 0;
for (int i = 0; i < count; i++) {
double oldRoot = resOld[i];
double newRoot = resNew[i];
double actRoot = roots[i];
dold += abs(oldRoot - actRoot);
dnew += abs(newRoot - actRoot);
}
oldDiffs.add(dold);
newDiffs.add(dnew);
if (dold < dnew) {
System.out.printf(r1 + ", " + r2 + ", " + r3 + "\n");
oldWins++;
} else if (dnew < dold) {
newWins++;
}
}
// if (m < count) {
// System.out.printf(r1 + ", " + r2 + ", " + r3 + "\n");
// }
if (m != count) {
System.out.printf(r1 + ", " + r2 + ", " + r3 + "\n");
System.out.print("There should be " + count + " roots. new method returned "+m+" roots: ");
for (int i = 0; i < m; i++) {
System.out.print(resNew[i]+", ");
}
System.out.println();
System.out.println();
} else {
System.out.print("Good: ");
for (int i = 0; i < m; i++) {
System.out.print(resNew[i]+", ");
}
System.out.println();
}
}
private static void printStats() {
double mean = 0;
for (Double d : oldDiffs) {
mean += d;
}
mean /= oldDiffs.size();
double stdDev = 0;
for (Double d : oldDiffs) {
stdDev += (d - mean) * (d - mean);
}
stdDev /= oldDiffs.size();
stdDev = Math.sqrt(stdDev);
System.out.printf("Old mean = %f, stdDev = %f, old wins =%d%n", mean, stdDev, oldWins);
mean = 0;
for (Double d : newDiffs) {
mean += d;
}
mean /= newDiffs.size();
stdDev = 0;
for (Double d : newDiffs) {
stdDev += (d - mean) * (d - mean);
}
stdDev /= newDiffs.size();
stdDev = Math.sqrt(stdDev);
System.out.printf("New mean = %f, stdDev = %f, new wins =%d%n", mean, stdDev, newWins);
System.out.println();
System.out.println("oldTooMany : " + oldTooMany);
System.out.println("oldTooFew : " + oldTooFew );
System.out.println("newTooFew : " + newTooFew );
System.out.println("newTooMany : " + newTooMany);
}
public static void trySolve2(double r1, double r2) {
trySolve3(r1, r2, r2, 2);
trySolve3(r1, r1, r2, 2);
}
public static void trySolve1(double r) {
trySolve3(r, r, r, 1);
}
public static void trySolve(double r1, double r2, double r3) {
trySolve3(r1, r2, r3, 3);
trySolve2(r1, r2);
trySolve2(r2, r3);
trySolve2(r3, r1);
trySolve1(r1);
trySolve1(r2);
trySolve1(r3);
System.out.println("--------------------------------------------------------------");
}
private static void containsRectPerfTest(int n, int b) {
Random r = new Random();
boolean[] results = new boolean[n];
int w = 1000;
CubicCurve2D c = new CubicCurve2D.Double(0, 0, 0, w, w, w, w, 0);
// clear misses:
long clearMisses = System.nanoTime();
for (int i = 0; i < n; i++) {
results[i] = c.contains(i, -100, b, b); // should all be false
}
clearMisses = System.nanoTime() - clearMisses;
System.out.println("noopt: " + results[r.nextInt(n)]);
final double[] inputs = new double[n*2];
for (int i = 0; i < n; i++) {
Point2D p = evalCurve(c, ((double)i)/n);
inputs[2*i] = p.getX() - 3*r.nextDouble();
inputs[2*i + 1] = p.getY();
}
long notClear = System.nanoTime();
for (int i = 0; i < n; i++) {
results[i] = c.contains(inputs[2*i], inputs[2*i+1], b, b);
}
notClear = System.nanoTime() - notClear;
System.out.println("noopt: " + results[r.nextInt(n)]);
long clearInside = System.nanoTime();
for (int i = 0; i < n; i++) {
results[i] = c.contains( (((double)i)/n) * (w-4) + 2, 1, b, b);
}
clearInside = System.nanoTime() - clearInside;
System.out.println("noopt: " + results[r.nextInt(n)]);
System.out.printf("clearMissesTime for %d runs = %f%n", n, clearMisses/1000000000.0);
System.out.printf("notClearTime for %d runs = %f%n", n, notClear /1000000000.0);
System.out.printf("clearInsisdeTime for %d runs = %f%n", n, clearInside/1000000000.0);
}
private static void intersectRectPerfTest(int n, int b) {
Random r = new Random();
boolean[] results = new boolean[n];
int w = 1000;
CubicCurve2D c = new CubicCurve2D.Double(0, 0, 0, w, w, w, w, 0);
// clear misses:
long clearMisses = System.nanoTime();
for (int i = 0; i < n; i++) {
results[i] = c.intersects(i, -100, b, b); // should all be false
}
clearMisses = System.nanoTime() - clearMisses;
System.out.println("noopt: " + results[r.nextInt(n)]);
final double[] inputs = new double[n*2];
for (int i = 0; i < n; i++) {
Point2D p = evalCurve(c, ((double)i)/n);
inputs[2*i] = p.getX() - 3*r.nextDouble();
inputs[2*i + 1] = p.getY();
}
long notClear = System.nanoTime();
for (int i = 0; i < n; i++) {
results[i] = c.intersects(inputs[2*i], inputs[2*i+1], b, b);
}
notClear = System.nanoTime() - notClear;
System.out.println("noopt: " + results[r.nextInt(n)]);
long clearInside = System.nanoTime();
for (int i = 0; i < n; i++) {
results[i] = c.intersects( (((double)i)/n) * (w-4) + 2, 1, b, b);
}
clearInside = System.nanoTime() - clearInside;
System.out.println("noopt: " + results[r.nextInt(n)]);
System.out.printf("clearMissesTime for %d runs = %f%n", n, clearMisses/1000000000.0);
System.out.printf("notClearTime for %d runs = %f%n", n, notClear /1000000000.0);
System.out.printf("clearInsisdeTime for %d runs = %f%n", n, clearInside/1000000000.0);
}
public static void main(String argv[]) throws IOException {
// cubicSolvePerfTest(20000);
// int n = 1000000;
// containsRectPerfTest(n, 200);
// intersectRectPerfTest(n, 200);
out.println("iBT: " + intersectsBigTest(10000000));
out.println("cBT: " + containsBigTest(10000000));
}
private static boolean intersectsBigTest(int n) {
double b = 40;
CubicCurve2D cc = new CubicCurve2D.Double(0, 0, b, b, -b/2, b, b/2, 0);
Random r = new Random();
final double minx = -2;
final double miny = -1;
double[][] ins = new double[n][4];
for (int i = 0; i < n; i++) {
final double w = rd(r, 0, b/5);
final double h = rd(r, 0, b/5);
final double maxx = (b+1) - w;
final double maxy = (b+1) - h;
double x = rd(r, minx, maxx);
double y = rd(r, miny, maxy);
ins[i][0] = x;
ins[i][1] = y;
ins[i][2] = w;
ins[i][3] = h;
}
boolean[] results = new boolean[n];
long t1 = System.nanoTime();
for (int i = 0; i < n; i++) {
results[i] = cc.intersects(ins[i][0], ins[i][1], ins[i][2], ins[i][3]);
}
t1 = System.nanoTime() - t1;
out.println("iBT time = " + (t1 / 1000000000.0));
Area a = new Area(cc);
for (int i = 0; i < n; i++) {
if (a.intersects(ins[i][0], ins[i][1], ins[i][2], ins[i][3]) != results[i]) {
out.println("disagreement on : " + Arrays.toString(ins[i]));
return false;
}
}
return true;
}
private static boolean containsBigTest(int n) {
double b = 40;
CubicCurve2D cc = new CubicCurve2D.Double(0, 0, b, b, -b/2, b, b/2, 0);
Random r = new Random();
final double minx = -2;
final double miny = -1;
double[][] ins = new double[n][4];
for (int i = 0; i < n; i++) {
final double w = rd(r, 0, b/5);
final double h = rd(r, 0, b/5);
final double maxx = (b+1) - w;
final double maxy = (b+1) - h;
double x = rd(r, minx, maxx);
double y = rd(r, miny, maxy);
ins[i][0] = x;
ins[i][1] = y;
ins[i][2] = w;
ins[i][3] = h;
}
boolean[] results = new boolean[n];
long t1 = System.nanoTime();
for (int i = 0; i < n; i++) {
results[i] = cc.contains(ins[i][0], ins[i][1], ins[i][2], ins[i][3]);
}
t1 = System.nanoTime() - t1;
out.println("cBT time = " + (t1 / 1000000000.0));
Area a = new Area(cc);
for (int i = 0; i < n; i++) {
if (a.contains(ins[i][0], ins[i][1], ins[i][2], ins[i][3]) != results[i]) {
out.println("disagreement on : " + Arrays.toString(ins[i]));
return false;
}
}
return true;
}
// TODO: write huge brute force correctness test that uses Area's methods as a reference.
private static boolean containTest3(boolean disp) throws IOException {
float W = 150.0f;
float H = 150.0f;
float x1 = 0.0f, y1 = 0.0f;
CubicCurve2D.Float c = new CubicCurve2D.Float(x1,y1,W,y1,W,y1,W,H); // make well behaved curve
float w = (c.x2-c.x1), h = (c.y2-c.y1);
Rectangle2D r = new Rectangle2D.Float(c.x1+w/2,c.y1+w/4,w/4,h/4);
System.out.println("r2 is in? = " + c.contains(r)); // should return true, actually returns false....bug???
if (disp) {
int width = 1500;
int height = 1000;
Object[] a = getGraphics(width, height);
BufferedImage bi = (BufferedImage) a[0];
Graphics2D gi = (Graphics2D) a[1];
gi.translate(width/2, height/2);
gi.setStroke(new BasicStroke(1f/20));
gi.setColor(Color.blue);
gi.draw(new Line2D.Double(c.getX1(), c.getY1(), c.getX2(), c.getY2()));
gi.draw(c);
gi.setColor(Color.green);
gi.setColor(Color.red);
gi.draw(r);
display(bi);
}
return c.contains(r);
}
private static boolean intersectTest(boolean disp) throws IOException {
CubicCurve2D c = new CubicCurve2D.Double(50.0, 300.0,
150.0, 166.6666717529297,
238.0, 456.66668701171875,
350.0, 300.0);
Rectangle2D r = new Rectangle2D.Double(260, 300, 10, 10);
boolean shadowtouches = c.intersects(r);
boolean shadowinside = c.contains(r);
if (disp) {
int width = 1500;
int height = 1000;
Object[] a = getGraphics(width, height);
BufferedImage bi = (BufferedImage) a[0];
Graphics2D gi = (Graphics2D) a[1];
gi.translate(100, 100);
gi.setColor(Color.blue);
gi.draw(c);
gi.draw(new Line2D.Double(c.getX1(), c.getY1(), c.getX2(), c.getY2()));
gi.setColor(Color.red);
gi.draw(r);
display(bi);
}
if (!shadowtouches)
{
System.out.println(" shadowinside: "+shadowinside);
System.out.println(" shadowtouches: "+shadowtouches);
}
return shadowtouches;
}
private static boolean containTest(boolean disp) throws IOException {
// CubicCurve2D cc = new CubicCurve2D.Double(5, -2, -1, 3.55, -3, -6, 0, 3);
// Rectangle2D rect = new Rectangle2D.Double(-1, -0.25, 4, 0.05);
// CubicCurve2D cc = new CubicCurve2D.Double(0, 0, 4, -4, 3, -2, 4, 0);
CubicCurve2D cc = new CubicCurve2D.Double(0, 0, 4, -4, -2, -4, 2, 0);
Rectangle2D rect = new Rectangle2D.Double(0.75, -2.5, 0.5, 2);
if (disp) {
int width = 1500;
int height = 1000;
Object[] a = getGraphics(width, height);
BufferedImage bi = (BufferedImage) a[0];
Graphics2D gi = (Graphics2D) a[1];
gi.translate(width/2, height/2);
gi.scale(40, 40);
gi.setStroke(new BasicStroke(1f/20));
gi.setColor(Color.blue);
gi.draw(new Line2D.Double(cc.getX1(), cc.getY1(), cc.getX2(), cc.getY2()));
gi.setColor(Color.yellow);
gi.draw(new Line2D.Double(cc.getX1(), cc.getY1(), cc.getCtrlX1(), cc.getCtrlY1()));
gi.setColor(Color.cyan);
gi.draw(new Line2D.Double(cc.getCtrlX2(), cc.getCtrlY2(), cc.getX2(), cc.getY2()));
gi.setColor(Color.green);
gi.draw(cc);
gi.setColor(Color.red);
gi.draw(rect);
display(bi);
}
return cc.contains(rect);
}
private static Object[] getGraphics(int w, int h) {
BufferedImage bi = new BufferedImage(w, h, BufferedImage.TYPE_INT_RGB);
Graphics2D g2 = (Graphics2D) bi.getGraphics();
g2.setColor(Color.white);
g2.fillRect(0, 0, w, h) ;
g2.setColor(Color.black);
return new Object[] {bi, g2};
}
private static void display(final BufferedImage im) {
SwingUtilities.invokeLater(new Runnable() {
public void run() {
JFrame f = new JFrame();
f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
f.add(new Panel() {
public void paint(Graphics g) {
g.drawImage(im, 0, 0, null);
}
});
f.pack();
f.setVisible(true);
f.setSize(im.getWidth(), im.getHeight());
}
});
}
}
