Review URL: https://codereview.appspot.com/5576043

git-svn-id: http://skia.googlecode.com/svn/trunk@3087 2bbb7eff-a529-9590-31e7-b0007b416f81
This commit is contained in:
caryclark@google.com 2012-01-25 18:57:23 +00:00
Родитель 76bd2540b5
Коммит 27accef223
46 изменённых файлов: 2308 добавлений и 317 удалений

Просмотреть файл

@ -1,5 +1,5 @@
#include "CubicIntersection.h"
#include "CubicIntersection_Tests.h"
#include "Intersection_Tests.h"
#include "IntersectionUtilities.h"
const Cubic convex[] = {

Просмотреть файл

@ -2,7 +2,7 @@
//
#include <math.h>
#include "CubicIntersection.h"
#include "CubicUtilities.h"
#define TEST_ALTERNATIVES 0
#if TEST_ALTERNATIVES

Просмотреть файл

@ -1,8 +1,8 @@
#include "CubicIntersection.h"
#include "CubicIntersection_TestData.h"
#include "CubicIntersection_Tests.h"
#include "Intersection_Tests.h"
void BezierClip_Test() {
void CubicBezierClip_Test() {
for (size_t index = 0; index < tests_count; ++index) {
const Cubic& cubic1 = tests[index][0];
const Cubic& cubic2 = tests[index][1];

Просмотреть файл

@ -1,45 +0,0 @@
/*
Suppose two cubics are coincident. Then a third cubic exists described by two
of the four endpoints. The coincident endpoints must be on or inside the convex
hulls of both cubics.
The coincident control points, while unknown, must lie on the line segment from
the coincident end point to its original control point.
Given a cubic c1, c2, A, and D:
A = c1[0]*(1 - t1)*(1 - t1)*(1 - t1) + 3*c1[1]*t1*(1 - t1)*(1 - t1) + 3*c1[2]*t1*t1*(1 - t1) + c1[3]*t1*t1*t1
D = c2[0]*(1 - t2)*(1 - t2)*(1 - t2) + 3*c2[1]*t2*(1 - t2)*(1 - t2) + 3*c2[2]*t2*t2*(1 - t2) + c213]*t2*t2*t2
Assuming that A was originally from c2:
B = c2[0]*(1 - t2) + c2[1]*t2
C = c1[0]*(1 - t1) + c1[0]*t1
If both ends of the same cubic is contained in the convex hull of the other,
then, are there t values of the larger cubic that describe the end points; is
the mid-value of those ts on the smaller cubic, and, are the tangents at the
smaller ends the same as on both.
This all requires knowing the t values.
Maybe solving the cubic is possible. Given x, find t. Given t, find y.
see en.wikipedia.org/wiki/Cubic_polynomial
Another way of writing the solution may be obtained by noting that the proof of above formula shows that the product of the two cube roots is rational. This gives the following formula in which or stands for any choice of the square or cube root, if
If and , the sign of has to be chosen to have .
If and , the three roots are equal:
If Q = 0 and , the above expression for the roots is correct but misleading, hiding the fact that no radical is needed to represent the roots. In fact, in this case, there is a double root,
and a simple root
*/

Просмотреть файл

@ -9,23 +9,27 @@ void chop_at(const Cubic& src, CubicPair& dst, double t);
void chop_at(const Quadratic& src, QuadraticPair& dst, double t);
int convex_hull(const Cubic& cubic, char order[4]);
bool convex_x_hull(const Cubic& cubic, char connectTo0[2], char connectTo3[2]);
double cube_root(double x);
int cubic_roots(const double coeff[4], double tValues[3]);
bool implicit_matches(const Cubic& cubic1, const Cubic& cubic2);
bool implicit_matches(const Quadratic& quad1, const Quadratic& quad2);
int quadratic_roots(const double coeff[3], double tValues[2]);
void sub_divide(const Cubic& src, double t1, double t2, Cubic& dst);
void sub_divide(const Quadratic& src, double t1, double t2, Quadratic& dst);
void tangent(const Cubic& cubic, double t, _Point& result);
void tangent(const Quadratic& cubic, double t, _Point& result);
// utilities used only for unit testing
bool point_on_parameterized_curve(const Cubic& cubic, const _Point& point);
bool point_on_parameterized_curve(const Quadratic& quad, const _Point& point);
// main functions
enum ReduceOrder_Flags {
kReduceOrder_NoQuadraticsAllowed,
kReduceOrder_QuadraticsAllowed
};
int reduceOrder(const Cubic& cubic, Cubic& reduction, ReduceOrder_Flags );
int reduceOrder(const _Line& line, _Line& reduction);
int reduceOrder(const Quadratic& quad, Quadratic& reduction);
bool intersectStart(const Cubic& cubic1, const Cubic& cubic2, Intersections& );
//bool intersectStart(const Cubic& cubic1, const Cubic& cubic2, Intersections& );
bool intersectStartT(const Cubic& cubic1, const Cubic& cubic2, Intersections& );
bool intersectStart(const Cubic& cubic, const _Line& line, Intersections& );
bool intersectStart(const Quadratic& q1, const Quadratic& q2, Intersections& );
bool intersectStart(const Quadratic& quad, const _Line& line, Intersections& );

Просмотреть файл

@ -44,9 +44,11 @@ bool intersect(double minT1, double maxT1, double minT2, double maxT2) {
sub_divide(cubic1, minT1, maxT1, intersections.swapped() ? larger : smaller);
sub_divide(cubic2, minT2, maxT2, intersections.swapped() ? smaller : larger);
Cubic smallResult;
if (reduceOrder(smaller, smallResult, kReduceOrder_NoQuadraticsAllowed) <= 2) {
if (reduceOrder(smaller, smallResult,
kReduceOrder_NoQuadraticsAllowed) <= 2) {
Cubic largeResult;
if (reduceOrder(larger, largeResult, kReduceOrder_NoQuadraticsAllowed) <= 2) {
if (reduceOrder(larger, largeResult,
kReduceOrder_NoQuadraticsAllowed) <= 2) {
_Point pt;
const _Line& smallLine = (const _Line&) smallResult;
const _Line& largeLine = (const _Line&) largeResult;
@ -87,16 +89,23 @@ bool intersect(double minT1, double maxT1, double minT2, double maxT2) {
double newMinT1 = interp(minT1, maxT1, minT);
double newMaxT1 = interp(minT1, maxT1, maxT);
split = (newMaxT1 - newMinT1 > (maxT1 - minT1) * tClipLimit) << 1;
printf("%s d=%d s=%d new1=(%g,%g) old1=(%g,%g) split=%d\n", __FUNCTION__, depth,
splits, newMinT1, newMaxT1, minT1, maxT1, split);
#define VERBOSE 0
#if VERBOSE
printf("%s d=%d s=%d new1=(%g,%g) old1=(%g,%g) split=%d\n",
__FUNCTION__, depth, splits, newMinT1, newMaxT1, minT1, maxT1,
split);
#endif
minT1 = newMinT1;
maxT1 = newMaxT1;
} else {
double newMinT2 = interp(minT2, maxT2, minT);
double newMaxT2 = interp(minT2, maxT2, maxT);
split = newMaxT2 - newMinT2 > (maxT2 - minT2) * tClipLimit;
printf("%s d=%d s=%d new2=(%g,%g) old2=(%g,%g) split=%d\n", __FUNCTION__, depth,
splits, newMinT2, newMaxT2, minT2, maxT2, split);
#if VERBOSE
printf("%s d=%d s=%d new2=(%g,%g) old2=(%g,%g) split=%d\n",
__FUNCTION__, depth, splits, newMinT2, newMaxT2, minT2, maxT2,
split);
#endif
minT2 = newMinT2;
maxT2 = newMaxT2;
}

Просмотреть файл

@ -1,6 +1,6 @@
#include "CubicIntersection.h"
#include "CubicIntersection_TestData.h"
#include "CubicIntersection_Tests.h"
#include "Intersection_Tests.h"
#include "Intersections.h"
#include "TestUtilities.h"
@ -14,51 +14,37 @@ void CubicIntersection_Test() {
int order1 = reduceOrder(cubic1, reduce1, kReduceOrder_NoQuadraticsAllowed);
int order2 = reduceOrder(cubic2, reduce2, kReduceOrder_NoQuadraticsAllowed);
if (order1 < 4) {
printf("[%d] cubic1 order=%d\n", (int) index, order1);
printf("%s [%d] cubic1 order=%d\n", __FUNCTION__, (int) index, order1);
continue;
}
if (order2 < 4) {
printf("[%d] cubic2 order=%d\n", (int) index, order2);
printf("%s [%d] cubic2 order=%d\n", __FUNCTION__, (int) index, order2);
continue;
}
if (order1 == 4 && order2 == 4) {
Intersections tIntersections;
intersectStartT(reduce1, reduce2, tIntersections);
#ifdef COMPUTE_BY_CUBIC_SUBDIVISION
Intersections intersections;
intersectStart(reduce1, reduce2, intersections);
#endif
if (tIntersections.intersected()) {
for (int pt = 0; pt < tIntersections.used(); ++pt) {
#ifdef COMPUTE_BY_CUBIC_SUBDIVISION
double t1 = intersections.fT[0][pt];
double x1, y1;
xy_at_t(cubic1, t1, x1, y1);
double t2 = intersections.fT[1][pt];
double x2, y2;
xy_at_t(cubic2, t2, x2, y2);
if (!approximately_equal(x1, x2)) {
printf("%s [%d] (1) t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, pt, t1, x1, y1, t2, x2, y2);
}
if (!approximately_equal(y1, y2)) {
printf("%s [%d] (2) t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, pt, t1, x1, y1, t2, x2, y2);
}
#endif
double tt1 = tIntersections.fT[0][pt];
double tx1, ty1;
xy_at_t(cubic1, tt1, tx1, ty1);
double tt2 = tIntersections.fT[1][pt];
double tx2, ty2;
xy_at_t(cubic2, tt2, tx2, ty2);
if (!approximately_equal(tx1, tx2)) {
printf("%s [%d,%d] x!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
}
if (!approximately_equal(ty1, ty2)) {
printf("%s [%d,%d] y!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
}
}
if (implicit_matches(reduce1, reduce2)) {
printf("%s [%d] coincident\n", __FUNCTION__, (int) index);
continue;
}
Intersections tIntersections;
intersectStartT(reduce1, reduce2, tIntersections);
if (!tIntersections.intersected()) {
printf("%s [%d] no intersection\n", __FUNCTION__, (int) index);
continue;
}
for (int pt = 0; pt < tIntersections.used(); ++pt) {
double tt1 = tIntersections.fT[0][pt];
double tx1, ty1;
xy_at_t(cubic1, tt1, tx1, ty1);
double tt2 = tIntersections.fT[1][pt];
double tx2, ty2;
xy_at_t(cubic2, tt2, tx2, ty2);
if (!approximately_equal(tx1, tx2)) {
printf("%s [%d,%d] x!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
}
if (!approximately_equal(ty1, ty2)) {
printf("%s [%d,%d] y!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
}
}
}

Просмотреть файл

@ -167,9 +167,6 @@ const size_t notLines_count = sizeof(notLines) / sizeof(notLines[0]);
static const double E = PointEpsilon * 2;
static const double F = PointEpsilon * 3;
static const double H = PointEpsilon * 4;
static const double J = PointEpsilon * 5;
static const double K = PointEpsilon * 8; // INVESTIGATE: why are larger multiples necessary?
const Cubic modEpsilonLines[] = {
{{0, E}, {0, 0}, {0, 0}, {1, 0}}, // horizontal
@ -305,69 +302,3 @@ const Cubic negEpsilonLines[] = {
};
const size_t negEpsilonLines_count = sizeof(negEpsilonLines) / sizeof(negEpsilonLines[0]);
const Quadratic quadraticLines[] = {
{{0, 0}, {0, 0}, {1, 0}},
{{0, 0}, {1, 0}, {0, 0}},
{{1, 0}, {0, 0}, {0, 0}},
{{1, 0}, {2, 0}, {3, 0}},
{{0, 0}, {0, 0}, {0, 1}},
{{0, 0}, {0, 1}, {0, 0}},
{{0, 1}, {0, 0}, {0, 0}},
{{0, 1}, {0, 2}, {0, 3}},
{{0, 0}, {0, 0}, {1, 1}},
{{0, 0}, {1, 1}, {0, 0}},
{{1, 1}, {0, 0}, {0, 0}},
{{1, 1}, {2, 2}, {3, 3}},
{{1, 1}, {3, 3}, {3, 3}},
{{1, 1}, {1, 1}, {2, 2}},
{{1, 1}, {2, 2}, {1, 1}},
{{1, 1}, {1, 1}, {3, 3}},
{{1, 1}, {2, 2}, {4, 4}}, // no coincident
{{1, 1}, {3, 3}, {4, 4}},
{{1, 1}, {3, 3}, {2, 2}},
{{1, 1}, {4, 4}, {2, 2}},
{{1, 1}, {4, 4}, {3, 3}},
{{2, 2}, {1, 1}, {3, 3}},
{{2, 2}, {1, 1}, {4, 4}},
{{2, 2}, {3, 3}, {1, 1}},
{{2, 2}, {3, 3}, {4, 4}},
{{2, 2}, {4, 4}, {1, 1}},
{{2, 2}, {4, 4}, {3, 3}},
};
const size_t quadraticLines_count = sizeof(quadraticLines) / sizeof(quadraticLines[0]);
const Quadratic quadraticModEpsilonLines[] = {
{{0, F}, {0, 0}, {1, 0}},
{{0, 0}, {1, 0}, {0, F}},
{{1, 0}, {0, F}, {0, 0}},
{{1, H}, {2, 0}, {3, 0}},
{{F, 0}, {0, 0}, {0, 1}},
{{0, 0}, {0, 1}, {F, 0}},
{{0, 1}, {F, 0}, {0, 0}},
{{H, 1}, {0, 2}, {0, 3}},
{{0, F}, {0, 0}, {1, 1}},
{{0, 0}, {1, 1}, {F, 0}},
{{1, 1}, {F, 0}, {0, 0}},
{{1, 1+J}, {2, 2}, {3, 3}},
{{1, 1}, {3, 3}, {3+F, 3}},
{{1, 1}, {1+F, 1}, {2, 2}},
{{1, 1}, {2, 2}, {1, 1+F}},
{{1, 1}, {1, 1+F}, {3, 3}},
{{1+H, 1}, {2, 2}, {4, 4}}, // no coincident
{{1, 1+K}, {3, 3}, {4, 4}},
{{1, 1}, {3+F, 3}, {2, 2}},
{{1, 1}, {4, 4+F}, {2, 2}},
{{1, 1}, {4, 4}, {3+F, 3}},
{{2, 2}, {1, 1}, {3, 3+F}},
{{2+F, 2}, {1, 1}, {4, 4}},
{{2, 2+F}, {3, 3}, {1, 1}},
{{2, 2}, {3+F, 3}, {4, 4}},
{{2, 2}, {4, 4+F}, {1, 1}},
{{2, 2}, {4, 4}, {3+F, 3}},
};
const size_t quadraticModEpsilonLines_count = sizeof(quadraticModEpsilonLines) / sizeof(quadraticModEpsilonLines[0]);

Просмотреть файл

@ -12,8 +12,6 @@ extern const Cubic notLines[];
extern const Cubic modEpsilonLines[];
extern const Cubic lessEpsilonLines[];
extern const Cubic negEpsilonLines[];
extern const Quadratic quadraticLines[];
extern const Quadratic quadraticModEpsilonLines[];
extern const size_t pointDegenerates_count;
extern const size_t notPointDegenerates_count;
@ -24,5 +22,3 @@ extern const size_t notLines_count;
extern const size_t modEpsilonLines_count;
extern const size_t lessEpsilonLines_count;
extern const size_t negEpsilonLines_count;
extern const size_t quadraticLines_count;
extern const size_t quadraticModEpsilonLines_count;

Просмотреть файл

@ -1,15 +0,0 @@
#include "CubicIntersection_Tests.h"
void CubicIntersection_Tests() {
// tests are in dependency order
Inline_Tests();
ConvexHull_Test();
ConvexHull_X_Test();
LineParameter_Test();
LineIntersection_Test();
QuadraticCoincidence_Test();
CubicCoincidence_Test();
ReduceOrder_Test();
// BezierClip_Test();
CubicIntersection_Test();
}

Просмотреть файл

@ -1,11 +0,0 @@
void BezierClip_Test();
void ConvexHull_Test();
void ConvexHull_X_Test();
void CubicCoincidence_Test();
void CubicIntersection_Test();
void CubicIntersection_Tests();
void Inline_Tests();
void LineIntersection_Test();
void LineParameter_Test();
void QuadraticCoincidence_Test();
void ReduceOrder_Test();

Просмотреть файл

@ -226,7 +226,7 @@ static double (*calc_proc[])(double a, double b, double c, double d,
/* Control points to parametric coefficients
s = 1 - t
Attt + 3Btt2 + 3Ctss + Dsss ==
Attt + 3Btts + 3Ctss + Dsss ==
Attt + 3B(1 - t)tt + 3C(1 - t)(t - tt) + D(1 - t)(1 - 2t + tt) ==
Attt + 3B(tt - ttt) + 3C(t - tt - tt + ttt) + D(1-2t+tt-t+2tt-ttt) ==
Attt + 3Btt - 3Bttt + 3Ct - 6Ctt + 3Cttt + D - 3Dt + 3Dtt - Dttt ==
@ -236,6 +236,23 @@ static double (*calc_proc[])(double a, double b, double c, double d,
c = 3*C - 3*D
d = D
*/
/* http://www.algorithmist.net/bezier3.html
p = 3 * A
q = 3 * B
r = 3 * C
a = A
b = q - p
c = p - 2 * q + r
d = D - A + q - r
B(t) = a + t * (b + t * (c + t * d))
so
B(t) = a + t*b + t*t*(c + t*d)
= a + t*b + t*t*c + t*t*t*d
*/
static void set_abcd(const double* cubic, double& a, double& b, double& c,
double& d) {
a = cubic[0]; // a = A
@ -251,22 +268,47 @@ static void calc_bc(const double d, double& b, double& c) {
b -= c; // b = 3*B - 6*C + 3*D
}
static void alt_set_abcd(const double* cubic, double& a, double& b, double& c,
double& d) {
a = cubic[0];
double p = 3 * a;
double q = 3 * cubic[2];
double r = 3 * cubic[4];
b = q - p;
c = p - 2 * q + r;
d = cubic[6] - a + q - r;
}
const bool try_alt = true;
bool implicit_matches(const Cubic& one, const Cubic& two) {
double p1[coeff_count]; // a'xxx , b'xxy , c'xyy , d'xx , e'xy , f'yy, etc.
double p2[coeff_count];
double a1, b1, c1, d1;
set_abcd(&one[0].x, a1, b1, c1, d1);
if (try_alt)
alt_set_abcd(&one[0].x, a1, b1, c1, d1);
else
set_abcd(&one[0].x, a1, b1, c1, d1);
double e1, f1, g1, h1;
set_abcd(&one[0].y, e1, f1, g1, h1);
if (try_alt)
alt_set_abcd(&one[0].y, e1, f1, g1, h1);
else
set_abcd(&one[0].y, e1, f1, g1, h1);
calc_ABCD(a1, e1, p1);
double a2, b2, c2, d2;
set_abcd(&two[0].x, a2, b2, c2, d2);
if (try_alt)
alt_set_abcd(&two[0].x, a2, b2, c2, d2);
else
set_abcd(&two[0].x, a2, b2, c2, d2);
double e2, f2, g2, h2;
set_abcd(&two[0].y, e2, f2, g2, h2);
if (try_alt)
alt_set_abcd(&two[0].y, e2, f2, g2, h2);
else
set_abcd(&two[0].y, e2, f2, g2, h2);
calc_ABCD(a2, e2, p2);
int first = 0;
for (int index = 0; index < coeff_count; ++index) {
if (index == xx_coeff) {
if (!try_alt && index == xx_coeff) {
calc_bc(d1, b1, c1);
calc_bc(h1, f1, g1);
calc_bc(d2, b2, c2);
@ -304,3 +346,7 @@ void tangent(const Cubic& cubic, double t, _Point& result) {
result.y = tangent(&cubic[0].y, t);
}
// unit test to return and validate parametric coefficients
#include "CubicParameterization_TestUtility.cpp"

Просмотреть файл

@ -4,7 +4,7 @@
* Resultant[a*t^3 + b*t^2 + c*t + d - x, e*t^3 + f*t^2 + g*t + h - y, t]
*/
const char result[] =
const char result1[] =
"-d^3 e^3 + c d^2 e^2 f - b d^2 e f^2 + a d^2 f^3 - c^2 d e^2 g + "
" 2 b d^2 e^2 g + b c d e f g - 3 a d^2 e f g - a c d f^2 g - "
" b^2 d e g^2 + 2 a c d e g^2 + a b d f g^2 - a^2 d g^3 + c^3 e^2 h - "
@ -31,7 +31,48 @@ const char result[] =
" 3 a^2 d e y^2 + a b^2 f y^2 - 2 a^2 c f y^2 - a^2 b g y^2 + "
" 3 a^3 h y^2 + 3 a^2 e x y^2 - a^3 y^3";
const size_t len = sizeof(result) - 1;
const size_t len1 = sizeof(result1) - 1;
/* Given:
* Expand[
* Det[{{a, b, c, (d - x), 0, 0},
* {0, a, b, c, (d - x), 0},
* {0, 0, a, b, c, (d - x)},
* {e, f, g, (h - y), 0, 0},
* {0, e, f, g, (h - y), 0},
* {0, 0, e, f, g, (h - y)}}]]
*/
// result1 and result2 are the same. 102 factors:
const char result2[] =
"-d^3 e^3 + c d^2 e^2 f - b d^2 e f^2 + a d^2 f^3 - c^2 d e^2 g + "
" 2 b d^2 e^2 g + b c d e f g - 3 a d^2 e f g - a c d f^2 g - "
" b^2 d e g^2 + 2 a c d e g^2 + a b d f g^2 - a^2 d g^3 + c^3 e^2 h - "
" 3 b c d e^2 h + 3 a d^2 e^2 h - b c^2 e f h + 2 b^2 d e f h + "
" a c d e f h + a c^2 f^2 h - 2 a b d f^2 h + b^2 c e g h - "
" 2 a c^2 e g h - a b d e g h - a b c f g h + 3 a^2 d f g h + "
" a^2 c g^2 h - b^3 e h^2 + 3 a b c e h^2 - 3 a^2 d e h^2 + "
" a b^2 f h^2 - 2 a^2 c f h^2 - a^2 b g h^2 + a^3 h^3 + 3 d^2 e^3 x - "
" 2 c d e^2 f x + 2 b d e f^2 x - 2 a d f^3 x + c^2 e^2 g x - "
" 4 b d e^2 g x - b c e f g x + 6 a d e f g x + a c f^2 g x + "
" b^2 e g^2 x - 2 a c e g^2 x - a b f g^2 x + a^2 g^3 x + "
" 3 b c e^2 h x - 6 a d e^2 h x - 2 b^2 e f h x - a c e f h x + "
" 2 a b f^2 h x + a b e g h x - 3 a^2 f g h x + 3 a^2 e h^2 x - "
" 3 d e^3 x^2 + c e^2 f x^2 - b e f^2 x^2 + a f^3 x^2 + "
" 2 b e^2 g x^2 - 3 a e f g x^2 + 3 a e^2 h x^2 + e^3 x^3 - "
" c^3 e^2 y + 3 b c d e^2 y - 3 a d^2 e^2 y + b c^2 e f y - "
" 2 b^2 d e f y - a c d e f y - a c^2 f^2 y + 2 a b d f^2 y - "
" b^2 c e g y + 2 a c^2 e g y + a b d e g y + a b c f g y - "
" 3 a^2 d f g y - a^2 c g^2 y + 2 b^3 e h y - 6 a b c e h y + "
" 6 a^2 d e h y - 2 a b^2 f h y + 4 a^2 c f h y + 2 a^2 b g h y - "
" 3 a^3 h^2 y - 3 b c e^2 x y + 6 a d e^2 x y + 2 b^2 e f x y + "
" a c e f x y - 2 a b f^2 x y - a b e g x y + 3 a^2 f g x y - "
" 6 a^2 e h x y - 3 a e^2 x^2 y - b^3 e y^2 + 3 a b c e y^2 - "
" 3 a^2 d e y^2 + a b^2 f y^2 - 2 a^2 c f y^2 - a^2 b g y^2 + "
" 3 a^3 h y^2 + 3 a^2 e x y^2 - a^3 y^3";
const size_t len2 = sizeof(result2) - 1;
const int factors = 8;
struct coeff {
@ -56,91 +97,93 @@ enum {
typedef std::vector<coeff> coeffs;
typedef std::vector<coeffs> n_coeffs;
static char skipSpace(size_t& index) {
static char skipSpace(const char* str, size_t& index) {
do {
++index;
} while (result[index] == ' ');
return result[index];
} while (str[index] == ' ');
return str[index];
}
static char backSkipSpace(size_t& end) {
while (result[end - 1] == ' ') {
static char backSkipSpace(const char* str, size_t& end) {
while (str[end - 1] == ' ') {
--end;
}
return result[end - 1];
return str[end - 1];
}
static void match(coeffs& co, const char pattern[]) {
static void match(const char* str, size_t len, coeffs& co, const char pattern[]) {
size_t patternLen = strlen(pattern);
size_t index = 0;
while (index < len) {
char ch = result[index];
char ch = str[index];
if (ch != '-' && ch != '+') {
printf("missing sign\n");
}
size_t end = index + 1;
while (result[end] != '+' && result[end] != '-' && ++end < len) {
while (str[end] != '+' && str[end] != '-' && ++end < len) {
;
}
backSkipSpace(end);
backSkipSpace(str, end);
size_t idx = index;
index = end;
skipSpace(index);
if (!strncmp(&result[end - patternLen], pattern, patternLen) == 0) {
skipSpace(str, index);
if (!strncmp(&str[end - patternLen], pattern, patternLen) == 0) {
continue;
}
size_t endCoeff = end - patternLen;
char last = backSkipSpace(endCoeff);
char last = backSkipSpace(str, endCoeff);
if (last == '2' || last == '3') {
last = result[endCoeff - 3]; // skip ^2
last = str[endCoeff - 3]; // skip ^2
}
if (last == 'x' || last == 'y') {
continue;
}
coeff c;
c.s = result[idx] == '-' ? -1 : 1;
c.s = str[idx] == '-' ? -1 : 1;
bzero(c.n, sizeof(c.n));
ch = skipSpace(idx);
ch = skipSpace(str, idx);
if (ch >= '2' && ch <= '6') {
c.s *= ch - '0';
ch = skipSpace(idx);
ch = skipSpace(str, idx);
}
while (idx < endCoeff) {
char x = result[idx];
char x = str[idx];
if (x < 'a' || x > 'a' + factors) {
printf("expected factor\n");
}
idx++;
int pow = 1;
if (result[idx] == '^') {
if (str[idx] == '^') {
idx++;
char exp = result[idx];
char exp = str[idx];
if (exp < '2' || exp > '3') {
printf("expected exponent\n");
}
pow = exp - '0';
}
skipSpace(idx);
skipSpace(str, idx);
c.n[x - 'a'] = pow;
}
co.push_back(c);
}
}
void cubecode_test();
void cubecode_test(int test);
void cubecode_test() {
void cubecode_test(int test) {
const char* str = test ? result2 : result1;
size_t len = strlen(str);
n_coeffs c(coeff_count);
match(c[xxx_coeff], "x^3"); // 1 factor
match(c[xxy_coeff], "x^2 y"); // 1 factor
match(c[xyy_coeff], "x y^2"); // 1 factor
match(c[yyy_coeff], "y^3"); // 1 factor
match(c[xx_coeff], "x^2"); // 7 factors
match(c[xy_coeff], "x y"); // 8 factors
match(c[yy_coeff], "y^2"); // 7 factors
match(c[x_coeff], "x"); // 21 factors
match(c[y_coeff], "y"); // 21 factors
match(c[c_coeff], ""); // 34 factors
match(str, len, c[xxx_coeff], "x^3"); // 1 factor
match(str, len, c[xxy_coeff], "x^2 y"); // 1 factor
match(str, len, c[xyy_coeff], "x y^2"); // 1 factor
match(str, len, c[yyy_coeff], "y^3"); // 1 factor
match(str, len, c[xx_coeff], "x^2"); // 7 factors
match(str, len, c[xy_coeff], "x y"); // 8 factors
match(str, len, c[yy_coeff], "y^2"); // 7 factors
match(str, len, c[x_coeff], "x"); // 21 factors
match(str, len, c[y_coeff], "y"); // 21 factors
match(str, len, c[c_coeff], ""); // 34 factors : total 102
#define COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS 0
#define WRITE_AS_NONOPTIMIZED_C_CODE 0
#if COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS

Просмотреть файл

@ -1,5 +1,5 @@
#include "CubicIntersection.h"
#include "CubicIntersection_Tests.h"
#include "Intersection_Tests.h"
#include "TestUtilities.h"
const Quadratic quadratics[] = {
@ -42,3 +42,61 @@ void CubicCoincidence_Test() {
}
}
}
// pairs of coincident cubics
// The on curve points of each cubic should be on both parameterized cubics.
const Cubic cubics[] = {
{
{1, -1},
{.333, 1},
{-.333, -1},
{-1, 1}
},
{
{-1, 1},
{-.333, -1},
{.333, 1},
{1, -1}
},
{
{0, 2},
{0, 1},
{1, 0},
{2, 0}
}, {
{2, 0},
{1, 0},
{0, 1},
{0, 2}
},
{
{315.74799999999999, 312.83999999999997},
{312.64400000000001, 318.13400000000001},
{305.83600000000001, 319.90899999999999},
{300.54199999999997, 316.80399999999997}
}, {
{317.12200000000001, 309.05000000000001},
{316.11200000000002, 315.10199999999998},
{310.38499999999999, 319.19},
{304.33199999999999, 318.17899999999997}
}
};
const size_t cubics_count = sizeof(cubics) / sizeof(cubics[0]);
int firstCubicParameterizationTest = 0;
void CubicParameterization_Test() {
for (size_t index = firstCubicParameterizationTest; index < cubics_count; ++index) {
for (size_t inner = 0; inner < 4; inner += 3) {
if (!point_on_parameterized_curve(cubics[index], cubics[index][inner])) {
printf("%s [%zu,%zu] 1 parameterization failed\n",
__FUNCTION__, index, inner);
}
if (!point_on_parameterized_curve(cubics[index], cubics[index ^ 1][inner])) {
printf("%s [%zu,%zu] 2 parameterization failed\n",
__FUNCTION__, index, inner);
}
}
}
}

Просмотреть файл

@ -0,0 +1,39 @@
// included by CubicParameterization.cpp
// accesses internal functions to validate parameterized coefficients
static void parameter_coeffs(const Cubic& cubic, double coeffs[coeff_count]) {
double ax, bx, cx, dx;
if (try_alt)
alt_set_abcd(&cubic[0].x, ax, bx, cx, dx);
else
set_abcd(&cubic[0].x, ax, bx, cx, dx);
double ay, by, cy, dy;
if (try_alt)
alt_set_abcd(&cubic[0].y, ay, by, cy, dy);
else
set_abcd(&cubic[0].y, ay, by, cy, dy);
calc_ABCD(ax, ay, coeffs);
if (!try_alt) calc_bc(dx, bx, cx);
if (!try_alt) calc_bc(dy, by, cy);
for (int index = xx_coeff; index < coeff_count; ++index) {
int procIndex = index - xx_coeff;
coeffs[index] = (*calc_proc[procIndex])(ax, bx, cx, dx, ay, by, cy, dy);
}
}
bool point_on_parameterized_curve(const Cubic& cubic, const _Point& point) {
double coeffs[coeff_count];
parameter_coeffs(cubic, coeffs);
double xxx = coeffs[xxx_coeff] * point.x * point.x * point.x;
double xxy = coeffs[xxy_coeff] * point.x * point.x * point.y;
double xyy = coeffs[xyy_coeff] * point.x * point.y * point.y;
double yyy = coeffs[yyy_coeff] * point.y * point.y * point.y;
double xx = coeffs[ xx_coeff] * point.x * point.x;
double xy = coeffs[ xy_coeff] * point.x * point.y;
double yy = coeffs[ yy_coeff] * point.y * point.y;
double x = coeffs[ x_coeff] * point.x;
double y = coeffs[ y_coeff] * point.y;
double c = coeffs[ c_coeff];
double sum = xxx + xxy + xyy + yyy + xx + xy + yy + x + y + c;
return approximately_zero(sum);
}

Просмотреть файл

@ -0,0 +1,138 @@
#include "CubicIntersection.h"
#include "CubicIntersection_TestData.h"
#include "Intersection_Tests.h"
#include "QuadraticIntersection_TestData.h"
#include "TestUtilities.h"
void CubicReduceOrder_Test() {
size_t index;
Cubic reduce;
int order;
enum {
RunAll,
RunPointDegenerates,
RunNotPointDegenerates,
RunLines,
RunNotLines,
RunModEpsilonLines,
RunLessEpsilonLines,
RunNegEpsilonLines,
RunQuadraticLines,
RunQuadraticModLines,
RunComputedLines,
RunNone
} run = RunAll;
int firstTestIndex = 0;
#if 0
run = RunComputedLines;
firstTestIndex = 18;
#endif
int firstPointDegeneratesTest = run == RunAll ? 0 : run == RunPointDegenerates ? firstTestIndex : INT_MAX;
int firstNotPointDegeneratesTest = run == RunAll ? 0 : run == RunNotPointDegenerates ? firstTestIndex : INT_MAX;
int firstLinesTest = run == RunAll ? 0 : run == RunLines ? firstTestIndex : INT_MAX;
int firstNotLinesTest = run == RunAll ? 0 : run == RunNotLines ? firstTestIndex : INT_MAX;
int firstModEpsilonTest = run == RunAll ? 0 : run == RunModEpsilonLines ? firstTestIndex : INT_MAX;
int firstLessEpsilonTest = run == RunAll ? 0 : run == RunLessEpsilonLines ? firstTestIndex : INT_MAX;
int firstNegEpsilonTest = run == RunAll ? 0 : run == RunNegEpsilonLines ? firstTestIndex : INT_MAX;
int firstQuadraticLineTest = run == RunAll ? 0 : run == RunQuadraticLines ? firstTestIndex : INT_MAX;
int firstQuadraticModLineTest = run == RunAll ? 0 : run == RunQuadraticModLines ? firstTestIndex : INT_MAX;
int firstComputedLinesTest = run == RunAll ? 0 : run == RunComputedLines ? firstTestIndex : INT_MAX;
for (index = firstPointDegeneratesTest; index < pointDegenerates_count; ++index) {
const Cubic& cubic = pointDegenerates[index];
order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
if (order != 1) {
printf("[%d] pointDegenerates order=%d\n", (int) index, order);
}
}
for (index = firstNotPointDegeneratesTest; index < notPointDegenerates_count; ++index) {
const Cubic& cubic = notPointDegenerates[index];
order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
if (order == 1) {
printf("[%d] notPointDegenerates order=%d\n", (int) index, order);
}
}
for (index = firstLinesTest; index < lines_count; ++index) {
const Cubic& cubic = lines[index];
order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
if (order != 2) {
printf("[%d] lines order=%d\n", (int) index, order);
}
}
for (index = firstNotLinesTest; index < notLines_count; ++index) {
const Cubic& cubic = notLines[index];
order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
if (order == 2) {
printf("[%d] notLines order=%d\n", (int) index, order);
}
}
for (index = firstModEpsilonTest; index < modEpsilonLines_count; ++index) {
const Cubic& cubic = modEpsilonLines[index];
order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
if (order == 2) {
printf("[%d] line mod by epsilon order=%d\n", (int) index, order);
}
}
for (index = firstLessEpsilonTest; index < lessEpsilonLines_count; ++index) {
const Cubic& cubic = lessEpsilonLines[index];
order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
if (order != 2) {
printf("[%d] line less by epsilon/2 order=%d\n", (int) index, order);
}
}
for (index = firstNegEpsilonTest; index < negEpsilonLines_count; ++index) {
const Cubic& cubic = negEpsilonLines[index];
order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
if (order != 2) {
printf("[%d] line neg by epsilon/2 order=%d\n", (int) index, order);
}
}
for (index = firstQuadraticLineTest; index < quadraticLines_count; ++index) {
const Quadratic& quad = quadraticLines[index];
Cubic cubic;
quad_to_cubic(quad, cubic);
order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
if (order != 2) {
printf("[%d] line quad order=%d\n", (int) index, order);
}
}
for (index = firstQuadraticModLineTest; index < quadraticModEpsilonLines_count; ++index) {
const Quadratic& quad = quadraticModEpsilonLines[index];
Cubic cubic;
quad_to_cubic(quad, cubic);
order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
if (order != 3) {
printf("[%d] line mod quad order=%d\n", (int) index, order);
}
}
// test if computed line end points are valid
for (index = firstComputedLinesTest; index < lines_count; ++index) {
const Cubic& cubic = lines[index];
bool controlsInside = controls_inside(cubic);
order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
if (reduce[0].x == reduce[1].x && reduce[0].y == reduce[1].y) {
printf("[%d] line computed ends match order=%d\n", (int) index, order);
}
if (controlsInside) {
if ( reduce[0].x != cubic[0].x && reduce[0].x != cubic[3].x
|| reduce[0].y != cubic[0].y && reduce[0].y != cubic[3].y
|| reduce[1].x != cubic[0].x && reduce[1].x != cubic[3].x
|| reduce[1].y != cubic[0].y && reduce[1].y != cubic[3].y) {
printf("[%d] line computed ends order=%d\n", (int) index, order);
}
} else {
// binary search for extrema, compare against actual results
// while a control point is outside of bounding box formed by end points, split
_Rect bounds = {DBL_MAX, DBL_MAX, -DBL_MAX, -DBL_MAX};
find_tight_bounds(cubic, bounds);
if ( !approximately_equal(reduce[0].x, bounds.left) && !approximately_equal(reduce[0].x, bounds.right)
|| !approximately_equal(reduce[0].y, bounds.top) && !approximately_equal(reduce[0].y, bounds.bottom)
|| !approximately_equal(reduce[1].x, bounds.left) && !approximately_equal(reduce[1].x, bounds.right)
|| !approximately_equal(reduce[1].y, bounds.top) && !approximately_equal(reduce[1].y, bounds.bottom)) {
printf("[%d] line computed tight bounds order=%d\n", (int) index, order);
}
}
}
}

Просмотреть файл

@ -1,40 +1,32 @@
#include "CubicIntersection.h"
//http://planetmath.org/encyclopedia/CubicEquation.html
/* the roots of x^3 + ax^2 + bx + c are
j = -2a^3 + 9ab - 27c
k = sqrt((2a^3 - 9ab + 27c)^2 + 4(-a^2 + 3b)^3)
t1 = -a/3 + cuberoot((j + k) / 54) + cuberoot((j - k) / 54)
t2 = -a/3 - ( 1 + i*cuberoot(3))/2 * cuberoot((j + k) / 54)
+ (-1 + i*cuberoot(3))/2 * cuberoot((j - k) / 54)
t3 = -a/3 + (-1 + i*cuberoot(3))/2 * cuberoot((j + k) / 54)
- ( 1 + i*cuberoot(3))/2 * cuberoot((j - k) / 54)
*/
#include "CubicUtilities.h"
#include "DataTypes.h"
#include "QuadraticUtilities.h"
const double PI = 4 * atan(1);
static bool is_unit_interval(double x) {
return x > 0 && x < 1;
}
const double PI = 4 * atan(1);
// from SkGeometry.cpp
int cubic_roots(const double coeff[4], double tValues[3]) {
if (approximately_zero(coeff[0])) // we're just a quadratic
{
return quadratic_roots(&coeff[1], tValues);
// from SkGeometry.cpp (and Numeric Solutions, 5.6)
int cubicRoots(double A, double B, double C, double D, double t[3]) {
if (approximately_zero(A)) { // we're just a quadratic
return quadraticRoots(B, C, D, t);
}
double a, b, c;
{
double invA = 1 / A;
a = B * invA;
b = C * invA;
c = D * invA;
}
double inva = 1 / coeff[0];
double a = coeff[1] * inva;
double b = coeff[2] * inva;
double c = coeff[3] * inva;
double a2 = a * a;
double Q = (a2 - b * 3) / 9;
double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
double Q3 = Q * Q * Q;
double R2MinusQ3 = R * R - Q3;
double adiv3 = a / 3;
double* roots = tValues;
double* roots = t;
double r;
if (R2MinusQ3 < 0) // we have 3 real roots
@ -68,5 +60,5 @@ int cubic_roots(const double coeff[4], double tValues[3]) {
if (is_unit_interval(r))
*roots++ = r;
}
return (int)(roots - tValues);
return (int)(roots - t);
}

Просмотреть файл

@ -0,0 +1,2 @@
double cube_root(double x);
int cubicRoots(double A, double B, double C, double D, double t[3]);

Просмотреть файл

@ -20,8 +20,7 @@
#include "SkTouchGesture.h"
#include "SkTypeface.h"
#include "CubicIntersection_Tests.h"
#include "CubicIntersection_TestData.h"
#include "Intersection_Tests.h"
extern void CreateSweep(SkBitmap* , float width);
extern void CreateHorz(SkBitmap* );
@ -162,16 +161,11 @@ SkOSWindow* create_sk_window(void* hwnd, int argc, char** argv) {
return new EdgeWindow(hwnd);
}
void cubecode_test();
void application_init() {
unsigned foo = 4;
SkGraphics::Init();
SkEvent::Init();
cubecode_test();
convert_testx();
CubicIntersection_Tests();
Intersection_Tests();
SkAntiEdge_Test();
}

Просмотреть файл

@ -1,5 +1,5 @@
#include "CubicIntersection.h"
#include "CubicIntersection_Tests.h"
#include "Intersection_Tests.h"
#include "IntersectionUtilities.h"
static void assert_that(int x, int y, const char* s) {

Просмотреть файл

@ -0,0 +1,29 @@
#include "CubicIntersection_TestData.h"
#include "Intersection_Tests.h"
void cubecode_test(int test);
void Intersection_Tests() {
cubecode_test(1);
convert_testx();
// tests are in dependency / complexity order
Inline_Tests();
ConvexHull_Test();
ConvexHull_X_Test();
LineParameter_Test();
LineIntersection_Test();
LineQuadraticIntersection_Test();
LineCubicIntersection_Test();
QuadraticCoincidence_Test();
QuadraticReduceOrder_Test();
QuadraticBezierClip_Test();
QuadraticIntersection_Test();
CubicParameterization_Test();
CubicCoincidence_Test();
CubicReduceOrder_Test();
CubicBezierClip_Test();
CubicIntersection_Test();
}

Просмотреть файл

@ -0,0 +1,18 @@
void ConvexHull_Test();
void ConvexHull_X_Test();
void CubicBezierClip_Test();
void CubicCoincidence_Test();
void CubicIntersection_Test();
void CubicParameterization_Test();
void CubicReduceOrder_Test();
void Inline_Tests();
void Intersection_Tests();
void LineCubicIntersection_Test();
void LineIntersection_Test();
void LineParameter_Test();
void LineQuadraticIntersection_Test();
void QuadraticBezierClip_Test();
void QuadraticCoincidence_Test();
void QuadraticIntersection_Test();
void QuadraticParameterization_Test();
void QuadraticReduceOrder_Test();

Просмотреть файл

@ -0,0 +1,139 @@
#include "CubicIntersection.h"
#include "CubicUtilities.h"
#include "Intersections.h"
#include "LineUtilities.h"
/*
Find the interection of a line and cubic by solving for valid t values.
Analogous to line-quadratic intersection, solve line-cubic intersection by
representing the cubic as:
x = a(1-t)^3 + 2b(1-t)^2t + c(1-t)t^2 + dt^3
y = e(1-t)^3 + 2f(1-t)^2t + g(1-t)t^2 + ht^3
and the line as:
y = i*x + j (if the line is more horizontal)
or:
x = i*y + j (if the line is more vertical)
Then using Mathematica, solve for the values of t where the cubic intersects the
line:
(in) Resultant[
a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - x,
e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - i*x - j, x]
(out) -e + j +
3 e t - 3 f t -
3 e t^2 + 6 f t^2 - 3 g t^2 +
e t^3 - 3 f t^3 + 3 g t^3 - h t^3 +
i ( a -
3 a t + 3 b t +
3 a t^2 - 6 b t^2 + 3 c t^2 -
a t^3 + 3 b t^3 - 3 c t^3 + d t^3 )
if i goes to infinity, we can rewrite the line in terms of x. Mathematica:
(in) Resultant[
a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - i*y - j,
e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
(out) a - j -
3 a t + 3 b t +
3 a t^2 - 6 b t^2 + 3 c t^2 -
a t^3 + 3 b t^3 - 3 c t^3 + d t^3 -
i ( e -
3 e t + 3 f t +
3 e t^2 - 6 f t^2 + 3 g t^2 -
e t^3 + 3 f t^3 - 3 g t^3 + h t^3 )
Solving this with Mathematica produces an expression with hundreds of terms;
instead, use Numeric Solutions recipe to solve the cubic.
The near-horizontal case, in terms of: Ax^3 + Bx^2 + Cx + D == 0
A = (-(-e + 3*f - 3*g + h) + i*(-a + 3*b - 3*c + d) )
B = 3*(-( e - 2*f + g ) + i*( a - 2*b + c ) )
C = 3*(-(-e + f ) + i*(-a + b ) )
D = (-( e ) + i*( a ) + j )
The near-vertical case, in terms of: Ax^3 + Bx^2 + Cx + D == 0
A = ( (-a + 3*b - 3*c + d) - i*(-e + 3*f - 3*g + h) )
B = 3*( ( a - 2*b + c ) - i*( e - 2*f + g ) )
C = 3*( (-a + b ) - i*(-e + f ) )
D = ( ( a ) - i*( e ) - j )
*/
class LineCubicIntersections : public Intersections {
public:
LineCubicIntersections(const Cubic& c, const _Line& l, Intersections& i)
: cubic(c)
, line(l)
, intersections(i) {
}
bool intersect() {
double slope;
double axisIntercept;
moreHorizontal = implicitLine(line, slope, axisIntercept);
double A = cubic[3].x; // d
double B = cubic[2].x * 3; // 3*c
double C = cubic[1].x * 3; // 3*b
double D = cubic[0].x; // a
A -= D - C + B; // A = -a + 3*b - 3*c + d
B += 3 * D - 2 * C; // B = 3*a - 6*b + 3*c
C -= 3 * D; // C = -3*a + 3*b
double E = cubic[3].y; // h
double F = cubic[2].y * 3; // 3*g
double G = cubic[1].y * 3; // 3*f
double H = cubic[0].y; // e
E -= H - G + F; // E = -e + 3*f - 3*g + h
F += 3 * H - 2 * G; // F = 3*e - 6*f + 3*g
G -= 3 * H; // G = -3*e + 3*f
if (moreHorizontal) {
A = A * slope - E;
B = B * slope - F;
C = C * slope - G;
D = D * slope - H + axisIntercept;
} else {
A = A - E * slope;
B = B - F * slope;
C = C - G * slope;
D = D - H * slope - axisIntercept;
}
double t[3];
int roots = cubicRoots(A, B, C, D, t);
for (int x = 0; x < roots; ++x) {
intersections.add(t[x], findLineT(t[x]));
}
return roots > 0;
}
protected:
double findLineT(double t) {
const double* cPtr;
const double* lPtr;
if (moreHorizontal) {
cPtr = &cubic[0].x;
lPtr = &line[0].x;
} else {
cPtr = &cubic[0].y;
lPtr = &line[0].y;
}
double s = 1 - t;
double cubicVal = cPtr[0] * s * s * s + 3 * cPtr[2] * s * s * t
+ 3 * cPtr[4] * s * t * t + cPtr[6] * t * t * t;
return (cubicVal - lPtr[0]) / (lPtr[2] - lPtr[0]);
}
private:
const Cubic& cubic;
const _Line& line;
Intersections& intersections;
bool moreHorizontal;
};
bool intersectStart(const Cubic& cubic, const _Line& line, Intersections& i) {
LineCubicIntersections c(cubic, line, i);
return c.intersect();
}

Просмотреть файл

@ -0,0 +1,55 @@
#include "CubicIntersection.h"
#include "Intersection_Tests.h"
#include "Intersections.h"
#include "LineUtilities.h"
#include "TestUtilities.h"
struct lineCubic {
Cubic cubic;
_Line line;
} lineCubicTests[] = {
{{{0, 0}, {0, 1}, {0, 1}, {1, 1}}, {{0, 1}, {1, 0}}}
};
size_t lineCubicTests_count = sizeof(lineCubicTests) / sizeof(lineCubicTests[0]);
const int firstLineCubicIntersectionTest = 0;
void LineCubicIntersection_Test() {
for (size_t index = firstLineCubicIntersectionTest; index < lineCubicTests_count; ++index) {
const Cubic& cubic = lineCubicTests[index].cubic;
const _Line& line = lineCubicTests[index].line;
Cubic reduce1;
_Line reduce2;
int order1 = reduceOrder(cubic, reduce1, kReduceOrder_NoQuadraticsAllowed);
int order2 = reduceOrder(line, reduce2);
if (order1 < 4) {
printf("[%d] cubic order=%d\n", (int) index, order1);
}
if (order2 < 2) {
printf("[%d] line order=%d\n", (int) index, order2);
}
if (order1 == 4 && order2 == 2) {
Intersections intersections;
intersectStart(reduce1, reduce2, intersections);
if (intersections.intersected()) {
for (int pt = 0; pt < intersections.used(); ++pt) {
double tt1 = intersections.fT[0][pt];
double tx1, ty1;
xy_at_t(cubic, tt1, tx1, ty1);
double tt2 = intersections.fT[1][pt];
double tx2, ty2;
xy_at_t(line, tt2, tx2, ty2);
if (!approximately_equal(tx1, tx2)) {
printf("%s [%d,%d] x!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
}
if (!approximately_equal(ty1, ty2)) {
printf("%s [%d,%d] y!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
}
}
}
}
}
}

Просмотреть файл

@ -65,6 +65,8 @@ bool lineIntersect(const _Line& a, const _Line& b, _Point* result) {
return no_intersection(result);
}
/* Are the line coincident? See if they overlap */
// FIXME: allow returning range of coincidence, instead of or in
// addition to midpoint
paramsA.lineEndPoints(a);
double b0dist = paramsA.pointDistance(b[0]);
bool b0on = approximately_zero_squared(b0dist);

Просмотреть файл

@ -1,6 +1,7 @@
#include "CubicIntersection_Tests.h"
#include "LineIntersection.h"
// FIXME: add tests for intersecting, non-intersecting, degenerate, coincident
const _Line tests[][2] = {
{{{166.86950047022856, 112.69654129527828}, {166.86948801592692, 112.69655741235339}},
{{166.86960700313026, 112.6965477747386}, {166.86925794355412, 112.69656471103423}}},
@ -16,6 +17,9 @@ void LineIntersection_Test() {
const _Line& line2 = tests[index][1];
_Point result;
lineIntersect(line1, line2, &result);
printf("%s (%g,%g)\n", __FUNCTION__, result.x, result.y);
// FIXME: validate results
// see if result is between start and end of both lines
// see if result is on both lines
// printf("%s (%g,%g)\n", __FUNCTION__, result.x, result.y);
}
}

Просмотреть файл

@ -0,0 +1,163 @@
#include "CubicIntersection.h"
#include "Intersections.h"
#include "LineUtilities.h"
#include "QuadraticUtilities.h"
/*
Find the interection of a line and quadratic by solving for valid t values.
From http://stackoverflow.com/questions/1853637/how-to-find-the-mathematical-function-defining-a-bezier-curve
"A Bezier curve is a parametric function. A quadratic Bezier curve (i.e. three
control points) can be expressed as: F(t) = A(1 - t)^2 + B(1 - t)t + Ct^2 where
A, B and C are points and t goes from zero to one.
This will give you two equations:
x = a(1 - t)^2 + b(1 - t)t + ct^2
y = d(1 - t)^2 + e(1 - t)t + ft^2
If you add for instance the line equation (y = kx + m) to that, you'll end up
with three equations and three unknowns (x, y and t)."
Similar to above, the quadratic is represented as
x = a(1-t)^2 + 2b(1-t)t + ct^2
y = d(1-t)^2 + 2e(1-t)t + ft^2
and the line as
y = g*x + h
Using Mathematica, solve for the values of t where the quadratic intersects the
line:
(in) t1 = Resultant[a*(1 - t)^2 + 2*b*(1 - t)*t + c*t^2 - x,
d*(1 - t)^2 + 2*e*(1 - t)*t + f*t^2 - g*x - h, x]
(out) -d + h + 2 d t - 2 e t - d t^2 + 2 e t^2 - f t^2 +
g (a - 2 a t + 2 b t + a t^2 - 2 b t^2 + c t^2)
(in) Solve[t1 == 0, t]
(out) {
{t -> (-2 d + 2 e + 2 a g - 2 b g -
Sqrt[(2 d - 2 e - 2 a g + 2 b g)^2 -
4 (-d + 2 e - f + a g - 2 b g + c g) (-d + a g + h)]) /
(2 (-d + 2 e - f + a g - 2 b g + c g))
},
{t -> (-2 d + 2 e + 2 a g - 2 b g +
Sqrt[(2 d - 2 e - 2 a g + 2 b g)^2 -
4 (-d + 2 e - f + a g - 2 b g + c g) (-d + a g + h)]) /
(2 (-d + 2 e - f + a g - 2 b g + c g))
}
}
Numeric Solutions (5.6) suggests to solve the quadratic by computing
Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C))
and using the roots
t1 = Q / A
t2 = C / Q
Using the results above (when the line tends towards horizontal)
A = (-(d - 2*e + f) + g*(a - 2*b + c) )
B = 2*( (d - e ) - g*(a - b ) )
C = (-(d ) + g*(a ) + h )
If g goes to infinity, we can rewrite the line in terms of x.
x = g'*y + h'
And solve accordingly in Mathematica:
(in) t2 = Resultant[a*(1 - t)^2 + 2*b*(1 - t)*t + c*t^2 - g'*y - h',
d*(1 - t)^2 + 2*e*(1 - t)*t + f*t^2 - y, y]
(out) a - h' - 2 a t + 2 b t + a t^2 - 2 b t^2 + c t^2 -
g' (d - 2 d t + 2 e t + d t^2 - 2 e t^2 + f t^2)
(in) Solve[t2 == 0, t]
(out) {
{t -> (2 a - 2 b - 2 d g' + 2 e g' -
Sqrt[(-2 a + 2 b + 2 d g' - 2 e g')^2 -
4 (a - 2 b + c - d g' + 2 e g' - f g') (a - d g' - h')]) /
(2 (a - 2 b + c - d g' + 2 e g' - f g'))
},
{t -> (2 a - 2 b - 2 d g' + 2 e g' +
Sqrt[(-2 a + 2 b + 2 d g' - 2 e g')^2 -
4 (a - 2 b + c - d g' + 2 e g' - f g') (a - d g' - h')])/
(2 (a - 2 b + c - d g' + 2 e g' - f g'))
}
}
Thus, if the slope of the line tends towards vertical, we use:
A = ( (a - 2*b + c) - g'*(d - 2*e + f) )
B = 2*(-(a - b ) + g'*(d - e ) )
C = ( (a ) - g'*(d ) - h' )
*/
class LineQuadraticIntersections : public Intersections {
public:
LineQuadraticIntersections(const Quadratic& q, const _Line& l, Intersections& i)
: quad(q)
, line(l)
, intersections(i) {
}
bool intersect() {
double slope;
double axisIntercept;
moreHorizontal = implicitLine(line, slope, axisIntercept);
double A = quad[2].x; // c
double B = quad[1].x; // b
double C = quad[0].x; // a
A += C - 2 * B; // A = a - 2*b + c
B -= C; // B = -(a - b)
double D = quad[2].y; // f
double E = quad[1].y; // e
double F = quad[0].y; // d
D += F - 2 * E; // D = d - 2*e + f
E -= F; // E = -(d - e)
if (moreHorizontal) {
A = A * slope - D;
B = B * slope - E;
C = C * slope - F + axisIntercept;
} else {
A = A - D * slope;
B = B - E * slope;
C = C - F * slope - axisIntercept;
}
double t[2];
int roots = quadraticRoots(A, B, C, t);
for (int x = 0; x < roots; ++x) {
intersections.add(t[x], findLineT(t[x]));
}
return roots > 0;
}
protected:
double findLineT(double t) {
const double* qPtr;
const double* lPtr;
if (moreHorizontal) {
qPtr = &quad[0].x;
lPtr = &line[0].x;
} else {
qPtr = &quad[0].y;
lPtr = &line[0].y;
}
double s = 1 - t;
double quadVal = qPtr[0] * s * s + 2 * qPtr[2] * s * t + qPtr[4] * t * t;
return (quadVal - lPtr[0]) / (lPtr[2] - lPtr[0]);
}
private:
const Quadratic& quad;
const _Line& line;
Intersections& intersections;
bool moreHorizontal;
};
bool intersectStart(const Quadratic& quad, const _Line& line, Intersections& i) {
LineQuadraticIntersections q(quad, line, i);
return q.intersect();
}

Просмотреть файл

@ -0,0 +1,55 @@
#include "CubicIntersection.h"
#include "Intersection_Tests.h"
#include "Intersections.h"
#include "LineUtilities.h"
#include "TestUtilities.h"
struct lineQuad {
Quadratic quad;
_Line line;
} lineQuadTests[] = {
{{{0, 0}, {0, 1}, {1, 1}}, {{0, 1}, {1, 0}}}
};
size_t lineQuadTests_count = sizeof(lineQuadTests) / sizeof(lineQuadTests[0]);
const int firstLineQuadIntersectionTest = 0;
void LineQuadraticIntersection_Test() {
for (size_t index = firstLineQuadIntersectionTest; index < lineQuadTests_count; ++index) {
const Quadratic& quad = lineQuadTests[index].quad;
const _Line& line = lineQuadTests[index].line;
Quadratic reduce1;
_Line reduce2;
int order1 = reduceOrder(quad, reduce1);
int order2 = reduceOrder(line, reduce2);
if (order1 < 3) {
printf("[%d] quad order=%d\n", (int) index, order1);
}
if (order2 < 2) {
printf("[%d] line order=%d\n", (int) index, order2);
}
if (order1 == 3 && order2 == 2) {
Intersections intersections;
intersectStart(reduce1, reduce2, intersections);
if (intersections.intersected()) {
for (int pt = 0; pt < intersections.used(); ++pt) {
double tt1 = intersections.fT[0][pt];
double tx1, ty1;
xy_at_t(quad, tt1, tx1, ty1);
double tt2 = intersections.fT[1][pt];
double tx2, ty2;
xy_at_t(line, tt2, tx2, ty2);
if (!approximately_equal(tx1, tx2)) {
printf("%s [%d,%d] x!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
}
if (!approximately_equal(ty1, ty2)) {
printf("%s [%d,%d] y!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
}
}
}
}
}
}

Просмотреть файл

@ -0,0 +1,22 @@
#include "LineUtilities.h"
bool implicitLine(const _Line& line, double& slope, double& axisIntercept) {
double lineDx = line[1].x - line[0].x;
double lineDy = line[1].y - line[0].y;
bool moreHorizontal = fabs(lineDx) > fabs(lineDy);
if (moreHorizontal) {
slope = lineDy / lineDx;
axisIntercept = line[0].y - slope * line[0].x;
} else {
slope = lineDx / lineDy;
axisIntercept = line[0].x - slope * line[0].y;
}
return moreHorizontal;
}
int reduceOrder(const _Line& line, _Line& reduced) {
reduced[0] = line[0];
int different = line[0] != line[1];
reduced[1] = line[different];
return 1 + different;
}

Просмотреть файл

@ -0,0 +1,4 @@
#include "DataTypes.h"
bool implicitLine(const _Line& line, double& slope, double& axisIntercept);
int reduceOrder(const _Line& line, _Line& reduced);

Просмотреть файл

@ -48,10 +48,14 @@ bool bezier_clip(const Quadratic& q1, const Quadratic& q2, double& minT, double&
maxT = 1;
}
// Find the intersection of distance convex hull and fat line.
x_at(distance2y[0], distance2y[1], top, bottom, flags, minT, maxT);
x_at(distance2y[1], distance2y[2], top, bottom, flags, minT, maxT);
x_at(distance2y[2], distance2y[0], top, bottom, flags, minT, maxT);
int idx = 0;
do {
int next = idx + 1;
if (next == 3) {
next = 0;
}
x_at(distance2y[idx], distance2y[next], top, bottom, flags, minT, maxT);
idx = next;
} while (idx);
return minT < maxT; // returns false if distance shows no intersection
}

Просмотреть файл

@ -1,10 +1,24 @@
/*
* QuadraticBezierClip_Test.cpp
* edge
*
* Created by Cary Clark on 1/10/12.
* Copyright 2012 __MyCompanyName__. All rights reserved.
*
*/
#include "CubicIntersection.h"
#include "Intersection_Tests.h"
#include "QuadraticIntersection_TestData.h"
void QuadraticBezierClip_Test() {
for (size_t index = 0; index < quadraticTests_count; ++index) {
const Quadratic& quad1 = quadraticTests[index][0];
const Quadratic& quad2 = quadraticTests[index][1];
Quadratic reduce1, reduce2;
int order1 = reduceOrder(quad1, reduce1);
int order2 = reduceOrder(quad2, reduce2);
if (order1 < 3) {
printf("%s [%d] quad1 order=%d\n", __FUNCTION__, (int)index, order1);
}
if (order2 < 3) {
printf("%s [%d] quad2 order=%d\n", __FUNCTION__, (int)index, order2);
}
if (order1 == 3 && order2 == 3) {
double minT = 0;
double maxT = 1;
bezier_clip(reduce1, reduce2, minT, maxT);
}
}
}

Просмотреть файл

@ -185,4 +185,3 @@ principal value is always non-negative.
The other square root is simply 1 times the principal square root; in other
words, the two square roots of a number sum to 0.
*/

Просмотреть файл

@ -0,0 +1,46 @@
#include "CubicIntersection.h"
#include "Intersection_Tests.h"
#include "Intersections.h"
#include "QuadraticIntersection_TestData.h"
#include "TestUtilities.h"
const int firstQuadIntersectionTest = 9;
void QuadraticIntersection_Test() {
for (size_t index = firstQuadIntersectionTest; index < quadraticTests_count; ++index) {
const Quadratic& quad1 = quadraticTests[index][0];
const Quadratic& quad2 = quadraticTests[index][1];
Quadratic reduce1, reduce2;
int order1 = reduceOrder(quad1, reduce1);
int order2 = reduceOrder(quad2, reduce2);
if (order1 < 3) {
printf("[%d] quad1 order=%d\n", (int) index, order1);
}
if (order2 < 3) {
printf("[%d] quad2 order=%d\n", (int) index, order2);
}
if (order1 == 3 && order2 == 3) {
Intersections intersections;
intersectStart(reduce1, reduce2, intersections);
if (intersections.intersected()) {
for (int pt = 0; pt < intersections.used(); ++pt) {
double tt1 = intersections.fT[0][pt];
double tx1, ty1;
xy_at_t(quad1, tt1, tx1, ty1);
double tt2 = intersections.fT[1][pt];
double tx2, ty2;
xy_at_t(quad2, tt2, tx2, ty2);
if (!approximately_equal(tx1, tx2)) {
printf("%s [%d,%d] x!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
}
if (!approximately_equal(ty1, ty2)) {
printf("%s [%d,%d] y!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
}
}
}
}
}
}

Просмотреть файл

@ -9,3 +9,92 @@
#include "QuadraticIntersection_TestData.h"
const Quadratic quadraticLines[] = {
{{0, 0}, {0, 0}, {1, 0}},
{{0, 0}, {1, 0}, {0, 0}},
{{1, 0}, {0, 0}, {0, 0}},
{{1, 0}, {2, 0}, {3, 0}},
{{0, 0}, {0, 0}, {0, 1}},
{{0, 0}, {0, 1}, {0, 0}},
{{0, 1}, {0, 0}, {0, 0}},
{{0, 1}, {0, 2}, {0, 3}},
{{0, 0}, {0, 0}, {1, 1}},
{{0, 0}, {1, 1}, {0, 0}},
{{1, 1}, {0, 0}, {0, 0}},
{{1, 1}, {2, 2}, {3, 3}},
{{1, 1}, {3, 3}, {3, 3}},
{{1, 1}, {1, 1}, {2, 2}},
{{1, 1}, {2, 2}, {1, 1}},
{{1, 1}, {1, 1}, {3, 3}},
{{1, 1}, {2, 2}, {4, 4}}, // no coincident
{{1, 1}, {3, 3}, {4, 4}},
{{1, 1}, {3, 3}, {2, 2}},
{{1, 1}, {4, 4}, {2, 2}},
{{1, 1}, {4, 4}, {3, 3}},
{{2, 2}, {1, 1}, {3, 3}},
{{2, 2}, {1, 1}, {4, 4}},
{{2, 2}, {3, 3}, {1, 1}},
{{2, 2}, {3, 3}, {4, 4}},
{{2, 2}, {4, 4}, {1, 1}},
{{2, 2}, {4, 4}, {3, 3}},
};
const size_t quadraticLines_count = sizeof(quadraticLines) / sizeof(quadraticLines[0]);
static const double F = PointEpsilon * 3;
static const double H = PointEpsilon * 4;
static const double J = PointEpsilon * 5;
static const double K = PointEpsilon * 8; // INVESTIGATE: why are larger multiples necessary?
const Quadratic quadraticModEpsilonLines[] = {
{{0, F}, {0, 0}, {1, 0}},
{{0, 0}, {1, 0}, {0, F}},
{{1, 0}, {0, F}, {0, 0}},
{{1, H}, {2, 0}, {3, 0}},
{{F, 0}, {0, 0}, {0, 1}},
{{0, 0}, {0, 1}, {F, 0}},
{{0, 1}, {F, 0}, {0, 0}},
{{H, 1}, {0, 2}, {0, 3}},
{{0, F}, {0, 0}, {1, 1}},
{{0, 0}, {1, 1}, {F, 0}},
{{1, 1}, {F, 0}, {0, 0}},
{{1, 1+J}, {2, 2}, {3, 3}},
{{1, 1}, {3, 3}, {3+F, 3}},
{{1, 1}, {1+F, 1}, {2, 2}},
{{1, 1}, {2, 2}, {1, 1+F}},
{{1, 1}, {1, 1+F}, {3, 3}},
{{1+H, 1}, {2, 2}, {4, 4}}, // no coincident
{{1, 1+K}, {3, 3}, {4, 4}},
{{1, 1}, {3+F, 3}, {2, 2}},
{{1, 1}, {4, 4+F}, {2, 2}},
{{1, 1}, {4, 4}, {3+F, 3}},
{{2, 2}, {1, 1}, {3, 3+F}},
{{2+F, 2}, {1, 1}, {4, 4}},
{{2, 2+F}, {3, 3}, {1, 1}},
{{2, 2}, {3+F, 3}, {4, 4}},
{{2, 2}, {4, 4+F}, {1, 1}},
{{2, 2}, {4, 4}, {3+F, 3}},
};
const size_t quadraticModEpsilonLines_count = sizeof(quadraticModEpsilonLines) / sizeof(quadraticModEpsilonLines[0]);
const Quadratic quadraticTests[][2] = {
{ // one intersection
{{0, 0},
{0, 1},
{1, 1}},
{{0, 1},
{0, 0},
{1, 0}}
},
{ // four intersections
{{1, 0},
{2, 6},
{3, 0}},
{{0, 1},
{6, 2},
{0, 3}}
}
};
const size_t quadraticTests_count = sizeof(quadraticTests) / sizeof(quadraticTests[0]);

Просмотреть файл

@ -1,9 +1,9 @@
/*
* QuadraticIntersection_TestData.h
* edge
*
* Created by Cary Clark on 1/10/12.
* Copyright 2012 __MyCompanyName__. All rights reserved.
*
*/
#include "DataTypes.h"
extern const Quadratic quadraticLines[];
extern const Quadratic quadraticModEpsilonLines[];
extern const Quadratic quadraticTests[][2];
extern const size_t quadraticLines_count;
extern const size_t quadraticModEpsilonLines_count;
extern const size_t quadraticTests_count;

Просмотреть файл

@ -130,3 +130,5 @@ void tangent(const Quadratic& quadratic, double t, _Point& result) {
result.y = tangent(&quadratic[0].y, t);
}
// unit test to return and validate parametric coefficients
#include "QuadraticParameterization_TestUtility.cpp"

Просмотреть файл

@ -1,5 +1,5 @@
#include "CubicIntersection.h"
#include "CubicIntersection_Tests.h"
#include "Intersection_Tests.h"
const Quadratic quadratics[] = {
{{0, 0}, {1, 0}, {1, 1}},
@ -18,20 +18,23 @@ void QuadraticCoincidence_Test() {
chop_at(test, split, 0.5);
Quadratic midThird;
sub_divide(test, 1.0/3, 2.0/3, midThird);
if (!implicit_matches(test, split.first())) {
printf("%s-1 %d", __FUNCTION__, (int)index);
}
if (!implicit_matches(test, split.second())) {
printf("%s-2 %d", __FUNCTION__, (int)index);
}
if (!implicit_matches(midThird, split.first())) {
printf("%s-3 %d", __FUNCTION__, (int)index);
}
if (!implicit_matches(midThird, split.second())) {
printf("%s-4 %d", __FUNCTION__, (int)index);
}
if (!implicit_matches(split.first(), split.second())) {
printf("%s-5 %d", __FUNCTION__, (int)index);
const Quadratic* quads[] = {
&test, &midThird, &split.first(), &split.second()
};
size_t quadsCount = sizeof(quads) / sizeof(quads[0]);
for (size_t one = 0; one < quadsCount; ++one) {
for (size_t two = 0; two < quadsCount; ++two) {
for (size_t inner = 0; inner < 3; inner += 2) {
if (!point_on_parameterized_curve(*quads[one], (*quads[two])[inner])) {
printf("%s %zu [%zu,%zu] %zu parameterization failed\n",
__FUNCTION__, index, one, two, inner);
}
}
if (!implicit_matches(*quads[one], *quads[two])) {
printf("%s %zu [%zu,%zu] coincidence failed\n", __FUNCTION__,
index, one, two);
}
}
}
}
}

Просмотреть файл

@ -0,0 +1,15 @@
// included by QuadraticParameterization.cpp
// accesses internal functions to validate parameterized coefficients
bool point_on_parameterized_curve(const Quadratic& quad, const _Point& point) {
double coeffs[coeff_count];
implicit_coefficients(quad, coeffs);
double xx = coeffs[ xx_coeff] * point.x * point.x;
double xy = coeffs[ xy_coeff] * point.x * point.y;
double yy = coeffs[ yy_coeff] * point.y * point.y;
double x = coeffs[ x_coeff] * point.x;
double y = coeffs[ y_coeff] * point.y;
double c = coeffs[ c_coeff];
double sum = xx + xy + yy + x + y + c;
return approximately_zero(sum);
}

Просмотреть файл

@ -0,0 +1,38 @@
#include "CubicIntersection.h"
#include "Intersection_Tests.h"
#include "QuadraticIntersection_TestData.h"
#include "TestUtilities.h"
void QuadraticReduceOrder_Test() {
size_t index;
Quadratic reduce;
int order;
enum {
RunAll,
RunQuadraticLines,
RunQuadraticModLines,
RunNone
} run = RunAll;
int firstTestIndex = 0;
#if 0
run = RunQuadraticLines;
firstTestIndex = 1;
#endif
int firstQuadraticLineTest = run == RunAll ? 0 : run == RunQuadraticLines ? firstTestIndex : INT_MAX;
int firstQuadraticModLineTest = run == RunAll ? 0 : run == RunQuadraticModLines ? firstTestIndex : INT_MAX;
for (index = firstQuadraticLineTest; index < quadraticLines_count; ++index) {
const Quadratic& quad = quadraticLines[index];
order = reduceOrder(quad, reduce);
if (order != 2) {
printf("[%d] line quad order=%d\n", (int) index, order);
}
}
for (index = firstQuadraticModLineTest; index < quadraticModEpsilonLines_count; ++index) {
const Quadratic& quad = quadraticModEpsilonLines[index];
order = reduceOrder(quad, reduce);
if (order != 3) {
printf("[%d] line mod quad order=%d\n", (int) index, order);
}
}
}

Просмотреть файл

@ -0,0 +1,21 @@
#include "QuadraticUtilities.h"
int quadraticRoots(double A, double B, double C, double t[2]) {
B *= 2;
double square = B * B - 4 * A * C;
if (square < 0) {
return 0;
}
double squareRt = sqrt(square);
double Q = (B + (B < 0 ? -squareRt : squareRt)) / -2;
int foundRoots = 0;
if ((Q <= A) ^ (Q < 0)) {
t[foundRoots++] = Q / A;
}
if ((C <= Q) ^ (C < 0)) {
t[foundRoots++] = C / Q;
}
return foundRoots;
}

Просмотреть файл

@ -4,7 +4,7 @@
* b = 2*B - 2*C
* c = C
*/
static void set_abc(const double* quad, double& a, double& b, double& c) {
inline void set_abc(const double* quad, double& a, double& b, double& c) {
a = quad[0]; // a = A
b = 2 * quad[2]; // b = 2*B
c = quad[4]; // c = C
@ -12,3 +12,5 @@ static void set_abc(const double* quad, double& a, double& b, double& c) {
a -= b; // a = A - 2*B + C
b -= c; // b = 2*B - 2*C
}
int quadraticRoots(double A, double B, double C, double t[2]);

Просмотреть файл

@ -1,6 +1,7 @@
#include "CubicIntersection.h"
#include "CubicIntersection_TestData.h"
#include "CubicIntersection_Tests.h"
#include "QuadraticIntersection_TestData.h"
#include "TestUtilities.h"
void ReduceOrder_Test() {

Просмотреть файл

@ -67,3 +67,19 @@ void xy_at_t(const Cubic& cubic, double t, double& x, double& y) {
x = a * cubic[0].x + b * cubic[1].x + c * cubic[2].x + d * cubic[3].x;
y = a * cubic[0].y + b * cubic[1].y + c * cubic[2].y + d * cubic[3].y;
}
void xy_at_t(const _Line& line, double t, double& x, double& y) {
double one_t = 1 - t;
x = one_t * line[0].x + t * line[1].x;
y = one_t * line[0].y + t * line[1].y;
}
void xy_at_t(const Quadratic& quad, double t, double& x, double& y) {
double one_t = 1 - t;
double a = one_t * one_t;
double b = 2 * one_t * t;
double c = t * t;
x = a * quad[0].x + b * quad[1].x + c * quad[2].x;
y = a * quad[0].y + b * quad[1].y + c * quad[2].y;
}

Просмотреть файл

@ -4,3 +4,5 @@ bool controls_inside(const Cubic& );
void find_tight_bounds(const Cubic& , _Rect& );
void quad_to_cubic(const Quadratic& , Cubic& );
void xy_at_t(const Cubic& , double t, double& x, double& y);
void xy_at_t(const _Line& , double t, double& x, double& y);
void xy_at_t(const Quadratic& , double t, double& x, double& y);

Разница между файлами не показана из-за своего большого размера Загрузить разницу