diff --git a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
index a144a66..fc9de10 100644
--- a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
+++ b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
@@ -20,6 +20,7 @@
 import java.util.Arrays;
 import java.util.Collection;
 import java.util.Collections;
+import java.util.Iterator;
 import java.util.List;
 import java.util.stream.Collectors;
 import java.util.stream.Stream;
@@ -47,6 +48,12 @@
     /** Constant containing the area of half of the spherical space. */
     private static final double HALF_SIZE = PlaneAngleRadians.TWO_PI;
 
+    /** Empirically determined threshold for computing the weighted centroid vector using the
+     * triangle fan approach. Areas with boundary sizes under this value use the triangle fan
+     * method to increase centroid accuracy.
+     */
+    private static final double TRIANGLE_FAN_CENTROID_COMPUTE_THRESHOLD = 1e-2;
+
     /** Construct an instance from its boundaries. Callers are responsible for ensuring
      * that the given path represents the boundary of a convex area. No validation is
      * performed.
@@ -132,34 +139,35 @@
         return weighted == null ? null : Point2S.from(weighted);
     }
 
-    /** Returns the weighted vector for the centroid. This vector is computed by scaling the
-     * pole vector of the great circle of each boundary arc by the size of the arc and summing
-     * the results. By combining the weighted centroid vectors of multiple areas, a single
-     * centroid can be computed for the whole group.
+    /** Return the weighted centroid vector of the area. The returned vector points in the direction of the
+     * centroid point on the surface of the unit sphere with the length of the vector proportional to the
+     * effective mass of the area at the centroid. By adding the weighted centroid vectors of multiple
+     * convex areas, a single centroid can be computed for the combined area.
      * @return weighted centroid vector.
      * @see <a href="https://archive.org/details/centroidinertiat00broc">
      *  <em>The Centroid and Inertia Tensor for a Spherical Triangle</em> - John E. Brock</a>
      */
     Vector3D getWeightedCentroidVector() {
         final List<GreatArc> arcs = getBoundaries();
-        switch (arcs.size()) {
+        final int numBoundaries = arcs.size();
+
+        switch (numBoundaries) {
         case 0:
             // full space; no centroid
             return null;
         case 1:
-            // hemisphere; centroid is the pole of the hemisphere
-            final GreatArc singleArc = arcs.get(0);
-            return singleArc.getCircle().getPole().withNorm(singleArc.getSize());
+            // hemisphere
+            return computeHemisphereWeightedCentroidVector(arcs.get(0));
+        case 2:
+            // lune
+            return computeLuneWeightedCentroidVector(arcs.get(0), arcs.get(1));
         default:
-            // 2 or more sides; use an extension of the approach outlined here:
-            // https://archive.org/details/centroidinertiat00broc
-            // In short, the centroid is the sum of the pole vectors of each side
-            // multiplied by their arc lengths.
-            Vector3D centroid = Vector3D.ZERO;
-            for (final GreatArc arc : getBoundaries()) {
-                centroid = centroid.add(arc.getCircle().getPole().withNorm(arc.getSize()));
+            // triangle or other convex polygon
+            if (getBoundarySize() < TRIANGLE_FAN_CENTROID_COMPUTE_THRESHOLD) {
+                return computeTriangleFanWeightedCentroidVector(arcs);
             }
-            return centroid;
+
+            return computeArcPoleWeightedCentroidVector(arcs);
         }
     }
 
@@ -313,4 +321,131 @@
                 full() :
                 new ConvexArea2S(arcs);
     }
+
+    /** Compute the weighted centroid vector for the hemisphere formed by the given arc.
+     * @param arc arc defining the hemisphere
+     * @return the weighted centroid vector for the hemisphere
+     * @see #getWeightedCentroidVector()
+     */
+    private static Vector3D computeHemisphereWeightedCentroidVector(final GreatArc arc) {
+        return arc.getCircle().getPole().withNorm(HALF_SIZE);
+    }
+
+    /** Compute the weighted centroid vector for the lune formed by the given arcs.
+     * @param a first arc for the lune
+     * @param b second arc for the lune
+     * @return the weighted centroid vector for the lune
+     * @see #getWeightedCentroidVector()
+     */
+    private static Vector3D computeLuneWeightedCentroidVector(final GreatArc a, final GreatArc b) {
+        final Point2S aMid = a.getCentroid();
+        final Point2S bMid = b.getCentroid();
+
+        // compute the centroid vector as the exact center of the lune to avoid inaccurate
+        // results with very small regions
+        final Vector3D.Unit centroid = aMid.slerp(bMid, 0.5).getVector();
+
+        // compute the weight using the reverse of the algorithm from computeArcPoleWeightedCentroidVector()
+        final double weight =
+            (a.getSize() * centroid.dot(a.getCircle().getPole())) +
+            (b.getSize() * centroid.dot(b.getCircle().getPole()));
+
+        return centroid.withNorm(weight);
+    }
+
+    /** Compute the weighted centroid vector for the triangle or polygon formed by the given arcs
+     * by adding together the arc pole vectors multiplied by their respective arc lengths. This
+     * algorithm is described in the paper <a href="https://archive.org/details/centroidinertiat00broc">
+     * <em>The Centroid and Inertia Tensor for a Spherical Triangle</em></a> by John E Brock.
+     *
+     * <p>Note: This algorithm works well in general but is susceptible to floating point errors
+     * on very small areas. In these cases, the computed centroid may not be in the expected location
+     * and may even be outside of the area. The {@link #computeTriangleFanWeightedCentroidVector(List)}
+     * method can produce more accurate results in these cases.</p>
+     * @param arcs boundary arcs for the area
+     * @return the weighted centroid vector for the area
+     * @see #computeTriangleFanWeightedCentroidVector(List)
+     */
+    private static Vector3D computeArcPoleWeightedCentroidVector(final List<GreatArc> arcs) {
+        Vector3D centroid = Vector3D.ZERO;
+
+        Vector3D arcContribution;
+        for (final GreatArc arc : arcs) {
+            arcContribution = arc.getCircle().getPole().withNorm(arc.getSize());
+
+            centroid = centroid.add(arcContribution);
+        }
+
+        return centroid;
+    }
+
+    /** Compute the weighted centroid vector for the triangle or polygon formed by the given arcs
+     * using a triangle fan approach. This method is specifically designed for use with areas of very small size,
+     * where use of the standard algorithm from {@link ##computeArcPoleWeightedCentroidVector(List))} can produce
+     * inaccurate results. The algorithm proceeds as follows:
+     * <ol>
+     *  <li>The polygon is divided into spherical triangles using a triangle fan.</li>
+     *  <li>For each triangle, the vectors of the 3 spherical points are added together to approximate the direction
+     *      of the spherical centroid. This ensures that the computed centroid lies within the area.</li>
+     *  <li>The length of the weighted centroid vector is determined by computing the sum of the contributions that
+     *      each arc in the triangle would make to the centroid using the algorithm from
+     *      {@link ##computeArcPoleWeightedCentroidVector(List)}. This essentially performs part of that algorithm in
+     *      reverse: given a centroid direction, compute the contribution that each arc makes.</li>
+     *  <li>The sum of the weighted centroid vectors for each triangle is computed and returned.</li>
+     * </ol>
+     * @param arcs boundary arcs for the area; must contain at least 3 arcs
+     * @return the weighted centroid vector for the area
+     * @see ##computeArcPoleWeightedCentroidVector(List)
+     */
+    private static Vector3D computeTriangleFanWeightedCentroidVector(final List<GreatArc> arcs) {
+        final Iterator<GreatArc> arcIt = arcs.iterator();
+
+        final Point2S p0 = arcIt.next().getStartPoint();
+        final Vector3D.Unit v0 = p0.getVector();
+
+        Vector3D areaCentroid = Vector3D.ZERO;
+
+        GreatArc arc;
+        Point2S p1;
+        Point2S p2;
+        Vector3D.Unit v1;
+        Vector3D.Unit v2;
+        Vector3D.Unit triangleCentroid;
+        double triangleCentroidLen;
+        while (arcIt.hasNext()) {
+            arc = arcIt.next();
+
+            if (!arc.contains(p0)) {
+                p1 = arc.getStartPoint();
+                p2 = arc.getEndPoint();
+
+                v1 = p1.getVector();
+                v2 = p2.getVector();
+
+                triangleCentroid = v0.add(v1).add(v2).normalize();
+                triangleCentroidLen =
+                        computeArcCentroidContribution(v0, v1, triangleCentroid) +
+                        computeArcCentroidContribution(v1, v2, triangleCentroid) +
+                        computeArcCentroidContribution(v2, v0, triangleCentroid);
+
+                areaCentroid = areaCentroid.add(triangleCentroid.withNorm(triangleCentroidLen));
+            }
+        }
+
+        return areaCentroid;
+    }
+
+    /** Compute the contribution made by a single arc to a weighted centroid vector.
+     * @param a first point in the arc
+     * @param b second point in the arc
+     * @param triangleCentroid the centroid vector for the area
+     * @return the contribution made by the arc {@code ab} to the length of the weighted centroid vector
+     */
+    private static double computeArcCentroidContribution(final Vector3D.Unit a, final Vector3D.Unit b,
+            final Vector3D.Unit triangleCentroid) {
+        final double arcLength = a.angle(b);
+        final Vector3D.Unit planeNormal = a.cross(b).normalize();
+
+        return arcLength * triangleCentroid.dot(planeNormal);
+    }
 }
