2 This file is part of Kig, a KDE program for Interactive Geometry...
3 Copyright (C) 2002 Maurizio Paolini <paolini@dmf.unicatt.it>
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
23 #include "conic-common.h"
26 #include "kigtransform.h"
35 ConicCartesianData::ConicCartesianData(
36 const ConicPolarData
& polardata
39 double ec
= polardata
.ecostheta0
;
40 double es
= polardata
.esintheta0
;
41 double p
= polardata
.pdimen
;
42 double fx
= polardata
.focus1
.x
;
43 double fy
= polardata
.focus1
.y
;
52 f
+= a
*fx
*fx
+ b
*fy
*fy
+ c
*fx
*fy
- d
*fx
- e
*fy
;
64 ConicPolarData::ConicPolarData( const ConicCartesianData
& cartdata
)
66 double a
= cartdata
.coeffs
[0];
67 double b
= cartdata
.coeffs
[1];
68 double c
= cartdata
.coeffs
[2];
69 double d
= cartdata
.coeffs
[3];
70 double e
= cartdata
.coeffs
[4];
71 double f
= cartdata
.coeffs
[5];
73 // 1. Compute theta (tilt of conic)
74 double theta
= std::atan2(c
, b
- a
)/2;
76 // now I should possibly add pi/2...
77 double costheta
= std::cos(theta
);
78 double sintheta
= std::sin(theta
);
79 // compute new coefficients (c should now be zero)
80 double aa
= a
*costheta
*costheta
+ b
*sintheta
*sintheta
- c
*sintheta
*costheta
;
81 double bb
= a
*sintheta
*sintheta
+ b
*costheta
*costheta
+ c
*sintheta
*costheta
;
83 { // we have a hyperbola we need to check the correct orientation
84 double dd
= d
*costheta
- e
*sintheta
;
85 double ee
= d
*sintheta
+ e
*costheta
;
86 double xc
= - dd
/ ( 2*aa
);
87 double yc
= - ee
/ ( 2*bb
);
88 double ff
= f
+ aa
*xc
*xc
+ bb
*yc
*yc
+ dd
*xc
+ ee
*yc
;
89 if (ff
*aa
> 0) // wrong orientation
91 if (theta
> 0) theta
-= M_PI
/2;
93 costheta
= cos(theta
);
94 sintheta
= sin(theta
);
95 aa
= a
*costheta
*costheta
+ b
*sintheta
*sintheta
- c
*sintheta
*costheta
;
96 bb
= a
*sintheta
*sintheta
+ b
*costheta
*costheta
+ c
*sintheta
*costheta
;
101 if ( std::fabs (bb
) < std::fabs (aa
) )
103 if (theta
> 0) theta
-= M_PI
/2;
104 else theta
+= M_PI
/2;
105 costheta
= cos(theta
);
106 sintheta
= sin(theta
);
107 aa
= a
*costheta
*costheta
+ b
*sintheta
*sintheta
- c
*sintheta
*costheta
;
108 bb
= a
*sintheta
*sintheta
+ b
*costheta
*costheta
+ c
*sintheta
*costheta
;
112 double cc
= 2*(a
- b
)*sintheta
*costheta
+
113 c
*(costheta
*costheta
- sintheta
*sintheta
);
114 // cc should be zero!!! cout << "cc = " << cc << "\n";
115 double dd
= d
*costheta
- e
*sintheta
;
116 double ee
= d
*sintheta
+ e
*costheta
;
124 // now b cannot be zero (otherwise conic is degenerate)
132 // 2. compute y coordinate of focuses
138 e
+= 2*yf
; // this should be zero!
140 // now: a > 0 -> ellipse
142 // a < 0 -> hyperbola
144 double eccentricity
= sqrt(1.0 - a
);
146 double sqrtdiscrim
= sqrt(d
*d
- 4*a
*f
);
147 if (d
< 0.0) sqrtdiscrim
= -sqrtdiscrim
;
148 double xf
= (4*a
*f
- 4*f
- d
*d
)/(d
+ eccentricity
*sqrtdiscrim
) / 2;
150 // 3. the focus needs to be rotated back into position
151 focus1
.x
= xf
*costheta
+ yf
*sintheta
;
152 focus1
.y
= -xf
*sintheta
+ yf
*costheta
;
154 // 4. final touch: the pdimen
155 pdimen
= -sqrtdiscrim
/2;
157 ecostheta0
= eccentricity
*costheta
;
158 esintheta0
= -eccentricity
*sintheta
;
162 ecostheta0
= -ecostheta0
;
163 esintheta0
= -esintheta0
;
167 const ConicCartesianData
calcConicThroughPoints (
168 const std::vector
<Coordinate
>& points
,
169 const LinearConstraints c1
,
170 const LinearConstraints c2
,
171 const LinearConstraints c3
,
172 const LinearConstraints c4
,
173 const LinearConstraints c5
)
175 assert( 0 < points
.size() && points
.size() <= 5 );
176 // points is a vector of up to 5 points through which the conic is
178 // this routine should compute the coefficients in the cartesian equation
179 // a x^2 + b y^2 + c xy + d x + e y + f = 0
180 // they are defined up to a multiplicative factor.
181 // since we don't know (in advance) which one of them is nonzero, we
182 // simply keep all 6 parameters, obtaining a 5x6 linear system which
183 // we solve using gaussian elimination with complete pivoting
184 // If there are too few, then we choose some cool way to fill in the
185 // empty parts in the matrix according to the LinearConstraints
188 // 5 rows, 6 columns..
194 double *matrix
[5] = {row0
, row1
, row2
, row3
, row4
};
197 LinearConstraints constraints
[] = {c1
, c2
, c3
, c4
, c5
};
199 int numpoints
= points
.size();
200 int numconstraints
= 5;
202 // fill in the matrix elements
203 for ( int i
= 0; i
< numpoints
; ++i
)
205 double xi
= points
[i
].x
;
206 double yi
= points
[i
].y
;
207 matrix
[i
][0] = xi
*xi
;
208 matrix
[i
][1] = yi
*yi
;
209 matrix
[i
][2] = xi
*yi
;
215 for ( int i
= 0; i
< numconstraints
; i
++ )
217 if (numpoints
>= 5) break; // don't add constraints if we have enough
218 for (int j
= 0; j
< 6; ++j
) matrix
[numpoints
][j
] = 0.0;
219 // force the symmetry axes to be
220 // parallel to the coordinate system (zero tilt): c = 0
221 if (constraints
[i
] == zerotilt
) matrix
[numpoints
][2] = 1.0;
222 // force a parabula (if zerotilt): b = 0
223 if (constraints
[i
] == parabolaifzt
) matrix
[numpoints
][1] = 1.0;
224 // force a circle (if zerotilt): a = b
225 if (constraints
[i
] == circleifzt
) {
226 matrix
[numpoints
][0] = 1.0;
227 matrix
[numpoints
][1] = -1.0; }
228 // force an equilateral hyperbola: a + b = 0
229 if (constraints
[i
] == equilateral
) {
230 matrix
[numpoints
][0] = 1.0;
231 matrix
[numpoints
][1] = 1.0; }
232 // force symmetry about y-axis: d = 0
233 if (constraints
[i
] == ysymmetry
) matrix
[numpoints
][3] = 1.0;
234 // force symmetry about x-axis: e = 0
235 if (constraints
[i
] == xsymmetry
) matrix
[numpoints
][4] = 1.0;
237 if (constraints
[i
] != noconstraint
) ++numpoints
;
240 if ( ! GaussianElimination( matrix
, numpoints
, 6, scambio
) )
241 return ConicCartesianData::invalidData();
242 // fine della fase di eliminazione
243 BackwardSubstitution( matrix
, numpoints
, 6, scambio
, solution
);
245 // now solution should contain the correct coefficients..
246 return ConicCartesianData( solution
);
249 const ConicPolarData
calcConicBFFP(
250 const std::vector
<Coordinate
>& args
,
253 assert( args
.size() >= 2 && args
.size() <= 3 );
254 assert( type
== 1 || type
== -1 );
258 Coordinate f1
= args
[0];
259 Coordinate f2
= args
[1];
261 double eccentricity
, d1
, d2
, dl
;
263 Coordinate f2f1
= f2
- f1
;
264 double f2f1l
= f2f1
.length();
265 ret
.ecostheta0
= f2f1
.x
/ f2f1l
;
266 ret
.esintheta0
= f2f1
.y
/ f2f1l
;
268 if ( args
.size() == 3 )
271 d1
= ( d
- f1
).length();
272 d2
= ( d
- f2
).length();
273 dl
= fabs( d1
+ type
* d2
);
274 eccentricity
= f2f1l
/dl
;
278 if ( type
> 0 ) eccentricity
= 0.7; else eccentricity
= 2.0;
279 dl
= f2f1l
/eccentricity
;
282 double rhomax
= (dl
+ f2f1l
) /2.0;
284 ret
.ecostheta0
*= eccentricity
;
285 ret
.esintheta0
*= eccentricity
;
286 ret
.pdimen
= type
*(1 - eccentricity
)*rhomax
;
287 ret
.focus1
= type
== 1 ? f1
: f2
;
291 const LineData
calcConicPolarLine (
292 const ConicCartesianData
& data
,
293 const Coordinate
& cpole
,
298 double a
= data
.coeffs
[0];
299 double b
= data
.coeffs
[1];
300 double c
= data
.coeffs
[2];
301 double d
= data
.coeffs
[3];
302 double e
= data
.coeffs
[4];
303 double f
= data
.coeffs
[5];
305 double alpha
= 2*a
*x
+ c
*y
+ d
;
306 double beta
= c
*x
+ 2*b
*y
+ e
;
307 double gamma
= d
*x
+ e
*y
+ 2*f
;
309 double normsq
= alpha
*alpha
+ beta
*beta
;
311 if (normsq
< 1e-10) // line at infinity
318 Coordinate reta
= -gamma
/normsq
* Coordinate (alpha
, beta
);
319 Coordinate retb
= reta
+ Coordinate (-beta
, alpha
);
320 return LineData( reta
, retb
);
323 const Coordinate
calcConicPolarPoint (
324 const ConicCartesianData
& data
,
325 const LineData
& polar
)
327 Coordinate p1
= polar
.a
;
328 Coordinate p2
= polar
.b
;
330 double alpha
= p2
.y
- p1
.y
;
331 double beta
= p1
.x
- p2
.x
;
332 double gamma
= p1
.y
*p2
.x
- p1
.x
*p2
.y
;
334 double a11
= data
.coeffs
[0];
335 double a22
= data
.coeffs
[1];
336 double a12
= data
.coeffs
[2]/2.0;
337 double a13
= data
.coeffs
[3]/2.0;
338 double a23
= data
.coeffs
[4]/2.0;
339 double a33
= data
.coeffs
[5];
341 // double detA = a11*a22*a33 - a11*a23*a23 - a22*a13*a13 - a33*a12*a12 +
344 double a11inv
= a22
*a33
- a23
*a23
;
345 double a22inv
= a11
*a33
- a13
*a13
;
346 double a33inv
= a11
*a22
- a12
*a12
;
347 double a12inv
= a23
*a13
- a12
*a33
;
348 double a23inv
= a12
*a13
- a23
*a11
;
349 double a13inv
= a12
*a23
- a13
*a22
;
351 double x
= a11inv
*alpha
+ a12inv
*beta
+ a13inv
*gamma
;
352 double y
= a12inv
*alpha
+ a22inv
*beta
+ a23inv
*gamma
;
353 double z
= a13inv
*alpha
+ a23inv
*beta
+ a33inv
*gamma
;
355 if (fabs(z
) < 1e-10) // point at infinity
357 return Coordinate::invalidCoord();
362 return Coordinate (x
, y
);
365 const Coordinate
calcConicLineIntersect( const ConicCartesianData
& c
,
370 assert( which
== 1 || which
== -1 || which
== 0 );
372 double aa
= c
.coeffs
[0];
373 double bb
= c
.coeffs
[1];
374 double cc
= c
.coeffs
[2];
375 double dd
= c
.coeffs
[3];
376 double ee
= c
.coeffs
[4];
377 double ff
= c
.coeffs
[5];
381 double dx
= l
.b
.x
- l
.a
.x
;
382 double dy
= l
.b
.y
- l
.a
.y
;
384 double aaa
= aa
*dx
*dx
+ bb
*dy
*dy
+ cc
*dx
*dy
;
385 double bbb
= 2*aa
*x
*dx
+ 2*bb
*y
*dy
+ cc
*x
*dy
+ cc
*y
*dx
+ dd
*dx
+ ee
*dy
;
386 double ccc
= aa
*x
*x
+ bb
*y
*y
+ cc
*x
*y
+ dd
*x
+ ee
*y
+ ff
;
389 if ( which
== 0 ) /* i.e. we have a known intersection */
391 t
= - bbb
/aaa
- knownparam
;
392 return l
.a
+ t
*(l
.b
- l
.a
);
395 double discrim
= bbb
*bbb
- 4*aaa
*ccc
;
398 return Coordinate::invalidCoord();
404 t
= bbb
+ which
*sqrt(discrim
);
407 t
= -bbb
+ which
*sqrt(discrim
);
411 return l
.a
+ t
*(l
.b
- l
.a
);
415 ConicPolarData::ConicPolarData(
416 const Coordinate
& f
, double d
,
417 double ec
, double es
)
418 : focus1( f
), pdimen( d
), ecostheta0( ec
), esintheta0( es
)
422 ConicPolarData::ConicPolarData()
423 : focus1(), pdimen( 0 ), ecostheta0( 0 ), esintheta0( 0 )
427 const ConicPolarData
calcConicBDFP(
428 const LineData
& directrix
,
429 const Coordinate
& cfocus
,
430 const Coordinate
& cpoint
)
434 Coordinate ba
= directrix
.dir();
435 double bal
= ba
.length();
436 ret
.ecostheta0
= -ba
.y
/ bal
;
437 ret
.esintheta0
= ba
.x
/ bal
;
439 Coordinate pa
= cpoint
- directrix
.a
;
441 double distpf
= (cpoint
- cfocus
).length();
442 double distpd
= ( pa
.y
*ba
.x
- pa
.x
*ba
.y
)/bal
;
444 double eccentricity
= distpf
/distpd
;
445 ret
.ecostheta0
*= eccentricity
;
446 ret
.esintheta0
*= eccentricity
;
448 Coordinate fa
= cfocus
- directrix
.a
;
449 ret
.pdimen
= ( fa
.y
*ba
.x
- fa
.x
*ba
.y
)/bal
;
450 ret
.pdimen
*= eccentricity
;
456 ConicCartesianData::ConicCartesianData( const double incoeffs
[6] )
458 std::copy( incoeffs
, incoeffs
+ 6, coeffs
);
461 const LineData
calcConicAsymptote(
462 const ConicCartesianData data
,
463 int which
, bool &valid
)
465 assert( which
== -1 || which
== 1 );
468 double a
=data
.coeffs
[0];
469 double b
=data
.coeffs
[1];
470 double c
=data
.coeffs
[2];
471 double d
=data
.coeffs
[3];
472 double e
=data
.coeffs
[4];
474 double normabc
= a
*a
+ b
*b
+ c
*c
;
475 double delta
= c
*c
- 4*a
*b
;
476 if (fabs(delta
) < 1e-6*normabc
) { valid
= false; return ret
; }
478 double yc
= (2*a
*e
- c
*d
)/delta
;
479 double xc
= (2*b
*d
- c
*e
)/delta
;
480 // let c be nonnegative; we no longer need d, e, f.
494 double sqrtdelta
= sqrt(delta
);
495 Coordinate displacement
;
497 displacement
= Coordinate(-2*b
, c
+ sqrtdelta
);
499 displacement
= Coordinate(c
+ sqrtdelta
, -2*a
);
500 ret
.a
= Coordinate(xc
, yc
);
501 ret
.b
= ret
.a
+ displacement
;
505 const ConicCartesianData
calcConicByAsymptotes(
506 const LineData
& line1
,
507 const LineData
& line2
,
508 const Coordinate
& p
)
510 Coordinate p1
= line1
.a
;
511 Coordinate p2
= line1
.b
;
515 double c1
= p1
.x
*p2
.y
- p2
.x
*p1
.y
;
516 double b1
= p2
.x
- p1
.x
;
517 double a1
= p1
.y
- p2
.y
;
522 double c2
= p1
.x
*p2
.y
- p2
.x
*p1
.y
;
523 double b2
= p2
.x
- p1
.x
;
524 double a2
= p1
.y
- p2
.y
;
528 double c
= a1
*b2
+ a2
*b1
;
529 double d
= a1
*c2
+ a2
*c1
;
530 double e
= b1
*c2
+ c1
*b2
;
532 double f
= a
*x
*x
+ b
*y
*y
+ c
*x
*y
+ d
*x
+ e
*y
;
535 return ConicCartesianData( a
, b
, c
, d
, e
, f
);
538 const LineData
calcConicRadical( const ConicCartesianData
& cequation1
,
539 const ConicCartesianData
& cequation2
,
540 int which
, int zeroindex
, bool& valid
)
542 assert( which
== 1 || which
== -1 );
543 assert( 0 < zeroindex
&& zeroindex
< 4 );
547 double a
= cequation1
.coeffs
[0];
548 double b
= cequation1
.coeffs
[1];
549 double c
= cequation1
.coeffs
[2];
550 double d
= cequation1
.coeffs
[3];
551 double e
= cequation1
.coeffs
[4];
552 double f
= cequation1
.coeffs
[5];
554 double a2
= cequation2
.coeffs
[0];
555 double b2
= cequation2
.coeffs
[1];
556 double c2
= cequation2
.coeffs
[2];
557 double d2
= cequation2
.coeffs
[3];
558 double e2
= cequation2
.coeffs
[4];
559 double f2
= cequation2
.coeffs
[5];
561 // background: the family of conics c + lambda*c2 has members that
562 // degenerate into a union of two lines. The values of lambda giving
563 // such degenerate conics is the solution of a third degree equation.
564 // The coefficients of such equation are given by:
565 // (Thanks to Dominique Devriese for the suggestion of this approach)
566 // domi: (And thanks to Maurizio for implementing it :)
568 double df
= 4*a
*b
*f
- a
*e
*e
- b
*d
*d
- c
*c
*f
+ c
*d
*e
;
569 double cf
= 4*a2
*b
*f
+ 4*a
*b2
*f
+ 4*a
*b
*f2
570 - 2*a
*e
*e2
- 2*b
*d
*d2
- 2*f
*c
*c2
571 - a2
*e
*e
- b2
*d
*d
- f2
*c
*c
572 + c2
*d
*e
+ c
*d2
*e
+ c
*d
*e2
;
573 double bf
= 4*a
*b2
*f2
+ 4*a2
*b
*f2
+ 4*a2
*b2
*f
574 - 2*a2
*e2
*e
- 2*b2
*d2
*d
- 2*f2
*c2
*c
575 - a
*e2
*e2
- b
*d2
*d2
- f
*c2
*c2
576 + c
*d2
*e2
+ c2
*d
*e2
+ c2
*d2
*e
;
577 double af
= 4*a2
*b2
*f2
- a2
*e2
*e2
- b2
*d2
*d2
- c2
*c2
*f2
+ c2
*d2
*e2
;
579 // assume both conics are nondegenerate, renormalize so that af = 1
584 af
= 1.0; // not needed, just for consistency
586 // computing the coefficients of the Sturm sequence
588 double p1a
= 2*bf
*bf
- 6*cf
;
589 double p1b
= bf
*cf
- 9*df
;
590 double p0a
= cf
*p1a
*p1a
+ p1b
*(3*p1b
- 2*bf
*p1a
);
591 double fval
, fpval
, lambda
;
593 if (p0a
< 0 && p1a
< 0)
600 lambda
= -bf
/3.0; //inflection point
601 double displace
= 1.0;
602 if (p1a
> 0) // with two stationary points
604 displace
+= sqrt(p1a
); // should be enough. The important
605 // thing is that it is larger than the
606 // semidistance between the stationary points
608 // compute the value at the inflection point using Horner scheme
609 fval
= bf
+ lambda
; // b + x
610 fval
= cf
+ lambda
*fval
; // c + xb + xx
611 fval
= df
+ lambda
*fval
; // d + xc + xxb + xxx
613 if (fabs(p0a
) < 1e-7)
614 { // this is the case if we intersect two vertical parabulas!
615 p0a
= 1e-7; // fall back to the one zero case
619 // we have three roots..
620 // we select the one we want ( defined by mzeroindex.. )
621 lambda
+= ( 2 - zeroindex
)* displace
;
625 // we have just one root
626 if( zeroindex
> 1 ) // cannot find second and third root
632 if (fval
> 0) // zero on the left
635 } else { // zero on the right
639 // p0a = 0 means we have a root with multiplicity two
643 // find a root of af*lambda^3 + bf*lambda^2 + cf*lambda + df = 0
644 // (use a Newton method starting from lambda = 0. Hope...)
650 const int maxiterations
= 30;
651 while (iterations
++ < maxiterations
) // using Newton, iterations should be very few
653 // compute value of function and polinomial
655 fval
= bf
+ lambda
*fval
; // b + xa
656 fpval
= fval
+ lambda
*fpval
; // b + 2xa
657 fval
= cf
+ lambda
*fval
; // c + xb + xxa
658 fpval
= fval
+ lambda
*fpval
; // c + 2xb + 3xxa
659 fval
= df
+ lambda
*fval
; // d + xc + xxb + xxxa
663 if (fabs(delta
) < 1e-6) break;
665 if (iterations
>= maxiterations
) { valid
= false; return ret
; }
667 // now we have the degenerate conic: a, b, c, d, e, f
677 // this is the determinant of the matrix of the new conic. It
678 // should be zero, for the new conic to be degenerate.
679 df
= 4*a
*b
*f
- a
*e
*e
- b
*d
*d
- c
*c
*f
+ c
*d
*e
;
681 //lets work in homogeneous coordinates...
683 double dis1
= e
*e
- 4*b
*f
;
684 double maxval
= fabs(dis1
);
686 double dis2
= d
*d
- 4*a
*f
;
687 if (fabs(dis2
) > maxval
)
692 double dis3
= c
*c
- 4*a
*b
;
693 if (fabs(dis3
) > maxval
)
698 // one of these must be nonzero (otherwise the matrix is ...)
699 // exchange elements so that the largest is the determinant of the
704 case 1: // exchange 1 <-> 3
705 temp
= a
; a
= f
; f
= temp
;
706 temp
= c
; c
= e
; e
= temp
;
707 temp
= dis1
; dis1
= dis3
; dis3
= temp
;
710 case 2: // exchange 2 <-> 3
711 temp
= b
; b
= f
; f
= temp
;
712 temp
= c
; c
= d
; d
= temp
;
713 temp
= dis2
; dis2
= dis3
; dis3
= temp
;
718 // this is the negative of the determinant of the top left of the
719 // matrix. If it is 0, then the conic is a parabola, if it is < 0,
720 // then the conic is an ellipse, if positive, the conic is a
721 // hyperbola. In this case, it should be positive, since we have a
722 // degenerate conic, which is a degenerate case of a hyperbola..
723 // note that it is negative if there is no valid conic to be
724 // found ( e.g. not enough intersections.. )
725 // double discrim = c*c - 4*a*b;
730 // i would put an assertion here, but well, i guess it doesn't
731 // really matter, and this prevents crashes if the math i still
732 // recall from high school happens to be wrong :)
737 double r
[3]; // direction of the null space
742 // now remember the switch:
745 case 1: // exchange 1 <-> 3
746 temp
= a
; a
= f
; f
= temp
;
747 temp
= c
; c
= e
; e
= temp
;
748 temp
= dis1
; dis1
= dis3
; dis3
= temp
;
749 temp
= r
[0]; r
[0] = r
[2]; r
[2] = temp
;
752 case 2: // exchange 2 <-> 3
753 temp
= b
; b
= f
; f
= temp
;
754 temp
= c
; c
= d
; d
= temp
;
755 temp
= dis2
; dis2
= dis3
; dis3
= temp
;
756 temp
= r
[1]; r
[1] = r
[2]; r
[2] = temp
;
760 // Computing a Householder reflection transformation that
761 // maps r onto [0, 0, k]
764 double rnormsq
= r
[0]*r
[0] + r
[1]*r
[1] + r
[2]*r
[2];
765 double k
= sqrt( rnormsq
);
766 if ( k
*r
[2] < 0) k
= -k
;
767 double wnorm
= sqrt( 2*rnormsq
+ 2*k
*r
[2] );
770 w
[2] = (r
[2] + k
)/wnorm
;
772 // matrix transformation using Householder matrix, the resulting
773 // matrix is zero on third row and column
774 // [q0,q1,q2]^t = A w
776 double q0
= a
*w
[0] + c
*w
[1]/2 + d
*w
[2]/2;
777 double q1
= b
*w
[1] + c
*w
[0]/2 + e
*w
[2]/2;
778 double alpha
= a
*w
[0]*w
[0] + b
*w
[1]*w
[1] + c
*w
[0]*w
[1] +
779 d
*w
[0]*w
[2] + e
*w
[1]*w
[2] + f
*w
[2]*w
[2];
780 double a00
= a
- 4*w
[0]*q0
+ 4*w
[0]*w
[0]*alpha
;
781 double a11
= b
- 4*w
[1]*q1
+ 4*w
[1]*w
[1]*alpha
;
782 double a01
= c
/2 - 2*w
[1]*q0
- 2*w
[0]*q1
+ 4*w
[0]*w
[1]*alpha
;
784 double dis
= a01
*a01
- a00
*a11
;
786 double sqrtdis
= sqrt( dis
);
790 px
= a01
+ which
*sqrtdis
;
794 py
= a01
- which
*sqrtdis
;
796 double p
[3]; // vector orthogonal to one of the two planes
797 double pscalw
= w
[0]*px
+ w
[1]*py
;
798 p
[0] = px
- 2*pscalw
*w
[0];
799 p
[1] = py
- 2*pscalw
*w
[1];
800 p
[2] = - 2*pscalw
*w
[2];
802 // "r" is the solution of the equation A*(x,y,z) = (0,0,0) where
803 // A is the matrix of the degenerate conic. This is what we
804 // called in the conic theory we saw in high school a "double
805 // point". It has the unique property that any line going through
806 // it is a "polar line" of the conic at hand. It only exists for
807 // degenerate conics. It has another unique property that if you
808 // take any other point on the conic, then the line between it and
809 // the double point is part of the conic.
810 // this is what we use here: we find the double point ( ret.a
811 // ), and then find another points on the conic.
813 ret
.a
= -p
[2]/(p
[0]*p
[0] + p
[1]*p
[1]) * Coordinate (p
[0],p
[1]);
814 ret
.b
= ret
.a
+ Coordinate (-p
[1],p
[0]);
820 const ConicCartesianData
calcConicTransformation (
821 const ConicCartesianData
& data
, const Transformation
& t
, bool& valid
)
826 a
[1][1] = data
.coeffs
[0];
827 a
[2][2] = data
.coeffs
[1];
828 a
[1][2] = a
[2][1] = data
.coeffs
[2]/2;
829 a
[0][1] = a
[1][0] = data
.coeffs
[3]/2;
830 a
[0][2] = a
[2][0] = data
.coeffs
[4]/2;
831 a
[0][0] = data
.coeffs
[5];
833 Transformation ti
= t
.inverse( valid
);
834 if ( ! valid
) return ConicCartesianData();
836 double supnorm
= 0.0;
837 for (int i
= 0; i
< 3; i
++)
839 for (int j
= 0; j
< 3; j
++)
842 for (int ii
= 0; ii
< 3; ii
++)
844 for (int jj
= 0; jj
< 3; jj
++)
846 b
[i
][j
] += a
[ii
][jj
]*ti
.data( ii
, i
)*ti
.data( jj
, j
);
849 if ( std::fabs( b
[i
][j
] ) > supnorm
) supnorm
= std::fabs( b
[i
][j
] );
853 for (int i
= 0; i
< 3; i
++)
855 for (int j
= 0; j
< 3; j
++)
861 return ConicCartesianData ( b
[1][1], b
[2][2], b
[1][2] + b
[2][1],
862 b
[0][1] + b
[1][0], b
[0][2] + b
[2][0], b
[0][0] );
865 ConicCartesianData::ConicCartesianData()
869 bool operator==( const ConicPolarData
& lhs
, const ConicPolarData
& rhs
)
871 return lhs
.focus1
== rhs
.focus1
&&
872 lhs
.pdimen
== rhs
.pdimen
&&
873 lhs
.ecostheta0
== rhs
.ecostheta0
&&
874 lhs
.esintheta0
== rhs
.esintheta0
;
877 ConicCartesianData
ConicCartesianData::invalidData()
879 ConicCartesianData ret
;
880 ret
.coeffs
[0] = double_inf
;
884 bool ConicCartesianData::valid() const
886 return finite( coeffs
[0] );