Skip to content

Commit fb4b5a2

Browse files
jmeidowjmeidow
committed
conics moved to geometry
1 parent 0cf9adf commit fb4b5a2

File tree

4 files changed

+243
-85
lines changed

4 files changed

+243
-85
lines changed

‎src/conics.h‎

Lines changed: 0 additions & 78 deletions
This file was deleted.

‎src/geometry/conics.h‎

Lines changed: 236 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,236 @@
1+
/*
2+
* This file is part of the GreasePad distribution (https://github.com/FraunhoferIOSB/GreasePad).
3+
* Copyright (c) 2022-2026 Jochen Meidow, Fraunhofer IOSB
4+
*
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 3 of the License, or
8+
* (at your option) any later version.
9+
*
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.
14+
*
15+
* You should have received a copy of the GNU General Public License
16+
* along with this program. If not, see <https://www.gnu.org/licenses/>.
17+
*/
18+
19+
#ifndef CONICS_H
20+
#define CONICS_H
21+
22+
#include <Eigen/Core>
23+
#include <Eigen/Eigenvalues>
24+
25+
#include <cassert>
26+
#include <cmath>
27+
#include <utility>
28+
29+
#include <QDebug>
30+
31+
#include "qassert.h"
32+
#include "qlogging.h"
33+
34+
#include "matfun.h"
35+
36+
37+
//! conics, rotations, bounding boxes, cross product
38+
namespace Geometry {
39+
40+
using Eigen::Matrix3d;
41+
using Eigen::Vector3d;
42+
using Eigen::Matrix2d;
43+
using Eigen::Vector2d;
44+
45+
using Matfun::skew;
46+
47+
48+
//! Base class for conics
49+
class ConicBase
50+
{
51+
public:
52+
explicit ConicBase(Matrix3d other) //!< Value constructor
53+
: CC(std::move(other))
54+
{
55+
constexpr double T_sym = 1e-6;
56+
assert( ( CC -CC.adjoint() ).norm() < T_sym );
57+
}
58+
59+
[[nodiscard]] Matrix3d C() const { return CC;} //!< getter
60+
[[nodiscard]] bool isCentral() const //!< Check if conic has a central point
61+
{
62+
return CC.topLeftCorner(2,2).determinant() !=0.;
63+
}
64+
65+
66+
[[nodiscard]] std::pair<Vector3d,Vector3d> intersect( const Vector3d & l ) const //!< Two intersection points with a straight line
67+
{
68+
Matrix3d const MM = skew(l);
69+
Matrix3d BB = MM.adjoint()*CC*MM;
70+
71+
int idx = 0; // [den,idx] = max( abs(l) );
72+
double const den = l.array().abs().maxCoeff(&idx);
73+
74+
// minors ...............................................
75+
double alpha = 0;
76+
switch (idx) {
77+
case 0:
78+
alpha = BB(1,1)*BB(2,2) -BB(2,1)*BB(1,2);
79+
break;
80+
case 1:
81+
alpha = BB(0,0)*BB(2,2) -BB(2,0)*BB(0,2);
82+
break;
83+
case 2:
84+
alpha = BB(0,0)*BB(1,1) -BB(1,0)*BB(0,1);
85+
break;
86+
default:
87+
Q_ASSERT_X( false, "ConicBase::intersect",
88+
"intersection of conic and straight line: index out of range");
89+
}
90+
91+
// intersection points .......................................
92+
Q_ASSERT( alpha <= 0 );
93+
Q_ASSERT( den > 0 );
94+
Matrix3d DD = BB +std::sqrt(-alpha)/den*MM;
95+
int r = 0;
96+
int c = 0;
97+
DD.array().abs().maxCoeff( &r, &c);
98+
Vector3d const p = DD.row(r);
99+
Vector3d const q = DD.col(c);
100+
101+
return {p,q};
102+
}
103+
104+
private:
105+
Matrix3d CC;
106+
};
107+
108+
109+
//! Ellipse
110+
class Ellipse : public ConicBase
111+
{
112+
public:
113+
explicit Ellipse(const Matrix3d &CC ) //!< Value constructor (uncertain point)
114+
: ConicBase(CC)
115+
{
116+
// check if matrix represents an ellipse
117+
assert( C().topLeftCorner(2,2).determinant() > 0.0 ); // PCV Table 5.8
118+
}
119+
[[nodiscard]] std::pair<Eigen::VectorXd,Eigen::VectorXd> poly( int N ) const //!< Get N points on ellipse
120+
{
121+
Matrix2d const Chh = C().topLeftCorner(2, 2);
122+
Vector2d const ch0 = C().topRightCorner(2, 1);
123+
Vector2d x0 = -Chh.ldlt().solve(ch0); // centre point
124+
double const c00q = C().coeff(2, 2) - ch0.dot(Chh.ldlt().solve(ch0));
125+
126+
assert( std::fabs(c00q)>0. );
127+
128+
Eigen::EigenSolver<Matrix2d> const eig(-Chh / c00q, true);
129+
Vector2d ev = eig.eigenvalues().real();
130+
Matrix2d const RR = eig.eigenvectors().real();
131+
132+
if ( ev(0)<0.0 ) {
133+
ev *= -1;
134+
}
135+
136+
const double two_pi = 2*3.14159;
137+
Eigen::VectorXd t = Eigen::VectorXd::LinSpaced( N, 0, two_pi);
138+
139+
Eigen::MatrixXd x(2,N);
140+
x.row(0) = t.array().sin()/std::sqrt(ev(0));
141+
x.row(1) = t.array().cos()/std::sqrt(ev(1));
142+
143+
x = RR*x; // rotation
144+
x.row(0).array() += x0(0); // translation
145+
x.row(1).array() += x0(1);
146+
147+
return { x.row(0), x.row(1)};
148+
}
149+
150+
[[nodiscard]] Vector3d polar( const Vector3d & x ) const { return C()*x; } //!< Compute polar l (straight line) for point x, i.e., l=C*x
151+
};
152+
153+
154+
//! Hyperbola
155+
class Hyperbola : public ConicBase
156+
{
157+
public:
158+
explicit Hyperbola(const Matrix3d &CC) //!< Value constructor (uncertain straight line)
159+
: ConicBase(CC)
160+
{
161+
// check if matrix represents a hyperbola
162+
assert( C().topLeftCorner(2,2).determinant() < 0.0); // PCV Table 5.8
163+
}
164+
165+
[[nodiscard]] Vector3d centerline() const //!< Get straight line
166+
{
167+
const Vector3d x0 = center();
168+
const double phi_ = angle_rad();
169+
const double nx = -std::sin(phi_);
170+
const double ny = std::cos(phi_);
171+
172+
assert( fabs(x0(2)) > 0 );
173+
return { nx, ny, -nx*x0(0)/x0(2) -ny*x0(1)/x0(2) };
174+
}
175+
176+
[[nodiscard]] std::pair<double,double> lengthsSemiAxes() const //!< Get lengths/2 of axes
177+
{
178+
auto ev = eigenvalues();
179+
double const Delta = C().determinant();
180+
double const D = C().topLeftCorner(2, 2).determinant();
181+
assert( -Delta/( ev.first*D ) >= 0. );
182+
assert( +Delta/( ev.second*D) >= 0. );
183+
double const a = std::sqrt(-Delta / (ev.first * D));
184+
double const b = std::sqrt(+Delta / (ev.second * D));
185+
186+
return {a,b};
187+
}
188+
189+
[[nodiscard]] double angle_rad() const //!< Get angle between straight line and x-axis in radians.
190+
{
191+
return atan2( 2*C().coeff(0,1), C().coeff(0,0)-C().coeff(1,1)) / 2;
192+
}
193+
194+
[[nodiscard]] double angle_deg() const //!< Get angle between straight line and x-axis in degerees.
195+
{
196+
constexpr double rho = 180./3.14159;
197+
return rho * angle_rad();
198+
}
199+
200+
[[nodiscard]] Vector3d center() const //!< Get center point of conic
201+
{
202+
Vector3d xh;
203+
if ( isCentral() ) {
204+
Matrix2d const C33 = C().topLeftCorner(2, 2);
205+
Vector2d const ch0 = C().topRightCorner(2, 1);
206+
Vector2d x0 = -C33.ldlt().solve(ch0);
207+
xh << x0(0), x0(1), 1.0;
208+
}
209+
else {
210+
qDebug() << "! no central conic"; // TODO(meijoc)
211+
xh << 0,0,0;
212+
}
213+
return xh;
214+
}
215+
216+
217+
private:
218+
[[nodiscard]] std::pair<double,double> eigenvalues() const
219+
{
220+
double const p = -C().topLeftCorner(2, 2).trace();
221+
double const q = C().topLeftCorner(2, 2).determinant();
222+
223+
double const ev0 = -p / 2 - std::sqrt(p * p / 4 - q);
224+
double const ev1 = -p / 2 + std::sqrt(p * p / 4 - q);
225+
226+
return {ev0,ev1};
227+
}
228+
229+
};
230+
231+
232+
233+
} // namespace Geometry
234+
235+
236+
#endif // CONICS_H