diff --git a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/GreatCircle.java b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/GreatCircle.java
index c1e830c..4847449 100644
--- a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/GreatCircle.java
+++ b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/GreatCircle.java
@@ -32,7 +32,8 @@
  * intersection of a sphere with a plane that passes through its center. It is
  * the largest diameter circle that can be drawn on the sphere and partitions the
  * sphere into two hemispheres. The vectors {@code u} and {@code v} lie in the great
- * circle plane, while the vector {@code w} (the pole) is perpendicular to it.
+ * circle plane, while the vector {@code w} (the pole) is perpendicular to it. The
+ * pole vector points toward the <em>minus</em> side of the hyperplane.
  *
  * <p>Instances of this class are guaranteed to be immutable.</p>
  * @see GreatCircles
diff --git a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
index 5a9d20e..2b8cf2e 100644
--- a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
+++ b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
@@ -171,16 +171,32 @@
         final DoublePrecisionContext precision = ((GreatArc) getRoot().getCut()).getPrecision();
 
         double sizeSum = 0;
-        Vector3D centroidVector = Vector3D.ZERO;
+        Vector3D centroidVectorSum = Vector3D.ZERO;
+
+        Vector3D centroidVector;
+        double maxCentroidVectorWeightSq = 0.0;
 
         for (final ConvexArea2S area : areas) {
             sizeSum += area.getSize();
-            centroidVector = centroidVector.add(area.getWeightedCentroidVector());
+
+            centroidVector = area.getWeightedCentroidVector();
+            maxCentroidVectorWeightSq = Math.max(maxCentroidVectorWeightSq, centroidVector.normSq());
+
+            centroidVectorSum = centroidVectorSum.add(centroidVector);
         }
 
