|
1
|
void QgsGeometryUtils::circleCenterRadius( const QgsPointV2& pt1, const QgsPointV2& pt2, const QgsPointV2& pt3, double& radius, double& centerX, double& centerY )
|
|
2
|
{
|
|
3
|
double dx21, dy21, dx31, dy31, h21, h31, d;
|
|
4
|
|
|
5
|
//closed circle
|
|
6
|
if ( qgsDoubleNear( pt1.x(), pt3.x() ) && qgsDoubleNear( pt1.y(), pt3.y() ) )
|
|
7
|
{
|
|
8
|
centerX = pt2.x();
|
|
9
|
centerY = pt2.y();
|
|
10
|
radius = sqrt( pow( pt2.x() - pt1.x(), 2.0 ) + pow( pt2.y() - pt1.y(), 2.0 ) );
|
|
11
|
return;
|
|
12
|
}
|
|
13
|
|
|
14
|
// Using cartesian circumcenter eguations from page https://en.wikipedia.org/wiki/Circumscribed_circle
|
|
15
|
dx21 = pt2.x() - pt1.x();
|
|
16
|
dy21 = pt2.y() - pt1.y();
|
|
17
|
dx31 = pt3.x() - pt1.x();
|
|
18
|
dy31 = pt3.y() - pt1.y();
|
|
19
|
|
|
20
|
h21 = pow( dx21, 2.0 ) + pow( dy21, 2.0 );
|
|
21
|
h31 = pow( dx31, 2.0 ) + pow( dy31, 2.0 );
|
|
22
|
|
|
23
|
// 2*Cross product, d<0 means clockwise and d>0 counterclockwise sweeping angle
|
|
24
|
d = 2 * (dx21 * dy31 - dx31 * dy21);
|
|
25
|
|
|
26
|
// Check colinearity, Cross product = 0
|
|
27
|
if ( qgsDoubleNear( fabs( d ), 0.0, 0.00000000001 ) )
|
|
28
|
{
|
|
29
|
radius = -1.0;
|
|
30
|
return;
|
|
31
|
}
|
|
32
|
|
|
33
|
// Calculate centroid coordinates and radius
|
|
34
|
centerX = pt1.x() + ( h21 * dy31 - h31 * dy21 ) / d;
|
|
35
|
centerY = pt1.y() - ( h21 * dx31 - h31 * dx21 ) / d;
|
|
36
|
radius = sqrt( pow( centerX - pt1.x(), 2.0 ) + pow( centerY - pt1.y(), 2.0 ) );
|
|
37
|
}
|