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