conicsec.cpp revision 76addc201c409e81eaaa73fe27cc0f79c4db097c
/*
* Authors:
* Nathan Hurst <njh@njhurst.com
*
* Copyright 2009 authors
*
* modify it either under the terms of the GNU Lesser General Public
* License version 2.1 as published by the Free Software Foundation
* (the "LGPL") or, at your option, under the terms of the Mozilla
* Public License Version 1.1 (the "MPL"). If you do not alter this
* notice, a recipient may use your version of this file under either
* the MPL or the LGPL.
*
* You should have received a copy of the LGPL along with this library
* in the file COPYING-LGPL-2.1; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
* You should have received a copy of the MPL along with this library
* in the file COPYING-MPL-1.1
*
* The contents of this file are subject to the Mozilla Public License
* Version 1.1 (the "License"); you may not use this file except in
* compliance with the License. You may obtain a copy of the License at
*
* This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
* OF ANY KIND, either express or implied. See the LGPL or the MPL for
* the specific language governing rights and limitations.
*/
// File: convert.h
#include <utility>
#include <sstream>
#include <stdexcept>
namespace Geom
{
if (seg) {
return *seg;
} else {
}
}
return a[0]*b[1] - a[1]*b[0];
}
template <typename T>
static T det(T a, T b, T c, T d) {
return a*d - b*c;
}
template <typename T>
return M[0][0]*M[1][1] - M[1][0]*M[0][1];
}
template <typename T>
M[2][1], M[2][2])
M[2][1], M[2][2])
M[1][1], M[1][2]));
}
}
/**
* Find the roots of (q2x + q1)x+q0 = 0
* Tries to be numerically robust.
*/
template <typename T>
if(q2 == 0) {
if(q1 == 0) { // zero or infinite roots
return r;
}
} else {
/*cout << q2 << ", "
<< q1 << ", "
<< q0 << "; "
<< desc << "\n";*/
if (desc < 0)
return r;
else if (desc == 0)
else {
}
}
return r;
}
public:
: std::runtime_error(s)
{ }
};
template <typename T>
{
std::ostringstream o;
if (!(o << x))
throw BadConversion("stringify(T)");
return o.str();
}
/* A G4 continuous cubic parametric approximation for rational quadratics.
See
An analysis of cubic approximation schemes for conic sections
Michael Floater
SINTEF
This is less accurate overall than some of his other schemes, but
produces very smooth joins and is still optimally h^-6
convergent.
*/
}
Point P,
try {
if(!oc) // what to do?
//assert(0);
// std::cout << "RatQuad::fromPointsTangents: triarea = " << triarea << std::endl;
if (triarea == 0)
{
}
{
}
// std::cout << "RatQuad::fromPointsTangents: tau0 = " << tau0 << std::endl;
// std::cout << "RatQuad::fromPointsTangents: tau1 = " << tau1 << std::endl;
// std::cout << "RatQuad::fromPointsTangents: tau2 = " << tau2 << std::endl;
// std::cout << "RatQuad::fromPointsTangents: w = " << w << std::endl;
} catch(Geom::InfiniteSolutions) {
}
}
}
}
return CubicBezier(P[0],
P[2]);
}
}
a.P[0] = P[0];
b.P[2] = P[2];
a.P[1] = (P[0]+w*P[1])/(1+w);
b.P[1] = (w*P[1]+P[2])/(1+w);
a.P[2] = b.P[0] = (0.5*a.P[1]+0.5*b.P[1]);
}
}
return out;
}
return res;
}
#if 0
double M[3][3] = {{c[0], c[1], c[3]},
{c[1], c[2], c[4]},
{c[3], c[4], c[5]}};
double D = det3(M);
if (c[0] == 0 && c[1] == 0 && c[2] == 0)
return "line";
if (descr < 0) {
if (c[0] == c[2] && c[1] == 0)
return res + "circle";
return res + "ellipse";
} else if (descr == 0) {
return res + "parabola";
} else if (descr > 0) {
if (c[0] + c[2] == 0) {
if (D == 0)
return res + "two lines";
return res + "rectangular hyperbola";
}
return res + "hyperbola";
}
return "no idea!";
}
#endif
//Point B0 = xC0.bottom();
//std::cout << determ << "\n";
(-A[1][0]*b[0] + A[0][0]*b[1]));
// Are these just the eigenvectors of A11?
//assert(fabs(b*b-c) > 1e-10);
//assert(fabs(b-d) > 1e-10);
//assert(fabs(b*b-c) > 1e-10);
//assert(fabs(b-d) > 1e-10);
} else {
//assert(fabs(b*b-c) > 1e-10);
//assert(fabs(b-d) > 1e-10);
}
}
}
} else {
// single or double line
// check for completely zero case (what to do?)
if(L2sq(g) == 0) {
trial_pt[0] += 1;
if(L2sq(g) == 0) {
if(L2sq(g) == 0) {
trial_pt[0] += 1;
if(L2sq(g) == 0) {
}
}
}
}
//std::cout << trial_pt << ", " << g << "\n";
/**
* At this point we have tried up to 4 points: 0,0, 1,0, 1,1, 2,1, 1.5,1.5
*
* No degenerate conic can pass through these points, so we can assume
* that we've found a perpendicular to the double line.
* Proof:
* any degenerate must consist of at most 2 lines. 1.5,0.5 is not on any pair of lines
* passing through the previous 4 trials.
*
* alternatively, there may be a way to determine this directly from xC0
*/
//std::cout << P0 << "\n";
// It's very likely that at least one of the conics is degenerate, this will hopefully pick the more generate of the two.
else
}
}
}
return res;
}
{(C.c[1])/2, C.c[2], (C.c[4])/2},
{(C.c[3])/2, (C.c[4])/2, C.c[5]}};
}
// You know, if either of the inputs are degenerate we should use them first!
}
}
D = det3(C);
}
// at this point we have a T and S and perhaps some roots that represent our degenerate conic
// Let's just pick one randomly (can we do better?)
//for(unsigned i = 0; i < rts.size(); i++) {
unsigned i = 0;
//::draw(cr, xC0, screen_rect); // degen
} else {
;//std::cout << D << "\n";
}
return res;
}
}
return xAx();//1., 0, 1., -2*(1+d)*p[0], -2*(1+d)*p[1], dot(p,p)+d*d);
}
}
double dist;
}
vv[0] = P[0]*P[0];
vv[3] = P[0];
}
}
return evaluate_at(P[0], P[1]);
}
}
double x = p[0];
double y = p[1];
c[1]*x + 2*c[2]*y + c[4]);
}
for(int i = 0; i < 6; i++) {
res.c[i] = c[i] - b.c[i];
}
return res;
}
for(int i = 0; i < 6; i++) {
res.c[i] = c[i] + b.c[i];
}
return res;
}
for(int i = 0; i < 5; i++) {
res.c[i] = c[i];
}
return res;
}
for(int i = 0; i < 6; i++) {
res.c[i] = c[i] * b;
}
return res;
}
}
}
return res;
}
}
}
C = crs[2];
A = crs[2];
}
if(!bisect_rts.empty()) {
int besti = -1;
for(unsigned i =0; i < bisect_rts.size(); i++) {
besti = i;
}
}
if(besti >= 0) {
return RatQuad::fromPointsTangents(A, C-A, B, C, A-C);
}
return rq;
//std::vector<SBasis> hrq = rq.homogenous();
/*SBasis vertex_poly = evaluate_at(hrq[0], hrq[1], hrq[2]);
std::vector<double> rts = roots(vertex_poly);
for(unsigned i = 0; i < rts.size(); i++) {
//draw_circ(cr, Point(rq.pointAt(rts[i])));
}*/
}
}
}
}
// Find the roots on line l
// form the quadratic Q(t) = 0 by composing l with xAx
double q1 = (2*c[0]*d[0]*o[0] +
c[1]*(d[0]*o[1]+d[1]*o[0]) +
2*c[2]*d[1]*o[1] +
c[3]*d[0] + c[4]*d[1]);
if(q2 == 0) {
if(q1 == 0) {
return r;
}
} else {
/*std::cout << q2 << ", "
<< q1 << ", "
<< q0 << "; "
<< desc << "\n";*/
if (desc < 0)
return r;
else if (desc == 0)
else {
double t;
if (q1 == 0)
{
t = -0.5 * desc;
}
else
{
}
}
}
return r;
}
}
double cx = -b*0.5/a;
return bnds;
}
c[1], 2*c[2],
0, 0);
return m;
}
} else {
}
}
double A[2][2] = {{2*c[0], c[1]},
{c[1], 2*c[2]}};
double b[2] = {-c[3], -c[4]};
return solve(A, b);
//return Point(-c[3], -c[4])*hessian().inverse();
}
if (c[0] == 0 && c[1] == 0 && c[2] == 0) {
for(int i = 1; i < 4; i++)
return ext;
}
double k = r[X].min();
k = r[X].max();
k = r[Y].min();
k = r[Y].max();
return ext;
}
/*
* helper functions
*/
bool at_infinity (Point const& p)
{
{
return true;
}
return false;
}
inline
{
}
/*
* Define a conic section by computing the one that fits better with
* N points.
*
* points: points to fit
*
* precondition: there must be at least 5 non-overlapping points
*/
{
if (sz < 5)
{
THROW_RANGEERROR("fitting error: too few points passed");
}
{
}
}
/*
* Define a section conic by providing the coordinates of one of its vertex,
* the major axis inclination angle and the coordinates of its foci
* with respect to the unidimensional system defined by the major axis with
* origin set at the provided vertex.
*
* _vertex : section conic vertex V
* _angle : section conic major axis angle
* _dist1: +/-distance btw V and nearest focus
* _dist2: +/-distance btw V and farest focus
*
* prerequisite: _dist1 <= _dist2
*/
{
{
{
return;
}
// y^2 - 4px == 0
return;
}
{
}
if (_dist1 < 0)
{
}
// ellipse and hyperbola
// std::cout << "cA: " << cA << std::endl;
// std::cout << "cC: " << cC << std::endl;
// std::cout << "cF: " << cF << std::endl;
double CxCx = C[X] * C[X];
double CxCy = C[X] * C[Y];
double CyCy = C[Y] * C[Y];
}
/*
* Define a conic section by providing one of its vertex and its foci.
*
* _vertex: section conic vertex
* _focus1: section conic focus
* _focus2: section conic focus
*/
{
if (at_infinity(_vertex))
{
THROW_RANGEERROR("case not handled: vertex at infinity");
}
if (at_infinity(_focus2))
{
if (at_infinity(_focus1))
{
THROW_RANGEERROR("case not handled: both focus at infinity");
}
return;
}
else if (at_infinity(_focus1))
{
return;
}
{
// std::cout << "t = " << t << std::endl;
// std::cout << "dist2 = " << dist2 << std::endl;
}
{
double dist1 = 0;
}
else
{
}
}
/*
* Define a conic section by passing a focus, the related directrix,
* and the eccentricity (e)
* (e < 1 -> ellipse; e = 1 -> parabola; e > 1 -> hyperbola)
*
* _focus: a focus of the conic section
* _directrix: the directrix related to the given focus
* _eccentricity: the eccentricity parameter of the conic section
*/
{
//std::cout << "O = " << O << std::endl;
coeff(1) = 0;
coeff(4) = 0;
coeff(5) = p * p;
//std::cout << "O = " << O << std::endl;
(*this) = translate (O);
}
/*
* Made up a degenerate conic section as a pair of lines
*
* l1, l2: lines that made up the conic section
*/
{
}
/*
* Return the section conic kind
*/
{
//double T3 = det(C);
//std::cout << "T3 = " << T3 << std::endl;
//std::cout << "sT3 = " << sT3 << std::endl;
//std::cout << "t2 = " << t2 << std::endl;
//std::cout << "t1 = " << t1 << std::endl;
//std::cout << "st2 = " << st2 << std::endl;
if (sT3 != 0)
{
if (st2 == 0)
{
return PARABOLA;
}
else if (st2 == 1)
{
{
//std::cout << "t1 * t1 - 4 * t2 = "
// << (t1 * t1 - 4 * t2) << std::endl;
//std::cout << "discr_sgn = " << discr_sgn << std::endl;
if (discr_sgn == 0)
{
return CIRCLE;
}
else
{
return REAL_ELLIPSE;
}
}
else // sT3 * st1 > 0
{
return IMAGINARY_ELLIPSE;
}
}
else // t2 < 0
{
if (st1 == 0)
{
return RECTANGULAR_HYPERBOLA;
}
else
{
return HYPERBOLA;
}
}
}
else // T3 == 0
{
if (st2 == 0)
{
//double T2 = NL::trace<2>(C);
//std::cout << "T2 = " << T2 << std::endl;
//std::cout << "sT2 = " << sT2 << std::endl;
if (sT2 == 0)
{
return DOUBLE_LINE;
}
if (sT2 == -1)
{
return TWO_REAL_PARALLEL_LINES;
}
else // T2 > 0
{
return TWO_IMAGINARY_PARALLEL_LINES;
}
}
else if (st2 == -1)
{
return TWO_REAL_CROSSING_LINES;
}
else // t2 > 0
{
return TWO_IMAGINARY_CROSSING_LINES;
}
}
return UNKNOWN;
}
/*
* Return a string representing the conic section kind
*/
{
switch (KIND)
{
case PARABOLA :
return "parabola";
case CIRCLE :
return "circle";
case REAL_ELLIPSE :
return "real ellispe";
case IMAGINARY_ELLIPSE :
return "imaginary ellispe";
case RECTANGULAR_HYPERBOLA :
return "rectangular hyperbola";
case HYPERBOLA :
return "hyperbola";
case DOUBLE_LINE :
return "double line";
case TWO_REAL_PARALLEL_LINES :
return "two real parallel lines";
case TWO_IMAGINARY_PARALLEL_LINES :
return "two imaginary parallel lines";
case TWO_REAL_CROSSING_LINES :
return "two real crossing lines";
case TWO_IMAGINARY_CROSSING_LINES :
return "two imaginary crossing lines";
default :
return "unknown";
}
}
/*
* Compute the solutions of the conic section algebraic equation with respect to
* one coordinate after substituting to the other coordinate the passed value
*
* sol: the computed solutions
* v: the provided value
* d: the index of the coordinate the passed value have to be substituted to
*/
{
if (d < 0 || d > Y)
{
THROW_RANGEERROR("dimension parameter out of range");
}
// p*t^2 + q*t + r = 0;
double p, q, r;
if (d == X)
{
p = coeff(2);
}
else
{
p = coeff(0);
}
if (p == 0)
{
if (q == 0) return;
double t = -r/q;
return;
}
if (q == 0)
{
if ((p > 0 && r > 0) || (p < 0 && r < 0)) return;
double t = -r / p;
return;
}
if (r == 0)
{
double t = -q/p;
return;
}
//std::cout << "p = " << p << ", q = " << q << ", r = " << r << std::endl;
double delta = q * q - 4 * p * r;
if (delta < 0) return;
if (delta == 0)
{
double t = -q / (2 * p);
return;
}
// else
}
/*
* Return the inclination angle of the major axis of the conic section
*/
double xAx::axis_angle() const
{
{
return l.angle();
}
double angle;
if (sgn_discr == 0)
{
//std::cout << "rotation_angle: sgn_discr = "
// << sgn_discr << std::endl;
}
else
{
angle /= 2;
}
//std::cout << "rotation_angle : angle = " << angle << std::endl;
return angle;
}
/*
* Translate the conic section by the given vector offset
*
* _offset: represent the vector offset
*/
{
return cs;
}
/*
* Rotate the conic section by the given angle wrt the point (0,0)
*
* angle: represent the rotation angle
*/
{
double cc = c * c;
double ss = s * s;
double cs = c * s;
// quadratic terms
// linear terms
return result;
}
/*
* Decompose a degenerate conic in two lines the conic section is made by.
* Return true if the decomposition is successfull, else if it fails.
*
* l1, l2: out parameters where the decomposed conic section is returned
*/
{
if (!is_quadratic() || !isDegenerate())
{
return false;
}
if (!D.is_zero()) // D == 0 <=> rank(C) < 2
{
//if (D.get<0,0>() < 0 || D.get<1,1>() < 0 || D.get<2,2>() < 0)
//{
//std::cout << "C: \n" << C << std::endl;
//std::cout << "D: \n" << D << std::endl;
/*
* This case should be impossible because any diagonal element
* of D is a square, but due to non exact aritmethic computation
* it can actually happen; however the algorithm seems to work
* correctly even if some diagonal term is negative, the only
* difference is that we should compute the absolute value of
* diagonal elements. So until we elaborate a better degenerate
* test it's better not rising exception when we have a negative
* diagonal element.
*/
//}
if (d[idx] == 0)
{
THROW_LOGICALERROR ("xAx::decompose: "
"rank 2 but adjoint with null diagonal");
}
M(1,2) += d[0]; M(2,1) -= d[0];
M(0,2) -= d[1]; M(2,0) += d[1];
M(0,1) += d[2]; M(1,0) -= d[2];
//std::cout << "C: \n" << C << std::endl;
//std::cout << "D: \n" << D << std::endl;
//std::cout << "d = " << d << std::endl;
//std::cout << "M = " << M << std::endl;
}
{
}
else
{
}
return true;
}
/*
* Return the rectangle that bound the conic section arc characterized by
* the passed points.
*
* P1: the initial point of the arc
* Q: the inner point of the arc
* P2: the final point of the arc
*
* prerequisite: the passed points must lie on the conic
*/
{
//std::cout << "BOUND: P1 = " << P1 << std::endl;
//std::cout << "BOUND: Q = " << Q << std::endl;
//std::cout << "BOUND: P2 = " << P2 << std::endl;
//std::cout << "BOUND: Qside = " << Qside << std::endl;
bool empty[2] = {false, false};
try // if the passed coefficients lead to an equation 0x + 0y + c == 0,
{ // with c != 0 the setCoefficients rise an exception
}
catch(Geom::LogicalError const &e)
{
empty[0] = true;
}
try
{
}
catch(Geom::LogicalError const &e)
{
empty[1] = true;
}
{
M.clear();
if (M.size() == 1)
{
{
//std::cout << "BOUND: M.size() == 1" << std::endl;
}
}
else if (M.size() == 2)
{
//std::cout << "BOUND: M.size() == 2" << std::endl;
swap (M[0], M[1]);
{
}
{
}
else
{
}
}
}
return B;
}
/*
* Return all points on the conic section nearest to the passed point "P".
*
* P: the point to compute the nearest one
*/
{
// TODO: manage the circle - centre case
// named C the conic we look for points (x,y) on C such that
// dot (grad (C(x,y)), rot90 (P -(x,y))) == 0; the set of points satisfying
// this equation is still a conic G, so the wanted points can be found by
// intersecting C with G
coeff(1),
{
{
idx = i;
}
}
{
if (i == idx) continue;
}
return points;
}
{
}
} // end namespace Geom
/*
Local Variables:
mode:c++
c-file-style:"stroustrup"
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
indent-tabs-mode:nil
fill-column:99
End:
*/
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :