qgscircle.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgscircle.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 *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17
18 #include "qgscircle.h"
19 #include "qgslinestring.h"
20 #include "qgsgeometryutils.h"
21 #include "qgstriangle.h"
22
23 #include <memory>
24
26  QgsEllipse( QgsPoint(), 0.0, 0.0, 0.0 )
27 {
28
29 }
30
31 QgsCircle::QgsCircle( const QgsPoint &center, double radius, double azimuth ) :
33 {
34
35 }
36
38 {
40  double azimuth = QgsGeometryUtils::lineAngle( pt1.x(), pt1.y(), pt2.x(), pt2.y() ) * 180.0 / M_PI;
41  double radius = pt1.distance( pt2 ) / 2.0;
42
44
45  return QgsCircle( center, radius, azimuth );
46 }
47
48 static bool isPerpendicular( const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double epsilon )
49 {
50  // check the given point are perpendicular to x or y axis
51
52  double yDelta_a = pt2.y() - pt1.y();
53  double xDelta_a = pt2.x() - pt1.x();
54  double yDelta_b = pt3.y() - pt2.y();
55  double xDelta_b = pt3.x() - pt2.x();
56
57  if ( ( std::fabs( xDelta_a ) <= epsilon ) && ( std::fabs( yDelta_b ) <= epsilon ) )
58  {
59  return false;
60  }
61
62  if ( std::fabs( yDelta_a ) <= epsilon )
63  {
64  return true;
65  }
66  else if ( std::fabs( yDelta_b ) <= epsilon )
67  {
68  return true;
69  }
70  else if ( std::fabs( xDelta_a ) <= epsilon )
71  {
72  return true;
73  }
74  else if ( std::fabs( xDelta_b ) <= epsilon )
75  {
76  return true;
77  }
78
79  return false;
80
81 }
82
83 QgsCircle QgsCircle::from3Points( const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double epsilon )
84 {
85  QgsPoint p1, p2, p3;
86
87  if ( !isPerpendicular( pt1, pt2, pt3, epsilon ) )
88  {
89  p1 = pt1;
90  p2 = pt2;
91  p3 = pt3;
92  }
93  else if ( !isPerpendicular( pt1, pt3, pt2, epsilon ) )
94  {
95  p1 = pt1;
96  p2 = pt3;
97  p3 = pt2;
98  }
99  else if ( !isPerpendicular( pt2, pt1, pt3, epsilon ) )
100  {
101  p1 = pt2;
102  p2 = pt1;
103  p3 = pt3;
104  }
105  else if ( !isPerpendicular( pt2, pt3, pt1, epsilon ) )
106  {
107  p1 = pt2;
108  p2 = pt3;
109  p3 = pt1;
110  }
111  else if ( !isPerpendicular( pt3, pt2, pt1, epsilon ) )
112  {
113  p1 = pt3;
114  p2 = pt2;
115  p3 = pt1;
116  }
117  else if ( !isPerpendicular( pt3, pt1, pt2, epsilon ) )
118  {
119  p1 = pt3;
120  p2 = pt1;
121  p3 = pt2;
122  }
123  else
124  {
125  return QgsCircle();
126  }
129  // Paul Bourke's algorithm
130  double yDelta_a = p2.y() - p1.y();
131  double xDelta_a = p2.x() - p1.x();
132  double yDelta_b = p3.y() - p2.y();
133  double xDelta_b = p3.x() - p2.x();
134
135  if ( qgsDoubleNear( xDelta_a, 0.0, epsilon ) || qgsDoubleNear( xDelta_b, 0.0, epsilon ) )
136  {
137  return QgsCircle();
138  }
139
140  double aSlope = yDelta_a / xDelta_a;
141  double bSlope = yDelta_b / xDelta_b;
142
143  // set z coordinate for center
144  QgsGeometryUtils::setZValueFromPoints( QgsPointSequence() << p1 << p2 << p3, center );
145
146  if ( ( std::fabs( xDelta_a ) <= epsilon ) && ( std::fabs( yDelta_b ) <= epsilon ) )
147  {
148  center.setX( 0.5 * ( p2.x() + p3.x() ) );
149  center.setY( 0.5 * ( p1.y() + p2.y() ) );
150  radius = center.distance( pt1 );
151
152  return QgsCircle( center, radius );
153  }
154
155  if ( std::fabs( aSlope - bSlope ) <= epsilon )
156  {
157  return QgsCircle();
158  }
159
160  center.setX(
161  ( aSlope * bSlope * ( p1.y() - p3.y() ) +
162  bSlope * ( p1.x() + p2.x() ) -
163  aSlope * ( p2.x() + p3.x() ) ) /
164  ( 2.0 * ( bSlope - aSlope ) )
165  );
166  center.setY(
167  -1.0 * ( center.x() - ( p1.x() + p2.x() ) / 2.0 ) /
168  aSlope + ( p1.y() + p2.y() ) / 2.0
169  );
170
171  radius = center.distance( p1 );
172
173  return QgsCircle( center, radius );
174 }
175
177 {
178  return QgsCircle( center, diameter / 2.0, azimuth );
179 }
180
182 {
183  double azimuth = QgsGeometryUtils::lineAngle( center.x(), center.y(), pt1.x(), pt1.y() ) * 180.0 / M_PI;
184
185  QgsPoint centerPt( center );
186  QgsGeometryUtils::setZValueFromPoints( QgsPointSequence() << center << pt1, centerPt );
187
188  return QgsCircle( centerPt, centerPt.distance( pt1 ), azimuth );
189 }
190
191 QgsCircle QgsCircle::from3Tangents( const QgsPoint &pt1_tg1, const QgsPoint &pt2_tg1, const QgsPoint &pt1_tg2, const QgsPoint &pt2_tg2, const QgsPoint &pt1_tg3, const QgsPoint &pt2_tg3, double epsilon )
192 {
193  QgsPoint p1, p2, p3;
194  bool isIntersect = false;
195  QgsGeometryUtils::segmentIntersection( pt1_tg1, pt2_tg1, pt1_tg2, pt2_tg2, p1, isIntersect, epsilon );
196  if ( !isIntersect )
197  return QgsCircle();
198  QgsGeometryUtils::segmentIntersection( pt1_tg1, pt2_tg1, pt1_tg3, pt2_tg3, p2, isIntersect, epsilon );
199  if ( !isIntersect )
200  return QgsCircle();
201  QgsGeometryUtils::segmentIntersection( pt1_tg2, pt2_tg2, pt1_tg3, pt2_tg3, p3, isIntersect, epsilon );
202  if ( !isIntersect )
203  return QgsCircle();
204
205  if ( p1.is3D() )
206  {
207  p1.convertTo( QgsWkbTypes::dropZ( p1.wkbType() ) );
208  }
209
210  if ( p2.is3D() )
211  {
212  p2.convertTo( QgsWkbTypes::dropZ( p2.wkbType() ) );
213  }
214
215  if ( p3.is3D() )
216  {
217  p3.convertTo( QgsWkbTypes::dropZ( p3.wkbType() ) );
218  }
219
220  return QgsTriangle( p1, p2, p3 ).inscribedCircle();
221 }
222
223 QgsCircle QgsCircle::minimalCircleFrom3Points( const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double epsilon )
224 {
225  double l1 = pt2.distance( pt3 );
226  double l2 = pt3.distance( pt1 );
227  double l3 = pt1.distance( pt2 );
228
229  if ( ( l1 * l1 ) - ( l2 * l2 + l3 * l3 ) >= epsilon )
230  return QgsCircle().from2Points( pt2, pt3 );
231  else if ( ( l2 * l2 ) - ( l1 * l1 + l3 * l3 ) >= epsilon )
232  return QgsCircle().from2Points( pt3, pt1 );
233  else if ( ( l3 * l3 ) - ( l1 * l1 + l2 * l2 ) >= epsilon )
234  return QgsCircle().from2Points( pt1, pt2 );
235  else
236  return QgsCircle().from3Points( pt1, pt2, pt3, epsilon );
237 }
238
239 int QgsCircle::intersections( const QgsCircle &other, QgsPoint &intersection1, QgsPoint &intersection2, bool useZ ) const
240 {
241  if ( useZ && mCenter.is3D() && other.center().is3D() && !qgsDoubleNear( mCenter.z(), other.center().z() ) )
242  return 0;
243
244  QgsPointXY int1, int2;
245
248  int1, int2 );
249  if ( res == 0 )
250  return 0;
251
252  intersection1 = QgsPoint( int1 );
253  intersection2 = QgsPoint( int2 );
254  if ( useZ && mCenter.is3D() )
255  {
258  }
259  return res;
260 }
261
262 bool QgsCircle::tangentToPoint( const QgsPointXY &p, QgsPointXY &pt1, QgsPointXY &pt2 ) const
263 {
265 }
266
267 int QgsCircle::outerTangents( const QgsCircle &other, QgsPointXY &line1P1, QgsPointXY &line1P2, QgsPointXY &line2P1, QgsPointXY &line2P2 ) const
268 {
270  QgsPointXY( other.center() ), other.radius(), line1P1, line1P2, line2P1, line2P2 );
271 }
272
274 {
275  double delta_x = std::fabs( pt1.x() - pt2.x() );
276  double delta_y = std::fabs( pt1.x() - pt2.y() );
277  if ( !qgsDoubleNear( delta_x, delta_y ) )
278  {
279  return QgsCircle();
280  }
281
283  QgsGeometryUtils::setZValueFromPoints( QgsPointSequence() << pt1 << pt2, center );
284
285  return QgsCircle( center, delta_x / 2.0, 0 );
286 }
287
288 double QgsCircle::area() const
289 {
290  return M_PI * mSemiMajorAxis * mSemiMajorAxis;
291 }
292
293 double QgsCircle::perimeter() const
294 {
295  return 2.0 * M_PI * mSemiMajorAxis;
296 }
297
299 {
300  mSemiMajorAxis = std::fabs( semiMajorAxis );
302 }
303
305 {
306  mSemiMajorAxis = std::fabs( semiMinorAxis );
308 }
309
311 {
313  quad.append( QgsPoint( mCenter.x(), mCenter.y() + mSemiMajorAxis, mCenter.z() ) );
314  quad.append( QgsPoint( mCenter.x() + mSemiMajorAxis, mCenter.y(), mCenter.z() ) );
315  quad.append( QgsPoint( mCenter.x(), mCenter.y() - mSemiMajorAxis, mCenter.z() ) );
316  quad.append( QgsPoint( mCenter.x() - mSemiMajorAxis, mCenter.y(), mCenter.z() ) );
317
319 }
320
322 {
323  std::unique_ptr<QgsCircularString> circString( new QgsCircularString() );
326  if ( oriented )
327  {
329  }
330  else
331  {
333  }
335  for ( QVector<QgsPoint>::const_iterator it = quad.constBegin(); it != quad.constEnd(); ++it )
336  {
337  points.append( *it );
338  }
339  circString->setPoints( points );
340
341  return circString.release();
342 }
343
344 bool QgsCircle::contains( const QgsPoint &point, double epsilon ) const
345 {
346  return ( mCenter.distance( point ) <= mSemiMajorAxis + epsilon );
347 }
348
350 {
352 }
353
354 QString QgsCircle::toString( int pointPrecision, int radiusPrecision, int azimuthPrecision ) const
355 {
356  QString rep;
357  if ( isEmpty() )
358  rep = QStringLiteral( "Empty" );
359  else
360  rep = QStringLiteral( "Circle (Center: %1, Radius: %2, Azimuth: %3)" )
361  .arg( mCenter.asWkt( pointPrecision ), 0, 's' )
362  .arg( qgsDoubleToString( mSemiMajorAxis, radiusPrecision ), 0, 'f' )
363  .arg( qgsDoubleToString( mAzimuth, azimuthPrecision ), 0, 'f' );
364
365  return rep;
366
367 }
virtual QgsPointSequence points(unsigned int segments=36) const
Returns a list of points with segmentation from segments.
Definition: qgsellipse.cpp:188
Circle geometry type.
Definition: qgscircle.h:43
A rectangle specified with double values.
Definition: qgsrectangle.h:40
double y
Definition: qgspoint.h:42
int intersections(const QgsCircle &other, QgsPoint &intersection1, QgsPoint &intersection2, bool useZ=false) const
Calculates the intersections points between this circle and an other circle.
Definition: qgscircle.cpp:239
static bool segmentIntersection(const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &q1, const QgsPoint &q2, QgsPoint &intersectionPoint, bool &isIntersection, double tolerance=1e-8, bool acceptImproperIntersection=false)
Compute the intersection between two segments.
Adds a z-dimension to the geometry, initialized to a preset value.
Definition: qgspoint.cpp:469
static double lineAngle(double x1, double y1, double x2, double y2)
Calculates the direction of line joining two points in radians, clockwise from the north direction...
QgsCircularString * toCircularString(bool oriented=false) const
Returns a circular string from the circle.
Definition: qgscircle.cpp:321
A class to represent a 2D point.
Definition: qgspointxy.h:43
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:278
QgsWkbTypes::Type wkbType() const
Returns the WKB type of the geometry.
Triangle geometry type.
Definition: qgstriangle.h:33
void setSemiMinorAxis(double semiMinorAxis) override
Inherited method.
Definition: qgscircle.cpp:304
static QgsCircle fromCenterDiameter(const QgsPoint &center, double diameter, double azimuth=0)
Constructs a circle by a center point and a diameter.
Definition: qgscircle.cpp:176
static int circleCircleOuterTangents(const QgsPointXY &center1, double radius1, const QgsPointXY &center2, double radius2, QgsPointXY &line1P1, QgsPointXY &line1P2, QgsPointXY &line2P1, QgsPointXY &line2P2)
Calculates the outer tangent points for two circles, centered at center1 and center2 and with radii o...
QgsCircle inscribedCircle() const
Inscribed circle of the triangle.
static QgsCircle fromExtent(const QgsPoint &pt1, const QgsPoint &pt2)
Constructs a circle by an extent (aka bounding box / QgsRectangle).
Definition: qgscircle.cpp:273
QgsPoint mCenter
Definition: qgsellipse.h:252
static QgsCircle from3Tangents(const QgsPoint &pt1_tg1, const QgsPoint &pt2_tg1, const QgsPoint &pt1_tg2, const QgsPoint &pt2_tg2, const QgsPoint &pt1_tg3, const QgsPoint &pt2_tg3, double epsilon=1E-8)
Constructs a circle by 3 tangents on the circle (aka inscribed circle of a triangle).
Definition: qgscircle.cpp:191
double mSemiMajorAxis
Definition: qgsellipse.h:253
static bool setZValueFromPoints(const QgsPointSequence &points, QgsPoint &point)
A Z dimension is added to point if one of the point in the list points is in 3D.
double mAzimuth
Definition: qgsellipse.h:255
double area() const override
The area of the ellipse.
Definition: qgscircle.cpp:288
bool tangentToPoint(const QgsPointXY &p, QgsPointXY &pt1, QgsPointXY &pt2) const
Calculates the tangent points between this circle and the point p.
Definition: qgscircle.cpp:262
virtual bool isEmpty() const
An ellipse is empty if axes are equal to 0.
Definition: qgsellipse.cpp:118
bool contains(const QgsPoint &point, double epsilon=1E-8) const
Returns true if the circle contains the point.
Definition: qgscircle.cpp:344
static QgsPoint midpoint(const QgsPoint &pt1, const QgsPoint &pt2)
Returns a middle point between points pt1 and pt2.
QString qgsDoubleToString(double a, int precision=17)
Returns a string representation of a double.
Definition: qgis.h:238
QString asWkt(int precision=17) const override
Returns a WKT representation of the geometry.
Definition: qgspoint.cpp:223
void setSemiMajorAxis(double semiMajorAxis) override
Inherited method.
Definition: qgscircle.cpp:298
double azimuth() const
Returns the azimuth.
Definition: qgsellipse.h:139
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
virtual QVector< QgsPoint > quadrant() const
The four quadrants of the ellipse.
Definition: qgsellipse.cpp:177
double mSemiMinorAxis
Definition: qgsellipse.h:254
void setX(double x)
Sets the point&#39;s x-coordinate.
Definition: qgspoint.h:213
double semiMajorAxis() const
Returns the semi-major axis.
Definition: qgsellipse.h:127
double distance(double x, double y) const
Returns the distance between this point and a specified x, y coordinate.
Definition: qgspoint.h:276
void setY(double y)
Sets the point&#39;s y-coordinate.
Definition: qgspoint.h:224
QVector< QgsPoint > QgsPointSequence
double perimeter() const override
The circumference of the ellipse using first approximation of Ramanujan.
Definition: qgscircle.cpp:293
static Type dropZ(Type type)
Drops the z dimension (if present) for a WKB type and returns the new type.
Definition: qgswkbtypes.h:1058
static QgsCircle from3Points(const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double epsilon=1E-8)
Constructs a circle by 3 points on the circle.
Definition: qgscircle.cpp:83
static QgsCircle minimalCircleFrom3Points(const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double epsilon=1E-8)
Constructs the smallest circle from 3 points.
Definition: qgscircle.cpp:223
static bool tangentPointAndCircle(const QgsPointXY &center, double radius, const QgsPointXY &p, QgsPointXY &pt1, QgsPointXY &pt2)
Calculates the tangent points between the circle with the specified center and radius and the point p...
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
bool convertTo(QgsWkbTypes::Type type) override
Converts the geometry to a specified type.
Definition: qgspoint.cpp:538
static int circleCircleIntersections(QgsPointXY center1, double radius1, QgsPointXY center2, double radius2, QgsPointXY &intersection1, QgsPointXY &intersection2)
Calculates the intersections points between the circle with center center1 and radius radius1 and the...
double z
Definition: qgspoint.h:43
Circular string geometry type.
static QgsCircle from2Points(const QgsPoint &pt1, const QgsPoint &pt2)
Constructs a circle by 2 points on the circle.
Definition: qgscircle.cpp:37
QgsRectangle boundingBox() const override
Returns the minimal bounding box for the ellipse.
Definition: qgscircle.cpp:349
int outerTangents(const QgsCircle &other, QgsPointXY &line1P1, QgsPointXY &line1P2, QgsPointXY &line2P1, QgsPointXY &line2P2) const
Calculates the outer tangent points between this circle and an other circle.
Definition: qgscircle.cpp:267
Returns the radius of the circle.
Definition: qgscircle.h:222
Ellipse geometry type.
Definition: qgsellipse.h:39
QgsPoint center() const
Returns the center point.
Definition: qgsellipse.h:121