[mapguide-commits] r10235 - in branches/4.0/MgDev: Common/Geometry/Spatial Server/src/UnitTesting
svn_mapguide at osgeo.org
svn_mapguide at osgeo.org
Fri Oct 30 04:56:40 PDT 2026
Author: jng
Date: 2026-10-30 04:56:40 -0700 (Fri, 30 Oct 2026)
New Revision: 10235
Modified:
branches/4.0/MgDev/Common/Geometry/Spatial/SpatialUtilityCircularArc.cpp
branches/4.0/MgDev/Common/Geometry/Spatial/SpatialUtilityCircularArc.h
branches/4.0/MgDev/Server/src/UnitTesting/TestGeometry.cpp
Log:
Compute the circular arc centre relative to the start position
MgSpatialUtilityCircularArc::ComputeCenterFromThreePositions() solved the
perpendicular-bisector system in absolute coordinates, so the right-hand side
of the first two equations was (||V2||^2 - ||V1||^2) / 2: a difference of two
numbers of order 1e12 whose real value is of the order of the arc's extent
squared. That costs most of the significant digits as soon as the arc is small
compared with the coordinates it sits at. The arc segment reported in #2257 -
three nearly collinear positions at coordinates of order 1e6 and only ~3e-3
apart, so a circle of radius ~1.199 - came out with its centre at (742508.013,
220472.437) and a radius of ~5.315, more than 2 units away from where it
belongs, and every interpolated position landed on that wrong circle instead of
on the arc.
Changes:
* Common/Geometry/Spatial/SpatialUtilityCircularArc.cpp: solve for the centre
relative to the start position (U = V - V1) rather than in absolute
coordinates. The matrix is already built from the position offsets, so only
the right-hand side changes - ||V2 - V1||^2 / 2 and ||V3 - V1||^2 / 2, with
Dxyz identically zero in that frame - and the solved offsets are translated
back by the start position at the end. Both sides of the equations are then
of the order of the arc's own extent. The LU decomposition, the
collinearity check and everything downstream of m_center are unchanged, as
is the circle case that does not reach this function. The same cancellation
also cost a radius-10 arc at 700000 about 6.8e-6 in its centre, which is
now exact to rounding.
* Common/Geometry/Spatial/SpatialUtilityCircularArc.h: export the constructor
and Linearize with MG_GEOMETRY_API. MgGeometry is a DLL on Windows and has
no module definition file, so without the export the import library has no
entry for them and MgUnitTesting, which lives in another DLL, would fail to
link; the macro expands to nothing elsewhere.
* Server/src/UnitTesting/TestGeometry.cpp: add "ArcSegment (ticket 2257)",
built from the three positions of the curve polygon attached to the ticket,
and "ArcSegment 3D (tilted circle)", a radius-5 circle tilted 30 degrees
out of the XY plane at (1e6, 2e6, 300), which covers the 3D branch the
change touches as well. Each computes the circle through its three
positions independently of the code under test and requires the centre and
radius to agree; the ticket case also requires every linearized position to
lie on that circle and the linearization to run from the start to the end
position.
The Geometry suite passes in Debug and Release (ubuntu22, x64) with the fix -
"All tests passed (300 assertions in 39 test cases)", no LeakSanitizer records
- and fails without it on exactly the two new cases: the #2257 centre is 2.17
units out in X, and the 3D one 1.97e-6 out in Y.
Fixes #2257
Modified: branches/4.0/MgDev/Common/Geometry/Spatial/SpatialUtilityCircularArc.cpp
===================================================================
--- branches/4.0/MgDev/Common/Geometry/Spatial/SpatialUtilityCircularArc.cpp 2026-10-30 11:08:53 UTC (rev 10234)
+++ branches/4.0/MgDev/Common/Geometry/Spatial/SpatialUtilityCircularArc.cpp 2026-10-30 11:56:40 UTC (rev 10235)
@@ -263,17 +263,17 @@
* given points *
* *
* Denote V, V1, V2, V3 the vectors of center, start, mid, end *
-* respectively, *
-* then V = (x,y,z) satisfies the following equtions: *
-* ( V - ( V1 + V2 ) / 2 ) * ( V2 - V1 ) = 0 *
-* ( V - ( V1 + V3 ) / 2 ) * ( V3 - V1 ) = 0 *
-* ( V - V1 ) * ( ( V2 - V1 ) X ( V3 - V1 ) ) = 0 *
+* respectively, and U = V - V1 the center relative to the start, *
+* then U = (x,y,z) satisfies the following equations: *
+* ( U - ( V2 - V1 ) / 2 ) * ( V2 - V1 ) = 0 *
+* ( U - ( V3 - V1 ) / 2 ) * ( V3 - V1 ) = 0 *
+* U * ( ( V2 - V1 ) X ( V3 - V1 ) ) = 0 *
* *
* i.e., *
* *
-* x21 * x + y21 * y + z21 * z = ( ||V2||^2 - ||V1||^2 ) / 2 *
-* x31 * x + y31 * y + z31 * z = ( ||V3||^2 - ||V1||^2 ) / 2 *
-* Dyz * x + Dzx * y + Dxy * z = Dxyz *
+* x21 * x + y21 * y + z21 * z = ||V2 - V1||^2 / 2 *
+* x31 * x + y31 * y + z31 * z = ||V3 - V1||^2 / 2 *
+* Dyz * x + Dzx * y + Dxy * z = 0 *
* *
* where *
* ( x21, y21, z21 ) = V2 - V1, *
@@ -282,14 +282,17 @@
* Dzx = z21 * x31 - z31 * x21, *
* Dxy = x21 * y31 - x31 * y21, *
* *
-* Dxyz = V1 * ( V2 X V3 ) *
+* ||V2 - V1||^2 = x21*x21 + y21*y21 + z21*z21, *
+* ||V3 - V1||^2 = x31*x31 + y31*y31 + z31*z31 *
* *
-* | x1 y1 z1 | *
-* = | x21 y21 z21 | = x1 * Dyz + y1 * Dzx + z1 * Dxy *
-* | x31 y31 z31 | *
+* The unknown solved for is U, not V: in absolute coordinates the *
+* right-hand side of the first two equations is the difference *
+* ( ||V2||^2 - ||V1||^2 ) / 2, of two numbers that grow with the *
+* square of the coordinates, so for coordinates much larger than the *
+* arc itself the subtraction costs most of the significant digits. *
+* Referred to the start position both sides stay of the order of the *
+* arc's own extent, and Dxyz is identically zero. *
* *
-* ||Vi||^2 = xi*xi + yi*yi + zi*zi *
-* *
************************************************************************/
void MgSpatialUtilityCircularArc::ComputeCenterFromThreePositions()
{
@@ -314,27 +317,12 @@
/* x21 * y31 - x31 * y21 */
a[8] = a[0] * a[4] - a[3] * a[1]; /* Dxy */
- tmp = m_start->GetX() * m_start->GetX() +
- m_start->GetY() * m_start->GetY() +
- (m_hasZ? (m_start->GetZ() * m_start->GetZ()) : 0.0);
+ /* Right-hand sides in the offsets the matrix already holds:
+ ||V2 - V1||^2 / 2 and ||V3 - V1||^2 / 2. Dxyz is zero in this frame. */
+ b[0] = 0.5 * (a[0] * a[0] + a[1] * a[1] + (m_hasZ? (a[2] * a[2]) : 0.0));
+ b[1] = 0.5 * (a[3] * a[3] + a[4] * a[4] + (m_hasZ? (a[5] * a[5]) : 0.0));
+ b[2] = 0.0;
- b[0] = m_mid->GetX() * m_mid->GetX() +
- m_mid->GetY() * m_mid->GetY() +
- (m_hasZ? (m_mid->GetZ() * m_mid->GetZ()) : 0.0);
-
- b[1] = m_end->GetX() * m_end->GetX() +
- m_end->GetY() * m_end->GetY() +
- (m_hasZ? (m_end->GetZ() * m_end->GetZ()) : 0.0);
-
- b[0] -= tmp;
- b[1] -= tmp;
- b[0] *= 0.5;
- b[1] *= 0.5;
- /* Dxyz = x1 * Dyz + y1 * Dzx + z1 * Dxy */
- b[2] = m_start->GetX() * a[6] +
- m_start->GetY() * a[7] +
- (m_hasZ? (m_start->GetZ() * a[8]) : 0.0);
-
/* Decompose a into LU */
if (MgMathUtility::LUDecompose(3, a, eindex ) == 0)
{
@@ -357,12 +345,14 @@
b[1] -= a[3] * b[0];
b[2] -= a[6] * b[0] + a[7] * b[1];
- /* Use UV = b to solve V, the center */
+ /* Use UV = b to solve U, the center relative to the start position */
double zz = m_hasZ? (b[2] / a[8]) : 0.0;
double yy = (b[1] - a[5] * (m_hasZ? zz : 0.0)) / a[4];
double xx = (b[0] - a[1] * yy - a[2] * (m_hasZ? zz : 0.0)) / a[0];
- m_center = factory.CreateCoordinateXYZ(xx, yy, zz);
+ m_center = factory.CreateCoordinateXYZ(m_start->GetX() + xx,
+ m_start->GetY() + yy,
+ m_hasZ? (m_start->GetZ() + zz) : 0.0);
}
Modified: branches/4.0/MgDev/Common/Geometry/Spatial/SpatialUtilityCircularArc.h
===================================================================
--- branches/4.0/MgDev/Common/Geometry/Spatial/SpatialUtilityCircularArc.h 2026-10-30 11:08:53 UTC (rev 10234)
+++ branches/4.0/MgDev/Common/Geometry/Spatial/SpatialUtilityCircularArc.h 2026-10-30 11:56:40 UTC (rev 10235)
@@ -24,6 +24,10 @@
/// \cond INTERNAL
///////////////////////////////////////////////////////////
// Circular arc with all possible parameters computed.
+//
+// The constructor and Linearize are exported because MgUnitTesting, which lives in
+// a different DLL, checks the centre and the linearization directly. MG_GEOMETRY_API
+// is empty off Windows, where every symbol is exported anyway.
struct MgSpatialUtilityCircularArc
{
Ptr<MgCoordinate> m_start;
@@ -42,7 +46,7 @@
/// Constructor. If the input parameters are valid points that do not describe
/// an arc, the constructor leave m_isValidArc==false, rather than throwing
/// an exception.
- MgSpatialUtilityCircularArc(
+ MG_GEOMETRY_API MgSpatialUtilityCircularArc(
MgCoordinate* start,
MgCoordinate* mid,
MgCoordinate* end);
@@ -61,7 +65,7 @@
/// If the arc is not valid, or otherwise cannot be linearized using
/// the given parameters, the arc's known points will simply be appended
/// as if they described a LineString.
- void Linearize(
+ MG_GEOMETRY_API void Linearize(
MgCoordinateCollection * distinctPositions,
double maxSpacing,
double maxOffset);
Modified: branches/4.0/MgDev/Server/src/UnitTesting/TestGeometry.cpp
===================================================================
--- branches/4.0/MgDev/Server/src/UnitTesting/TestGeometry.cpp 2026-10-30 11:08:53 UTC (rev 10234)
+++ branches/4.0/MgDev/Server/src/UnitTesting/TestGeometry.cpp 2026-10-30 11:56:40 UTC (rev 10235)
@@ -18,9 +18,12 @@
#include "MapGuideCommon.h"
#include "TestServiceFactory.h"
#include "FoundationDefs.h"
+#include "Spatial/SpatialUtilityCircularArc.h"
#include "catch.hpp"
#include "CatchHelperMacros.h"
+#include <math.h>
+
const STRING ArbitraryWkt = L"LOCAL_CS [ \"Non-Earth (Mile)\", LOCAL_DATUM [\"Local Datum\", 0], UNIT [\"Mile_US\", 1609.344], AXIS [\"X\", EAST], AXIS[\"Y\", NORTH]]";
const STRING GeographicWkt = L"GEOGCS[\"LL84\",DATUM[\"WGS84\",SPHEROID[\"WGS84\",6378137.000,298.25722293]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]]";
const STRING ProjectedWkt = L"PROJCS[\"GA-W\",GEOGCS[\"LL27\",DATUM[\"NAD27\",SPHEROID[\"CLRK66\",6378206.400,294.97869821]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"false_easting\",500000.000],PARAMETER[\"false_northing\",0.000],PARAMETER[\"scale_factor\",0.999900000000],PARAMETER[\"central_meridian\",-84.16666666666670],PARAMETER[\"latitude_of_origin\",30.00000000000000],UNIT[\"Foot_US\",0.30480060960122]]";
@@ -2382,6 +2385,95 @@
}
}
+TEST_CASE("ArcSegment (ticket 2257)", "[Geometry]")
+{
+ // The three positions of an arc segment from the curve polygon attached to trac ticket
+ // #2257. They are nearly collinear, so the circle through them has a radius of only
+ // ~1.199, and they sit at coordinates of order 1e6 while being only ~3e-3 apart. That
+ // combination makes the centre very sensitive to the precision of the calculation.
+ MgGeometryFactory factory;
+
+ Ptr<MgCoordinate> start = factory.CreateCoordinateXY(742505.210662424, 220467.920336341);
+ Ptr<MgCoordinate> control = factory.CreateCoordinateXY(742505.209330668, 220467.921167299);
+ Ptr<MgCoordinate> end = factory.CreateCoordinateXY(742505.208, 220467.922);
+
+ // Centre and radius of the circle through the three positions, computed here with the
+ // coordinates referred to the start position so that the squares stay small, as an
+ // independent reference.
+ double px = control->GetX() - start->GetX();
+ double py = control->GetY() - start->GetY();
+ double qx = end->GetX() - start->GetX();
+ double qy = end->GetY() - start->GetY();
+ double det = 2.0 * (px * qy - qx * py);
+ REQUIRE(det != 0.0);
+
+ double centerX = start->GetX() + (qy * (px*px + py*py) - py * (qx*qx + qy*qy)) / det;
+ double centerY = start->GetY() + (px * (qx*qx + qy*qy) - qx * (px*px + py*py)) / det;
+ double radius = sqrt((start->GetX() - centerX) * (start->GetX() - centerX) +
+ (start->GetY() - centerY) * (start->GetY() - centerY));
+
+ // Four orders of magnitude above the rounding error of the calculation, and four below
+ // the error made by solving for the centre in absolute coordinates.
+ const double tolerance = 1.0e-6;
+
+ MgSpatialUtilityCircularArc arc(start, control, end);
+ REQUIRE(arc.m_isValidArc);
+ REQUIRE(fabs(arc.m_center->GetX() - centerX) <= tolerance);
+ REQUIRE(fabs(arc.m_center->GetY() - centerY) <= tolerance);
+ REQUIRE(fabs(arc.m_radius - radius) <= tolerance);
+
+ // Linearize with a spacing short enough to interpolate rather than fall back to the
+ // known positions; every interpolated position has to lie on that circle.
+ Ptr<MgCoordinateCollection> positions = new MgCoordinateCollection();
+ arc.Linearize(positions, 5.0e-4, 0.0);
+ REQUIRE(positions->GetCount() > 2);
+
+ for (INT32 i = 0; i < positions->GetCount(); i++)
+ {
+ Ptr<MgCoordinate> position = positions->GetItem(i);
+ double dx = position->GetX() - centerX;
+ double dy = position->GetY() - centerY;
+ REQUIRE(fabs(sqrt(dx*dx + dy*dy) - radius) <= tolerance);
+ }
+
+ // The linearized arc runs from the start position to the end position.
+ Ptr<MgCoordinate> first = positions->GetItem(0);
+ Ptr<MgCoordinate> last = positions->GetItem(positions->GetCount() - 1);
+ REQUIRE(first->GetX() == start->GetX());
+ REQUIRE(first->GetY() == start->GetY());
+ REQUIRE(last->GetX() == end->GetX());
+ REQUIRE(last->GetY() == end->GetY());
+}
+
+TEST_CASE("ArcSegment 3D (tilted circle)", "[Geometry]")
+{
+ // A radius-5 circle tilted 30 degrees out of the XY plane and centred at
+ // (1000000, 2000000, 300), with the three positions at 0, 60 and 120 degrees around the
+ // plane spanned by (1,0,0) and (0, sqrt(3)/2, 1/2). The arc's own extent is a few units
+ // while the coordinates are of order 1e6.
+ MgGeometryFactory factory;
+
+ const double halfRoot3 = 0.86602540378443864676;
+ const double radius = 5.0;
+
+ Ptr<MgCoordinate> start = factory.CreateCoordinateXYZ(1000000.0 + radius,
+ 2000000.0,
+ 300.0);
+ Ptr<MgCoordinate> control = factory.CreateCoordinateXYZ(1000000.0 + 0.5 * radius,
+ 2000000.0 + halfRoot3 * halfRoot3 * radius,
+ 300.0 + halfRoot3 * 0.5 * radius);
+ Ptr<MgCoordinate> end = factory.CreateCoordinateXYZ(1000000.0 - 0.5 * radius,
+ 2000000.0 + halfRoot3 * halfRoot3 * radius,
+ 300.0 + halfRoot3 * 0.5 * radius);
+
+ MgSpatialUtilityCircularArc arc(start, control, end);
+ REQUIRE(arc.m_isValidArc);
+ REQUIRE(fabs(arc.m_center->GetX() - 1000000.0) <= 1.0e-6);
+ REQUIRE(fabs(arc.m_center->GetY() - 2000000.0) <= 1.0e-6);
+ REQUIRE(fabs(arc.m_center->GetZ() - 300.0) <= 1.0e-6);
+ REQUIRE(fabs(arc.m_radius - radius) <= 1.0e-6);
+}
+
TEST_CASE("Simplify with bad params", "[Geometry]")
{
try
@@ -2908,4 +3000,4 @@
{
throw;
}
-}
\ No newline at end of file
+}
More information about the mapguide-commits
mailing list