‎src/greasepad.pro‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,6 @@ contains(QMAKE_CC, gcc) {
2727

2828
SOURCES += \
2929
adjustment.cpp \
30-
conics.cpp \
3130
constraints.cpp \
3231
main.cpp \
3332
mainwindow.cpp \
@@ -49,10 +48,10 @@ SOURCES += \
4948
HEADERS += \
5049
adjustment.h \
5150
commands.h \
52-
conics.h \
5351
conncomp.h \
5452
constraints.h \
5553
geometry/aabb.h \
54+
geometry/conics.h \
5655
geometry/minrot.h \
5756
global.h \
5857
mainwindow.h \

‎src/qsegment.cpp‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/*
22
* This file is part of the GreasePad distribution (https://github.com/FraunhoferIOSB/GreasePad).
3-
* Copyright (c) 2022-2025 Jochen Meidow, Fraunhofer IOSB
3+
* Copyright (c) 2022-2026 Jochen Meidow, Fraunhofer IOSB
44
*
55
* This program is free software: you can redistribute it and/or modify
66
* it under the terms of the GNU General Public License as published by
@@ -16,7 +16,8 @@
1616
* along with this program. If not, see <https://www.gnu.org/licenses/>.
1717
*/
1818

19-
#include "conics.h"
19+
20+
#include "geometry/conics.h"
2021
#include "matfun.h"
2122
#include "qsegment.h"
2223
#include "upoint.h"
@@ -236,11 +237,11 @@ void QSegment::setShape( const uPoint &ux,
236237
// two ellipses ...............................................
237238
uPoint ux2 = x2.euclidean();
238239
Matrix3d CC = cof3( k2*ux2.Cov() -ux2.v()*ux2.v().adjoint() );
239-
Conic::Ellipse const ell_x(CC);
240+
Geometry::Ellipse const ell_x(CC);
240241

241242
ux2 = y2.euclidean();
242243
CC = cof3( k2*ux2.Cov() -ux2.v()*ux2.v().adjoint() );
243-
Conic::Ellipse const ell_y(CC);
244+
Geometry::Ellipse const ell_y(CC);
244245

245246
ellipse_.first = toPoly( ell_x.poly( nSupport ) );
246247
ellipse_.second = toPoly( ell_y.poly( nSupport ) );
@@ -251,7 +252,7 @@ void QSegment::setShape( const uPoint &ux,
251252
// uStraightLine ul = uStraightLine::cross(x2,y2);
252253
ul = ul.euclidean();
253254
CC = k2*ul.Cov() -ul.v()*ul.v().adjoint();
254-
Conic::Hyperbola const hyp( CC );
255+
Geometry::Hyperbola const hyp( CC );
255256

256257
// Two polar lines ............................................
257258
Vector3d const lx = ell_y.polar(x2.v()).normalized();

0 commit comments

Comments
 (0)