1 /***************************************************************************
2 qgsellipse.cpp
3 --------------
4 begin : March 2017
5 copyright : (C) 2017 by Loîc Bartoletti
6 email : lbartoletti at tuxfamily dot org
7 ***************************************************************************/
8
9 /***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17
18 #include "qgsunittypes.h"
19 #include "qgslinestring.h"
20 #include "qgsellipse.h"
21 #include "qgsgeometryutils.h"
22
23 #include <memory>
24 #include <limits>
25
normalizeAxis()26 void QgsEllipse::normalizeAxis()
27 {
28 mSemiMajorAxis = std::fabs( mSemiMajorAxis );
29 mSemiMinorAxis = std::fabs( mSemiMinorAxis );
30 if ( mSemiMajorAxis < mSemiMinorAxis )
31 {
32 std::swap( mSemiMajorAxis, mSemiMinorAxis );
33 mAzimuth = 180.0 / M_PI *
34 QgsGeometryUtils::normalizedAngle( M_PI / 180.0 * ( mAzimuth + 90 ) );
35 }
36 }
37
QgsEllipse(const QgsPoint & center,const double axis_a,const double axis_b,const double azimuth)38 QgsEllipse::QgsEllipse( const QgsPoint ¢er, const double axis_a, const double axis_b, const double azimuth )
39 : mCenter( center )
40 , mSemiMajorAxis( axis_a )
41 , mSemiMinorAxis( axis_b )
42 , mAzimuth( azimuth )
43 {
44 normalizeAxis();
45 }
46
fromFoci(const QgsPoint & pt1,const QgsPoint & pt2,const QgsPoint & pt3)47 QgsEllipse QgsEllipse::fromFoci( const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3 )
48 {
49 const double dist_p1p2 = pt1.distance( pt2 );
50 const double dist_p1p3 = pt1.distance( pt3 );
51 const double dist_p2p3 = pt2.distance( pt3 );
52
53 const double dist = dist_p1p3 + dist_p2p3;
54 const double azimuth = 180.0 / M_PI * QgsGeometryUtils::lineAngle( pt1.x(), pt1.y(), pt2.x(), pt2.y() );
55 QgsPoint center = QgsGeometryUtils::midpoint( pt1, pt2 );
56
57 const double axis_a = dist / 2.0;
58 const double axis_b = std::sqrt( std::pow( axis_a, 2.0 ) - std::pow( dist_p1p2 / 2.0, 2.0 ) );
59
60 QgsGeometryUtils::transferFirstZOrMValueToPoint( QgsPointSequence() << pt1 << pt2 << pt3, center );
61
62 return QgsEllipse( center, axis_a, axis_b, azimuth );
63 }
64
fromExtent(const QgsPoint & pt1,const QgsPoint & pt2)65 QgsEllipse QgsEllipse::fromExtent( const QgsPoint &pt1, const QgsPoint &pt2 )
66 {
67 QgsPoint center = QgsGeometryUtils::midpoint( pt1, pt2 );
68 const double axis_a = std::fabs( pt2.x() - pt1.x() ) / 2.0;
69 const double axis_b = std::fabs( pt2.y() - pt1.y() ) / 2.0;
70 const double azimuth = 90.0;
71
72 QgsGeometryUtils::transferFirstZOrMValueToPoint( QgsPointSequence() << pt1 << pt2, center );
73
74 return QgsEllipse( center, axis_a, axis_b, azimuth );
75 }
76
fromCenterPoint(const QgsPoint & center,const QgsPoint & pt1)77 QgsEllipse QgsEllipse::fromCenterPoint( const QgsPoint ¢er, const QgsPoint &pt1 )
78 {
79 const double axis_a = std::fabs( pt1.x() - center.x() );
80 const double axis_b = std::fabs( pt1.y() - center.y() );
81 const double azimuth = 90.0;
82
83 QgsPoint centerPt( center );
84 QgsGeometryUtils::transferFirstZOrMValueToPoint( QgsPointSequence() << center << pt1, centerPt );
85
86 return QgsEllipse( centerPt, axis_a, axis_b, azimuth );
87 }
88
fromCenter2Points(const QgsPoint & center,const QgsPoint & pt1,const QgsPoint & pt2)89 QgsEllipse QgsEllipse::fromCenter2Points( const QgsPoint ¢er, const QgsPoint &pt1, const QgsPoint &pt2 )
90 {
91 const double azimuth = 180.0 / M_PI * QgsGeometryUtils::lineAngle( center.x(), center.y(), pt1.x(), pt1.y() );
92 const double axis_a = center.distance( pt1 );
93
94 const double length = pt2.distance( QgsGeometryUtils::projectPointOnSegment( pt2, center, pt1 ) );
95 const QgsPoint pp = center.project( length, 90 + azimuth );
96 const double axis_b = center.distance( pp );
97
98 QgsPoint centerPt( center );
99 QgsGeometryUtils::transferFirstZOrMValueToPoint( QgsPointSequence() << center << pt1 << pt2, centerPt );
100
101 return QgsEllipse( centerPt, axis_a, axis_b, azimuth );
102 }
103
operator ==(const QgsEllipse & elp) const104 bool QgsEllipse::operator ==( const QgsEllipse &elp ) const
105 {
106 return ( ( mCenter == elp.mCenter ) &&
107 qgsDoubleNear( mSemiMajorAxis, elp.mSemiMajorAxis, 1E-8 ) &&
108 qgsDoubleNear( mSemiMinorAxis, elp.mSemiMinorAxis, 1E-8 ) &&
109 qgsDoubleNear( mAzimuth, elp.mAzimuth, 1E-8 )
110 );
111 }
112
operator !=(const QgsEllipse & elp) const113 bool QgsEllipse::operator !=( const QgsEllipse &elp ) const
114 {
115 return !operator==( elp );
116 }
117
isEmpty() const118 bool QgsEllipse::isEmpty() const
119 {
120 return ( qgsDoubleNear( mSemiMajorAxis, 0.0, 1E-8 ) ||
121 qgsDoubleNear( mSemiMinorAxis, 0.0, 1E-8 ) );
122 }
123
setSemiMajorAxis(const double axis_a)124 void QgsEllipse::setSemiMajorAxis( const double axis_a )
125 {
126 mSemiMajorAxis = axis_a;
127 normalizeAxis();
128 }
setSemiMinorAxis(const double axis_b)129 void QgsEllipse::setSemiMinorAxis( const double axis_b )
130 {
131 mSemiMinorAxis = axis_b;
132 normalizeAxis();
133 }
134
setAzimuth(const double azimuth)135 void QgsEllipse::setAzimuth( const double azimuth )
136 {
137 mAzimuth = 180.0 / M_PI *
138 QgsGeometryUtils::normalizedAngle( M_PI / 180.0 * azimuth );
139 }
140
focusDistance() const141 double QgsEllipse::focusDistance() const
142 {
143 return std::sqrt( mSemiMajorAxis * mSemiMajorAxis - mSemiMinorAxis * mSemiMinorAxis );
144 }
145
foci() const146 QVector<QgsPoint> QgsEllipse::foci() const
147 {
148 QVector<QgsPoint> f;
149 const double dist_focus = focusDistance();
150 f.append( mCenter.project( dist_focus, mAzimuth ) );
151 f.append( mCenter.project( -dist_focus, mAzimuth ) );
152
153 return f;
154 }
155
eccentricity() const156 double QgsEllipse::eccentricity() const
157 {
158 if ( isEmpty() )
159 {
160 return std::numeric_limits<double>::quiet_NaN();
161 }
162 return focusDistance() / mSemiMajorAxis;
163 }
164
area() const165 double QgsEllipse::area() const
166 {
167 return M_PI * mSemiMajorAxis * mSemiMinorAxis;
168 }
169
perimeter() const170 double QgsEllipse::perimeter() const
171 {
172 const double a = mSemiMajorAxis;
173 const double b = mSemiMinorAxis;
174 return M_PI * ( 3 * ( a + b ) - std::sqrt( 10 * a * b + 3 * ( a * a + b * b ) ) );
175 }
176
quadrant() const177 QVector<QgsPoint> QgsEllipse::quadrant() const
178 {
179 QVector<QgsPoint> quad;
180 quad.append( mCenter.project( mSemiMajorAxis, mAzimuth ) );
181 quad.append( mCenter.project( mSemiMinorAxis, mAzimuth + 90 ) );
182 quad.append( mCenter.project( -mSemiMajorAxis, mAzimuth ) );
183 quad.append( mCenter.project( -mSemiMinorAxis, mAzimuth + 90 ) );
184
185 return quad;
186 }
187
points(unsigned int segments) const188 QgsPointSequence QgsEllipse::points( unsigned int segments ) const
189 {
190 QgsPointSequence pts;
191
192 if ( segments < 3 )
193 {
194 return pts;
195 }
196
197
198 const QgsWkbTypes::Type pType( mCenter.wkbType() );
199 const double z = mCenter.z();
200 const double m = mCenter.m();
201
202 QVector<double> t;
203 t.reserve( segments );
204 const double azimuth = std::atan2( quadrant().at( 0 ).y() - mCenter.y(), quadrant().at( 0 ).x() - mCenter.x() );
205 for ( unsigned int i = 0; i < segments; ++i )
206 {
207 t.append( 2 * M_PI - ( ( 2 * M_PI ) / segments * i ) ); // Since the algorithm used rotates in the trigonometric direction (counterclockwise)
208 }
209
210 for ( QVector<double>::const_iterator it = t.constBegin(); it != t.constEnd(); ++it )
211 {
212 const double x = mCenter.x() +
213 mSemiMajorAxis * std::cos( *it ) * std::cos( azimuth ) -
214 mSemiMinorAxis * std::sin( *it ) * std::sin( azimuth );
215 const double y = mCenter.y() +
216 mSemiMajorAxis * std::cos( *it ) * std::sin( azimuth ) +
217 mSemiMinorAxis * std::sin( *it ) * std::cos( azimuth );
218 pts.push_back( QgsPoint( pType, x, y, z, m ) );
219 }
220
221 return pts;
222 }
223
toPolygon(unsigned int segments) const224 QgsPolygon *QgsEllipse::toPolygon( unsigned int segments ) const
225 {
226 std::unique_ptr<QgsPolygon> polygon( new QgsPolygon() );
227 if ( segments < 3 )
228 {
229 return polygon.release();
230 }
231
232 polygon->setExteriorRing( toLineString( segments ) );
233
234 return polygon.release();
235 }
236
toLineString(unsigned int segments) const237 QgsLineString *QgsEllipse::toLineString( unsigned int segments ) const
238 {
239 std::unique_ptr<QgsLineString> ext( new QgsLineString() );
240 if ( segments < 3 )
241 {
242 return ext.release();
243 }
244
245 QgsPointSequence pts;
246 pts = points( segments );
247 pts.append( pts.at( 0 ) ); // close linestring
248
249 ext->setPoints( pts );
250
251 return ext.release();
252 }
253
boundingBox() const254 QgsRectangle QgsEllipse::boundingBox() const
255 {
256 if ( isEmpty() )
257 {
258 return QgsRectangle();
259 }
260
261 const double angle = mAzimuth * M_PI / 180.0;
262
263 const double ux = mSemiMajorAxis * std::cos( angle );
264 const double uy = mSemiMinorAxis * std::sin( angle );
265 const double vx = mSemiMajorAxis * std::sin( angle );
266 const double vy = mSemiMinorAxis * std::cos( angle );
267
268 const double halfHeight = std::sqrt( ux * ux + uy * uy );
269 const double halfWidth = std::sqrt( vx * vx + vy * vy );
270
271 const QgsPointXY p1( mCenter.x() - halfWidth, mCenter.y() - halfHeight );
272 const QgsPointXY p2( mCenter.x() + halfWidth, mCenter.y() + halfHeight );
273
274 return QgsRectangle( p1, p2 );
275 }
276
toString(int pointPrecision,int axisPrecision,int azimuthPrecision) const277 QString QgsEllipse::toString( int pointPrecision, int axisPrecision, int azimuthPrecision ) const
278 {
279 QString rep;
280 if ( isEmpty() )
281 rep = QStringLiteral( "Empty" );
282 else
283 rep = QStringLiteral( "Ellipse (Center: %1, Semi-Major Axis: %2, Semi-Minor Axis: %3, Azimuth: %4)" )
284 .arg( mCenter.asWkt( pointPrecision ), 0, 's' )
285 .arg( qgsDoubleToString( mSemiMajorAxis, axisPrecision ), 0, 'f' )
286 .arg( qgsDoubleToString( mSemiMinorAxis, axisPrecision ), 0, 'f' )
287 .arg( qgsDoubleToString( mAzimuth, azimuthPrecision ), 0, 'f' );
288
289 return rep;
290 }
291
orientedBoundingBox() const292 QgsPolygon *QgsEllipse::orientedBoundingBox() const
293 {
294 std::unique_ptr<QgsPolygon> ombb( new QgsPolygon() );
295 if ( isEmpty() )
296 {
297 return ombb.release();
298 }
299
300 const QVector<QgsPoint> q = quadrant();
301
302 const QgsPoint p1 = q.at( 0 ).project( mSemiMinorAxis, mAzimuth - 90 );
303 const QgsPoint p2 = q.at( 0 ).project( mSemiMinorAxis, mAzimuth + 90 );
304 const QgsPoint p3 = q.at( 2 ).project( mSemiMinorAxis, mAzimuth + 90 );
305 const QgsPoint p4 = q.at( 2 ).project( mSemiMinorAxis, mAzimuth - 90 );
306
307 QgsLineString *ext = new QgsLineString();
308 ext->setPoints( QgsPointSequence() << p1 << p2 << p3 << p4 );
309
310 ombb->setExteriorRing( ext );
311
312 return ombb.release();
313 }
314