-        final Point2S centroid = centroidVector.eq(Vector3D.ZERO, precision) ?
-                null :
-                Point2S.from(centroidVector);
+        // Convert the weighted centroid vector to a point on the sphere surface. If the centroid vector
+        // length is less than the max length of the combined convex areas and the vector itself is
+        // equivalent to zero, then we know that there are opposing and approximately equal areas in the
+        // region, resulting in an indeterminate centroid. This would occur, for example, if there were
+        // equal areas around each pole.
+        final Point2S centroid;
+        if (centroidVectorSum.normSq() < maxCentroidVectorWeightSq &&
+                centroidVectorSum.eq(Vector3D.ZERO, precision)) {
+            centroid = null;
+        } else {
+            centroid = Point2S.from(centroidVectorSum);
+        }
 
         return new RegionSizeProperties<>(sizeSum, centroid);
     }
diff --git a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java
index 6683be5..ddcb505 100644
--- a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java
+++ b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java
@@ -538,6 +538,76 @@
     }
 
     @Test
+    public void testGetCentroid_diminishingLunes() {
+        // arrange
+        final double eps = 1e-14;
+        final DoublePrecisionContext precision = new EpsilonDoublePrecisionContext(eps);
+
+        final double centerAz = 1;
+        final double centerPolar = 0.5 * Math.PI;
+        final Point2S center = Point2S.of(centerAz, centerPolar);
+        final Point2S pole = Point2S.PLUS_K;
+
+        final double startOffset = PlaneAngleRadians.PI_OVER_TWO;
+        final double minOffset = 1e-14;
+
+        ConvexArea2S area;
+        Point2S p1;
+        Point2S p2;
+        Point2S centroid;
+        for (double offset = startOffset; offset > minOffset; offset *= 0.5) {
+            p1 = Point2S.of(centerAz - offset, centerPolar);
+            p2 = Point2S.of(centerAz + offset, centerPolar);
+
+            area = ConvexArea2S.fromBounds(
+                    GreatCircles.fromPoints(pole, p1, precision),
+                    GreatCircles.fromPoints(p2, pole, precision));
+
+            // act
+            centroid = area.getCentroid();
+
+            // assert
+            Assert.assertTrue(area.contains(centroid));
+            SphericalTestUtils.assertPointsEq(center, centroid, TEST_EPS);
+        }
+    }
+
+    @Test
+    public void testGetCentroid_diminishingSquares() {
+        // arrange
+        final double eps = 1e-14;
+        final DoublePrecisionContext precision = new EpsilonDoublePrecisionContext(eps);
+
+        final double centerAz = 1;
+        final double centerPolar = 0.5 * Math.PI;
+        final Point2S center = Point2S.of(centerAz, centerPolar);
+
+        final double minOffset = 1e-14;
+
+        ConvexArea2S area;
+        Point2S p1;
+        Point2S p2;
+        Point2S p3;
+        Point2S p4;
+        Point2S centroid;
+        for (double offset = 0.5; offset > minOffset; offset *= 0.5) {
+            p1 = Point2S.of(centerAz, centerPolar - offset);
+            p2 = Point2S.of(centerAz - offset, centerPolar);
+            p3 = Point2S.of(centerAz, centerPolar + offset);
+            p4 = Point2S.of(centerAz + offset, centerPolar);
+
+            area = ConvexArea2S.fromVertexLoop(Arrays.asList(p1, p2, p3, p4), precision);
+
+            // act
+            centroid = area.getCentroid();
+
+            // assert
+            Assert.assertTrue(area.contains(centroid));
+            SphericalTestUtils.assertPointsEq(center, centroid, TEST_EPS);
+        }
+    }
+
+    @Test
     public void testBoundaryStream() {
         // arrange
         final GreatCircle circle = GreatCircles.fromPole(Vector3D.Unit.PLUS_X, TEST_PRECISION);
@@ -823,8 +893,10 @@
             final ConvexArea2S plus = split.getPlus();
             final double plusSize = plus.getSize();
 
-            final Point2S computedCentroid = Point2S.from(minus.getWeightedCentroidVector()
-                    .add(plus.getWeightedCentroidVector()));
+            final Vector3D minusWeightedCentroid = minus.getWeightedCentroidVector();
+            final Vector3D plusWeightedCentroid = plus.getWeightedCentroidVector();
+
+            final Point2S computedCentroid = Point2S.from(minusWeightedCentroid.add(plusWeightedCentroid));
 
             Assert.assertEquals(size, minusSize + plusSize, TEST_EPS);
             SphericalTestUtils.assertPointsEq(centroid, computedCentroid, TEST_EPS);
diff --git a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
index fb0784c..1b8032f 100644
--- a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
+++ b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
@@ -21,6 +21,7 @@
 import java.util.Collections;
 import java.util.List;
 import java.util.stream.Collectors;
+import java.util.stream.Stream;
 
 import org.apache.commons.geometry.core.GeometryTestUtils;
 import org.apache.commons.geometry.core.RegionLocation;
@@ -361,7 +362,11 @@
         SphericalTestUtils.assertPointsEq(Point2S.MINUS_I,
                 tree.project(Point2S.of(PlaneAngleRadians.PI + 0.5, PlaneAngleRadians.PI_OVER_TWO)), TEST_EPS);
 
-        SphericalTestUtils.assertPointsEq(Point2S.PLUS_K, tree.project(tree.getCentroid()), TEST_EPS);
+        final Point2S centroid = tree.getCentroid();
+        SphericalTestUtils.assertPointsEq(Point2S.PLUS_K,
+                tree.project(centroid.slerp(Point2S.PLUS_K, 1e-10)), TEST_EPS);
+        SphericalTestUtils.assertPointsEq(Point2S.PLUS_J,
+                tree.project(centroid.slerp(Point2S.PLUS_J, 1e-10)), TEST_EPS);
     }
 
     @Test
@@ -560,6 +565,42 @@
     }
 
     @Test
+    public void testGeometricProperties_polygonWithHole_small() {
+        // arrange
+        final Point2S center = Point2S.of(0.5, 2);
+
+        final double outerRadius = 1e-5;
+        final double innerRadius = 1e-7;
+
+        final RegionBSPTree2S outer = buildDiamond(center, outerRadius);
+        final RegionBSPTree2S inner = buildDiamond(center, innerRadius);
+
+        // rotate the inner diamond a quarter turn to become a square
+        inner.transform(Transform2S.createRotation(center, 0.25 * Math.PI));
+
+        // act
+        final RegionBSPTree2S tree = RegionBSPTree2S.empty();
+        tree.difference(outer, inner);
+
+        // assert
+
+        // use Euclidean approximations of the area and boundary size since those will be more accurate
+        // at these sizes
+        final double area = (2 * outerRadius * outerRadius) - (2 * innerRadius * innerRadius);
+        Assert.assertEquals(area, tree.getSize(), TEST_EPS);
+
+        final double outerSideLength = Math.hypot(outerRadius, outerRadius);
+        final double innerSideLength = Math.hypot(innerRadius, innerRadius);
+        final double boundarySize = 4 * (outerSideLength + innerSideLength);
+        Assert.assertEquals(boundarySize, tree.getBoundarySize(), TEST_EPS);
+
+        SphericalTestUtils.assertPointsEq(center, tree.getCentroid(), TEST_EPS);
+        checkCentroidConsistency(tree);
+
+        SphericalTestUtils.checkClassify(tree, RegionLocation.OUTSIDE, center);
+    }
+
+    @Test
     public void testGeometricProperties_polygonWithHole_complex() {
         // arrange
         final Point2S center = Point2S.of(0.5, 2);
@@ -602,6 +643,59 @@
     }
 
     @Test
+    public void testGeometricProperties_smallRightTriangle() {
+        // arrange
+        final double azOffset = 1e-5;
+        final double polarOffset = 1e-6;
+
+        final double minAz = 0;
+        final double maxAz = minAz + azOffset;
+        final double maxPolar = PlaneAngleRadians.PI_OVER_TWO;
+        final double minPolar = maxPolar - polarOffset;
+
+        final Point2S p0 = Point2S.of(minAz, maxPolar);
+        final Point2S p1 = Point2S.of(maxAz, maxPolar);
+        final Point2S p2 = Point2S.of(maxAz, minPolar);
+
+        // act
+        final RegionBSPTree2S tree = GreatArcPath.fromVertexLoop(Arrays.asList(p0, p1, p2), TEST_PRECISION)
+                .toTree();
+
+        // assert
+
+        // use Euclidean approximations of the area and boundary size since those will be more accurate
+        // at these sizes
+        final double expectedArea = 0.5 * azOffset * polarOffset;
+        Assert.assertEquals(expectedArea, tree.getSize(), TEST_EPS);
+
+        final double expectedBoundarySize = azOffset + polarOffset + Math.hypot(azOffset, polarOffset);
+        Assert.assertEquals(expectedBoundarySize, tree.getBoundarySize(), TEST_EPS);
+
+        Assert.assertTrue(tree.contains(tree.getCentroid()));
+        checkCentroidConsistency(tree);
+
+        SphericalTestUtils.checkClassify(tree, RegionLocation.INSIDE,
+                tree.getCentroid(),
+                Point2S.of(minAz + (0.75 * azOffset), minPolar + (0.75 * polarOffset)));
+
+        SphericalTestUtils.checkClassify(tree, RegionLocation.BOUNDARY,
+                p0, p1, p2, p0.slerp(p1, 0.5), p1.slerp(p2, 0.5), p2.slerp(p0, 0.5));
+
+        final double midAz = minAz + (0.5 * azOffset);
+        final double pastMinAz = minAz - azOffset;
+        final double pastMaxAz = maxAz + azOffset;
+
+        final double midPolar = minPolar + (0.5 * polarOffset);
+        final double pastMinPolar = minPolar - polarOffset;
+        final double pastMaxPolar = maxPolar + polarOffset;
+
+        SphericalTestUtils.checkClassify(tree, RegionLocation.OUTSIDE,
+                tree.getCentroid().antipodal(),
+                Point2S.of(pastMinAz, midPolar), Point2S.of(pastMaxAz, midPolar),
+                Point2S.of(midAz, pastMinPolar), Point2S.of(midAz, pastMaxPolar));
+    }
+
+    @Test
     public void testGeometricProperties_equalAndOppositeRegions() {
         // arrange
         final Point2S center = Point2S.PLUS_I;
@@ -761,7 +855,7 @@
     @Test
     public void testGeographicMap() {
         // arrange
-        final RegionBSPTree2S continental = latLongToTree(new double[][] {
+        final RegionBSPTree2S continental = latLongToTree(TEST_PRECISION, new double[][] {
                 {51.14850,  2.51357}, {50.94660,  1.63900}, {50.12717,  1.33876}, {49.34737, -0.98946},
                 {49.77634, -1.93349}, {48.64442, -1.61651}, {48.90169, -3.29581}, {48.68416, -4.59234},
                 {47.95495, -4.49155}, {47.57032, -2.96327}, {46.01491, -1.19379}, {44.02261, -1.38422},
@@ -771,7 +865,7 @@
                 {46.27298,  6.02260}, {46.72577,  6.03738}, {47.62058,  7.46675}, {49.01778,  8.09927},
                 {49.20195,  6.65822}, {49.44266,  5.89775}, {49.98537,  4.79922}
             });
-        final RegionBSPTree2S corsica = latLongToTree(new double[][] {
+        final RegionBSPTree2S corsica = latLongToTree(TEST_PRECISION, new double[][] {
                 {42.15249,  9.56001}, {43.00998,  9.39000}, {42.62812,  8.74600}, {42.25651,  8.54421},
                 {41.58361,  8.77572}, {41.38000,  9.22975}
             });
@@ -797,11 +891,12 @@
         final int numPts = 200;
 
         // counterclockwise
-        final RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false);
-        SphericalTestUtils.assertPointsEq(center, ccw.getCentroid(), CENTROID_EPS);
+        final RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false, TEST_PRECISION);
+        SphericalTestUtils.assertPointsEq(center, ccw.getCentroid(), TEST_EPS);
 
         // clockwise; centroid should just be antipodal for the circle center
-        final RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true);
+        final RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true, TEST_PRECISION);
+
         SphericalTestUtils.assertPointsEq(center.antipodal(), cw.getCentroid(), CENTROID_EPS);
     }
 
@@ -815,10 +910,10 @@
         final double ccwArea = 4.0 * PlaneAngleRadians.PI * Math.pow(Math.sin(radius / 2.0), 2.0);
         final double cwArea = 4.0 * PlaneAngleRadians.PI - ccwArea;
 
-        final RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false);
+        final RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false, TEST_PRECISION);
         Assert.assertEquals("Counterclockwise size", ccwArea, ccw.getSize(), TEST_EPS);
 
-        final RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true);
+        final RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true, TEST_PRECISION);
         Assert.assertEquals("Clockwise size", cwArea, cw.getSize(), TEST_EPS);
     }
 
@@ -831,13 +926,59 @@
         // boundary size is independent from winding
         final double boundary = PlaneAngleRadians.TWO_PI * Math.sin(radius);
 
-        final RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false);
+        final RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false, TEST_PRECISION);
         Assert.assertEquals("Counterclockwise boundary size", boundary, ccw.getBoundarySize(), 1.0e-7);
 
-        final RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true);
+        final RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true, TEST_PRECISION);
         Assert.assertEquals("Clockwise boundary size", boundary, cw.getBoundarySize(), 1.0e-7);
     }
 
+    @Test
+    public void testSmallCircleToPolygon() {
+        // arrange
+        final double radius = 5.0e-8;
+        final Point2S center = Point2S.of(0.5, 1.5);
+        final int numPts = 100;
+
+        // act
+        final RegionBSPTree2S circle = circleToPolygon(center, radius, numPts, false, TEST_PRECISION);
+
+        // assert
+        // https://en.wikipedia.org/wiki/Spherical_cap
+        final double area = 4.0 * PlaneAngleRadians.PI * Math.pow(Math.sin(radius / 2.0), 2.0);
+        final double boundary = PlaneAngleRadians.TWO_PI * Math.sin(radius);
+
+        SphericalTestUtils.assertPointsEq(center, circle.getCentroid(), TEST_EPS);
+        Assert.assertEquals(area, circle.getSize(), TEST_EPS);
+        Assert.assertEquals(boundary, circle.getBoundarySize(), TEST_EPS);
+    }
+
+    @Test
+    public void testSmallGeographicalRectangle() {
+        // arrange
+        final double[][] vertices = {
+            {42.656216727628696, -70.61919768884546},
+            {42.65612858998112, -70.61938607250165},
+            {42.65579098923594, -70.61909615581666},
+            {42.655879126692355, -70.61890777301083}
+        };
+
+        // act
+        final RegionBSPTree2S rectangle = latLongToTree(TEST_PRECISION, vertices);
+
+        // assert
+        // approximate the centroid as average of vertices
+        final double avgLat = Stream.of(vertices).mapToDouble(v -> v[0]).average().getAsDouble();
+        final double avgLon = Stream.of(vertices).mapToDouble(v -> v[1]).average().getAsDouble();
+        final Point2S expectedCentroid = latLongToPoint(avgLat, avgLon);
+
+        SphericalTestUtils.assertPointsEq(expectedCentroid, rectangle.getCentroid(), TEST_EPS);
+
+        // expected results computed using GeographicLib (https://geographiclib.sourceforge.io/)
+        Assert.assertEquals(1.997213869978027E-11, rectangle.getSize(), TEST_EPS);
+        Assert.assertEquals(1.9669710464585642E-5, rectangle.getBoundarySize(), TEST_EPS);
+    }
+
     /**
      * Insert hyperplane convex subsets defining the positive quadrant area.
      * @param tree
@@ -875,8 +1016,8 @@
         }
     }
 
-    private static RegionBSPTree2S latLongToTree(final double[][] points) {
-        final GreatArcPath.Builder pathBuilder = GreatArcPath.builder(TEST_PRECISION);
+    private static RegionBSPTree2S latLongToTree(final DoublePrecisionContext precision, final double[][] points) {
+        final GreatArcPath.Builder pathBuilder = GreatArcPath.builder(precision);
 
         for (int i = 0; i < points.length; ++i) {
             pathBuilder.append(latLongToPoint(points[i][0], points[i][1]));
@@ -972,10 +1113,11 @@
         final double angleB = Math.asin(Math.sin(b) / sinC);
 
         // use Girard's theorem
-        return angleA + angleB  - (0.5 * PlaneAngleRadians.PI);
+        return angleA + angleB - PlaneAngleRadians.PI_OVER_TWO;
     }
 
-    private static RegionBSPTree2S circleToPolygon(final Point2S center, final double radius, final int numPts, final boolean clockwise) {
+    private static RegionBSPTree2S circleToPolygon(final Point2S center, final double radius, final int numPts,
+                                                   final boolean clockwise, final DoublePrecisionContext precision) {
         final List<Point2S> pts = new ArrayList<>(numPts);
 
         // get an arbitrary point on the circle boundary
@@ -990,6 +1132,6 @@
             pts.add(rotate.apply(pts.get(i - 1)));
         }
 
-        return GreatArcPath.fromVertexLoop(pts, TEST_PRECISION).toTree();
+        return GreatArcPath.fromVertexLoop(pts, precision).toTree();
     }
 }