View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.geometry.spherical.twod;
23  
24  import org.hipparchus.exception.LocalizedCoreFormats;
25  import org.hipparchus.exception.MathIllegalArgumentException;
26  import org.hipparchus.exception.MathIllegalStateException;
27  import org.hipparchus.exception.MathRuntimeException;
28  import org.hipparchus.geometry.LocalizedGeometryFormats;
29  import org.hipparchus.geometry.enclosing.EnclosingBall;
30  import org.hipparchus.geometry.euclidean.threed.Rotation;
31  import org.hipparchus.geometry.euclidean.threed.Vector3D;
32  import org.hipparchus.geometry.partitioning.Region;
33  import org.hipparchus.geometry.partitioning.Region.Location;
34  import org.hipparchus.geometry.partitioning.RegionFactory;
35  import org.hipparchus.geometry.spherical.oned.ArcsSet;
36  import org.hipparchus.geometry.spherical.oned.LimitAngle;
37  import org.hipparchus.geometry.spherical.oned.S1Point;
38  import org.hipparchus.geometry.spherical.oned.Sphere1D;
39  import org.hipparchus.geometry.spherical.oned.SubLimitAngle;
40  import org.hipparchus.random.UnitSphereRandomVectorGenerator;
41  import org.hipparchus.random.Well1024a;
42  import org.hipparchus.util.FastMath;
43  import org.hipparchus.util.MathUtils;
44  import org.junit.jupiter.api.Assertions;
45  import org.junit.jupiter.api.Test;
46  
47  import java.lang.reflect.InvocationTargetException;
48  import java.lang.reflect.Method;
49  import java.util.ArrayList;
50  import java.util.List;
51  import java.util.function.IntPredicate;
52  
53  import static org.junit.jupiter.api.Assertions.assertEquals;
54  import static org.junit.jupiter.api.Assertions.assertNotEquals;
55  import static org.junit.jupiter.api.Assertions.assertNull;
56  import static org.junit.jupiter.api.Assertions.assertSame;
57  import static org.junit.jupiter.api.Assertions.assertTrue;
58  import static org.junit.jupiter.api.Assertions.fail;
59  
60  class SphericalPolygonsSetTest {
61  
62      @Test
63      void testFullSphere() {
64          SphericalPolygonsSet full = new SphericalPolygonsSet(1.0e-10);
65          UnitSphereRandomVectorGenerator random =
66                  new UnitSphereRandomVectorGenerator(3, new Well1024a(0x852fd2a0ed8d2f6dL));
67          for (int i = 0; i < 1000; ++i) {
68              Vector3D v = new Vector3D(random.nextVector());
69              assertEquals(Location.INSIDE, full.checkPoint(new S2Point(v)));
70          }
71          assertEquals(4 * FastMath.PI, new SphericalPolygonsSet(0.01, new S2Point[0]).getSize(), 1.0e-10);
72          assertEquals(0, new SphericalPolygonsSet(0.01, new S2Point[0]).getBoundarySize(), 1.0e-10);
73          assertEquals(0, full.getBoundaryLoops().size());
74          assertTrue(full.getEnclosingCap().getRadius() > 0);
75          assertTrue(Double.isInfinite(full.getEnclosingCap().getRadius()));
76          assertEquals(Location.INSIDE, full.checkPoint(full.getInteriorPoint()));
77      }
78  
79      @Test
80      void testEmpty() {
81          SphericalPolygonsSet empty =
82              (SphericalPolygonsSet) new RegionFactory<Sphere2D, S2Point, Circle, SubCircle>().
83                      getComplement(new SphericalPolygonsSet(1.0e-10));
84          UnitSphereRandomVectorGenerator random =
85                  new UnitSphereRandomVectorGenerator(3, new Well1024a(0x76d9205d6167b6ddL));
86          for (int i = 0; i < 1000; ++i) {
87              Vector3D v = new Vector3D(random.nextVector());
88              assertEquals(Location.OUTSIDE, empty.checkPoint(new S2Point(v)));
89          }
90          assertEquals(0, empty.getSize(), 1.0e-10);
91          assertEquals(0, empty.getBoundarySize(), 1.0e-10);
92          assertEquals(0, empty.getBoundaryLoops().size());
93          assertTrue(empty.getEnclosingCap().getRadius() < 0);
94          assertTrue(Double.isInfinite(empty.getEnclosingCap().getRadius()));
95          assertNull(empty.getInteriorPoint());
96      }
97  
98      @Test
99      void testSouthHemisphere() {
100         double tol = 0.01;
101         double sinTol = FastMath.sin(tol);
102         SphericalPolygonsSet south = new SphericalPolygonsSet(Vector3D.MINUS_K, tol);
103         UnitSphereRandomVectorGenerator random =
104                 new UnitSphereRandomVectorGenerator(3, new Well1024a(0x6b9d4a6ad90d7b0bL));
105         for (int i = 0; i < 1000; ++i) {
106             Vector3D v = new Vector3D(random.nextVector());
107             if (v.getZ() < -sinTol) {
108                 assertEquals(Location.INSIDE, south.checkPoint(new S2Point(v)));
109             } else if (v.getZ() > sinTol) {
110                 assertEquals(Location.OUTSIDE, south.checkPoint(new S2Point(v)));
111             } else {
112                 assertEquals(Location.BOUNDARY, south.checkPoint(new S2Point(v)));
113             }
114         }
115         assertEquals(1, south.getBoundaryLoops().size());
116 
117         EnclosingBall<Sphere2D, S2Point> southCap = south.getEnclosingCap();
118         assertEquals(0.0, S2Point.MINUS_K.distance(southCap.getCenter()), 1.0e-10);
119         assertEquals(MathUtils.SEMI_PI, southCap.getRadius(), 1.0e-10);
120 
121         EnclosingBall<Sphere2D, S2Point> northCap =
122                 ((SphericalPolygonsSet) new RegionFactory<Sphere2D, S2Point, Circle, SubCircle>().
123                         getComplement(south)).getEnclosingCap();
124         assertEquals(0.0, S2Point.PLUS_K.distance(northCap.getCenter()), 1.0e-10);
125         assertEquals(MathUtils.SEMI_PI, northCap.getRadius(), 1.0e-10);
126 
127         assertEquals(0.0, S2Point.distance(S2Point.MINUS_K, south.getInteriorPoint()), 1.0e-15);
128         assertEquals(Location.INSIDE, south.checkPoint(south.getInteriorPoint()));
129 
130         final RegionFactory<Sphere2D, S2Point, Circle, SubCircle> factory = new RegionFactory<>();
131         final SphericalPolygonsSet reversed = (SphericalPolygonsSet) factory.getComplement(south);
132         assertEquals(0.0, S2Point.distance(S2Point.PLUS_K, reversed.getInteriorPoint()), 1.0e-15);
133         assertEquals(Location.INSIDE, reversed.checkPoint(reversed.getInteriorPoint()));
134 
135     }
136 
137     @Test
138     void testPositiveOctantByIntersection() {
139         double tol = 0.01;
140         double sinTol = FastMath.sin(tol);
141         RegionFactory<Sphere2D, S2Point, Circle, SubCircle> factory = new RegionFactory<>();
142         SphericalPolygonsSet plusX = new SphericalPolygonsSet(Vector3D.PLUS_I, tol);
143         SphericalPolygonsSet plusY = new SphericalPolygonsSet(Vector3D.PLUS_J, tol);
144         SphericalPolygonsSet plusZ = new SphericalPolygonsSet(Vector3D.PLUS_K, tol);
145         SphericalPolygonsSet octant =
146                 (SphericalPolygonsSet) factory.intersection(factory.intersection(plusX, plusY), plusZ);
147         UnitSphereRandomVectorGenerator random =
148                 new UnitSphereRandomVectorGenerator(3, new Well1024a(0x9c9802fde3cbcf25L));
149         for (int i = 0; i < 1000; ++i) {
150             Vector3D v = new Vector3D(random.nextVector());
151             if ((v.getX() > sinTol) && (v.getY() > sinTol) && (v.getZ() > sinTol)) {
152                 assertEquals(Location.INSIDE, octant.checkPoint(new S2Point(v)));
153             } else if ((v.getX() < -sinTol) || (v.getY() < -sinTol) || (v.getZ() < -sinTol)) {
154                 assertEquals(Location.OUTSIDE, octant.checkPoint(new S2Point(v)));
155             } else {
156                 assertEquals(Location.BOUNDARY, octant.checkPoint(new S2Point(v)));
157             }
158         }
159 
160         List<Vertex> loops = octant.getBoundaryLoops();
161         assertEquals(1, loops.size());
162         boolean xPFound = false;
163         boolean yPFound = false;
164         boolean zPFound = false;
165         boolean xVFound = false;
166         boolean yVFound = false;
167         boolean zVFound = false;
168         Vertex first = loops.get(0);
169         int count = 0;
170         for (Vertex v = first; count == 0 || v != first; v = v.getOutgoing().getEnd()) {
171             ++count;
172             Edge e = v.getIncoming();
173             assertSame(v, e.getStart().getOutgoing().getEnd());
174             xPFound = xPFound || e.getCircle().getPole().distance(Vector3D.PLUS_I) < 1.0e-10;
175             yPFound = yPFound || e.getCircle().getPole().distance(Vector3D.PLUS_J) < 1.0e-10;
176             zPFound = zPFound || e.getCircle().getPole().distance(Vector3D.PLUS_K) < 1.0e-10;
177             assertEquals(MathUtils.SEMI_PI, e.getLength(), 1.0e-10);
178             xVFound = xVFound || v.getLocation().getVector().distance(Vector3D.PLUS_I) < 1.0e-10;
179             yVFound = yVFound || v.getLocation().getVector().distance(Vector3D.PLUS_J) < 1.0e-10;
180             zVFound = zVFound || v.getLocation().getVector().distance(Vector3D.PLUS_K) < 1.0e-10;
181         }
182         assertTrue(xPFound);
183         assertTrue(yPFound);
184         assertTrue(zPFound);
185         assertTrue(xVFound);
186         assertTrue(yVFound);
187         assertTrue(zVFound);
188         assertEquals(3, count);
189 
190         assertEquals(0.0, octant.getBarycenter().distance(new S2Point(new Vector3D(1, 1, 1))), 1.0e-10);
191         assertEquals(MathUtils.SEMI_PI, octant.getSize(), 1.0e-10);
192 
193         EnclosingBall<Sphere2D, S2Point> cap = octant.getEnclosingCap();
194         assertEquals(0.0, octant.getBarycenter().distance(cap.getCenter()), 1.0e-10);
195         assertEquals(FastMath.acos(1.0 / FastMath.sqrt(3)), cap.getRadius(), 1.0e-10);
196 
197         EnclosingBall<Sphere2D, S2Point> reversedCap =
198                 ((SphericalPolygonsSet) factory.getComplement(octant)).getEnclosingCap();
199         assertEquals(0, reversedCap.getCenter().distance(new S2Point(new Vector3D(-1, -1, -1))), 1.0e-10);
200         assertEquals(FastMath.PI - FastMath.asin(1.0 / FastMath.sqrt(3)), reversedCap.getRadius(), 1.0e-10);
201 
202         assertEquals(Location.INSIDE, octant.checkPoint(octant.getInteriorPoint()));
203     }
204 
205     @Test
206     void testPositiveOctantByVertices() {
207         double tol = 0.01;
208         double sinTol = FastMath.sin(tol);
209         SphericalPolygonsSet octant = new SphericalPolygonsSet(tol, S2Point.PLUS_I, S2Point.PLUS_J, S2Point.PLUS_K);
210         UnitSphereRandomVectorGenerator random =
211                 new UnitSphereRandomVectorGenerator(3, new Well1024a(0xb8fc5acc91044308L));
212         for (int i = 0; i < 1000; ++i) {
213             Vector3D v = new Vector3D(random.nextVector());
214             if ((v.getX() > sinTol) && (v.getY() > sinTol) && (v.getZ() > sinTol)) {
215                 assertEquals(Location.INSIDE, octant.checkPoint(new S2Point(v)));
216             } else if ((v.getX() < -sinTol) || (v.getY() < -sinTol) || (v.getZ() < -sinTol)) {
217                 assertEquals(Location.OUTSIDE, octant.checkPoint(new S2Point(v)));
218             } else {
219                 assertEquals(Location.BOUNDARY, octant.checkPoint(new S2Point(v)));
220             }
221         }
222 
223         assertEquals(Location.INSIDE, octant.checkPoint(octant.getInteriorPoint()));
224 
225     }
226 
227     @Test
228     void testNonConvex() {
229         double tol = 0.01;
230         double sinTol = FastMath.sin(tol);
231         RegionFactory<Sphere2D, S2Point, Circle, SubCircle> factory = new RegionFactory<>();
232         SphericalPolygonsSet plusX = new SphericalPolygonsSet(Vector3D.PLUS_I, tol);
233         SphericalPolygonsSet plusY = new SphericalPolygonsSet(Vector3D.PLUS_J, tol);
234         SphericalPolygonsSet plusZ = new SphericalPolygonsSet(Vector3D.PLUS_K, tol);
235         SphericalPolygonsSet threeOctants =
236                 (SphericalPolygonsSet) factory.difference(plusZ, factory.intersection(plusX, plusY));
237 
238         UnitSphereRandomVectorGenerator random =
239                 new UnitSphereRandomVectorGenerator(3, new Well1024a(0x9c9802fde3cbcf25L));
240         for (int i = 0; i < 1000; ++i) {
241             Vector3D v = new Vector3D(random.nextVector());
242             if (((v.getX() < -sinTol) || (v.getY() < -sinTol)) && (v.getZ() > sinTol)) {
243                 assertEquals(Location.INSIDE, threeOctants.checkPoint(new S2Point(v)));
244             } else if (((v.getX() > sinTol) && (v.getY() > sinTol)) || (v.getZ() < -sinTol)) {
245                 assertEquals(Location.OUTSIDE, threeOctants.checkPoint(new S2Point(v)));
246             } else {
247                 assertEquals(Location.BOUNDARY, threeOctants.checkPoint(new S2Point(v)));
248             }
249         }
250 
251         List<Vertex> loops = threeOctants.getBoundaryLoops();
252         assertEquals(1, loops.size());
253         boolean xPFound = false;
254         boolean yPFound = false;
255         boolean zPFound = false;
256         boolean xVFound = false;
257         boolean yVFound = false;
258         boolean zVFound = false;
259         Vertex first = loops.get(0);
260         int count = 0;
261         double sumPoleX = 0;
262         double sumPoleY = 0;
263         double sumPoleZ = 0;
264         for (Vertex v = first; count == 0 || v != first; v = v.getOutgoing().getEnd()) {
265             ++count;
266             Edge e = v.getIncoming();
267             assertSame(v, e.getStart().getOutgoing().getEnd());
268             if (e.getCircle().getPole().distance(Vector3D.MINUS_I) < 1.0e-10) {
269                 xPFound = true;
270                 sumPoleX += e.getLength();
271             } else if (e.getCircle().getPole().distance(Vector3D.MINUS_J) < 1.0e-10) {
272                 yPFound = true;
273                 sumPoleY += e.getLength();
274             } else {
275                 assertEquals(0.0, e.getCircle().getPole().distance(Vector3D.PLUS_K), 1.0e-10);
276                 zPFound = true;
277                 sumPoleZ += e.getLength();
278             }
279             xVFound = xVFound || v.getLocation().getVector().distance(Vector3D.PLUS_I) < 1.0e-10;
280             yVFound = yVFound || v.getLocation().getVector().distance(Vector3D.PLUS_J) < 1.0e-10;
281             zVFound = zVFound || v.getLocation().getVector().distance(Vector3D.PLUS_K) < 1.0e-10;
282         }
283         assertTrue(xPFound);
284         assertTrue(yPFound);
285         assertTrue(zPFound);
286         assertTrue(xVFound);
287         assertTrue(yVFound);
288         assertTrue(zVFound);
289         assertEquals(MathUtils.SEMI_PI, sumPoleX, 1.0e-10);
290         assertEquals(MathUtils.SEMI_PI, sumPoleY, 1.0e-10);
291         assertEquals(1.5 * FastMath.PI, sumPoleZ, 1.0e-10);
292 
293         assertEquals(1.5 * FastMath.PI, threeOctants.getSize(), 1.0e-10);
294 
295         assertEquals(Location.INSIDE, threeOctants.checkPoint(threeOctants.getInteriorPoint()));
296     }
297 
298     @Test
299     void testModeratlyComplexShape() {
300         double tol = 0.01;
301         List<SubCircle> boundary = new ArrayList<>();
302         boundary.add(create(Vector3D.MINUS_J, Vector3D.PLUS_I,  Vector3D.PLUS_K,  tol, 0.0, MathUtils.SEMI_PI));
303         boundary.add(create(Vector3D.MINUS_I, Vector3D.PLUS_K,  Vector3D.PLUS_J,  tol, 0.0, MathUtils.SEMI_PI));
304         boundary.add(create(Vector3D.PLUS_K,  Vector3D.PLUS_J,  Vector3D.MINUS_I, tol, 0.0, MathUtils.SEMI_PI));
305         boundary.add(create(Vector3D.MINUS_J, Vector3D.MINUS_I, Vector3D.MINUS_K, tol, 0.0, MathUtils.SEMI_PI));
306         boundary.add(create(Vector3D.MINUS_I, Vector3D.MINUS_K, Vector3D.MINUS_J, tol, 0.0, MathUtils.SEMI_PI));
307         boundary.add(create(Vector3D.PLUS_K,  Vector3D.MINUS_J, Vector3D.PLUS_I,  tol, 0.0, MathUtils.SEMI_PI));
308         SphericalPolygonsSet polygon = new SphericalPolygonsSet(boundary, tol);
309 
310         assertEquals(Location.OUTSIDE, polygon.checkPoint(new S2Point(new Vector3D( 1,  1,  1).normalize())));
311         assertEquals(Location.INSIDE,  polygon.checkPoint(new S2Point(new Vector3D(-1,  1,  1).normalize())));
312         assertEquals(Location.INSIDE,  polygon.checkPoint(new S2Point(new Vector3D(-1, -1,  1).normalize())));
313         assertEquals(Location.INSIDE,  polygon.checkPoint(new S2Point(new Vector3D( 1, -1,  1).normalize())));
314         assertEquals(Location.OUTSIDE, polygon.checkPoint(new S2Point(new Vector3D( 1,  1, -1).normalize())));
315         assertEquals(Location.OUTSIDE, polygon.checkPoint(new S2Point(new Vector3D(-1,  1, -1).normalize())));
316         assertEquals(Location.INSIDE,  polygon.checkPoint(new S2Point(new Vector3D(-1, -1, -1).normalize())));
317         assertEquals(Location.OUTSIDE, polygon.checkPoint(new S2Point(new Vector3D( 1, -1, -1).normalize())));
318 
319         assertEquals(MathUtils.TWO_PI, polygon.getSize(), 1.0e-10);
320         assertEquals(3 * FastMath.PI, polygon.getBoundarySize(), 1.0e-10);
321 
322         List<Vertex> loops = polygon.getBoundaryLoops();
323         assertEquals(1, loops.size());
324         boolean pXFound = false;
325         boolean mXFound = false;
326         boolean pYFound = false;
327         boolean mYFound = false;
328         boolean pZFound = false;
329         boolean mZFound = false;
330         Vertex first = loops.get(0);
331         int count = 0;
332         for (Vertex v = first; count == 0 || v != first; v = v.getOutgoing().getEnd()) {
333             ++count;
334             Edge e = v.getIncoming();
335             assertSame(v, e.getStart().getOutgoing().getEnd());
336             pXFound = pXFound || v.getLocation().getVector().distance(Vector3D.PLUS_I)  < 1.0e-10;
337             mXFound = mXFound || v.getLocation().getVector().distance(Vector3D.MINUS_I) < 1.0e-10;
338             pYFound = pYFound || v.getLocation().getVector().distance(Vector3D.PLUS_J)  < 1.0e-10;
339             mYFound = mYFound || v.getLocation().getVector().distance(Vector3D.MINUS_J) < 1.0e-10;
340             pZFound = pZFound || v.getLocation().getVector().distance(Vector3D.PLUS_K)  < 1.0e-10;
341             mZFound = mZFound || v.getLocation().getVector().distance(Vector3D.MINUS_K) < 1.0e-10;
342             assertEquals(MathUtils.SEMI_PI, e.getLength(), 1.0e-10);
343         }
344         assertTrue(pXFound);
345         assertTrue(mXFound);
346         assertTrue(pYFound);
347         assertTrue(mYFound);
348         assertTrue(pZFound);
349         assertTrue(mZFound);
350         assertEquals(6, count);
351 
352         assertEquals(Location.INSIDE, polygon.checkPoint(polygon.getInteriorPoint()));
353 
354     }
355 
356     @Test
357     void testSeveralParts() {
358         double tol = 0.01;
359         double sinTol = FastMath.sin(tol);
360         List<SubCircle> boundary = new ArrayList<>();
361 
362         // first part: +X, +Y, +Z octant
363         boundary.add(create(Vector3D.PLUS_J,  Vector3D.PLUS_K,  Vector3D.PLUS_I,  tol, 0.0, MathUtils.SEMI_PI));
364         boundary.add(create(Vector3D.PLUS_K,  Vector3D.PLUS_I,  Vector3D.PLUS_J,  tol, 0.0, MathUtils.SEMI_PI));
365         boundary.add(create(Vector3D.PLUS_I,  Vector3D.PLUS_J,  Vector3D.PLUS_K,  tol, 0.0, MathUtils.SEMI_PI));
366 
367         // first part: -X, -Y, -Z octant
368         boundary.add(create(Vector3D.MINUS_J, Vector3D.MINUS_I, Vector3D.MINUS_K, tol, 0.0, MathUtils.SEMI_PI));
369         boundary.add(create(Vector3D.MINUS_I, Vector3D.MINUS_K, Vector3D.MINUS_J, tol, 0.0, MathUtils.SEMI_PI));
370         boundary.add(create(Vector3D.MINUS_K, Vector3D.MINUS_J, Vector3D.MINUS_I,  tol, 0.0, MathUtils.SEMI_PI));
371 
372         SphericalPolygonsSet polygon = new SphericalPolygonsSet(boundary, tol);
373 
374         UnitSphereRandomVectorGenerator random =
375                 new UnitSphereRandomVectorGenerator(3, new Well1024a(0xcc5ce49949e0d3ecL));
376         for (int i = 0; i < 1000; ++i) {
377             Vector3D v = new Vector3D(random.nextVector());
378             if ((v.getX() < -sinTol) && (v.getY() < -sinTol) && (v.getZ() < -sinTol)) {
379                 assertEquals(Location.INSIDE, polygon.checkPoint(new S2Point(v)));
380             } else if ((v.getX() < sinTol) && (v.getY() < sinTol) && (v.getZ() < sinTol)) {
381                 assertEquals(Location.BOUNDARY, polygon.checkPoint(new S2Point(v)));
382             } else if ((v.getX() > sinTol) && (v.getY() > sinTol) && (v.getZ() > sinTol)) {
383                 assertEquals(Location.INSIDE, polygon.checkPoint(new S2Point(v)));
384             } else if ((v.getX() > -sinTol) && (v.getY() > -sinTol) && (v.getZ() > -sinTol)) {
385                 assertEquals(Location.BOUNDARY, polygon.checkPoint(new S2Point(v)));
386             } else {
387                 assertEquals(Location.OUTSIDE, polygon.checkPoint(new S2Point(v)));
388             }
389         }
390 
391         assertEquals(FastMath.PI, polygon.getSize(), 1.0e-10);
392         assertEquals(3 * FastMath.PI, polygon.getBoundarySize(), 1.0e-10);
393 
394         // there should be two separate boundary loops
395         assertEquals(2, polygon.getBoundaryLoops().size());
396 
397         assertEquals(Location.INSIDE, polygon.checkPoint(polygon.getInteriorPoint()));
398 
399     }
400 
401     @Test
402     void testPartWithHole() {
403         double tol = 0.01;
404         double alpha = 0.7;
405         S2Point center = new S2Point(new Vector3D(1, 1, 1));
406         SphericalPolygonsSet hexa = new SphericalPolygonsSet(center.getVector(), Vector3D.PLUS_K, alpha, 6, tol);
407         SphericalPolygonsSet hole  = new SphericalPolygonsSet(tol,
408                                                               new S2Point(FastMath.PI / 6, FastMath.PI / 3),
409                                                               new S2Point(FastMath.PI / 3, FastMath.PI / 3),
410                                                               new S2Point(FastMath.PI / 4, FastMath.PI / 6));
411         SphericalPolygonsSet hexaWithHole =
412                 (SphericalPolygonsSet) new RegionFactory<Sphere2D, S2Point, Circle, SubCircle>().
413                         difference(hexa, hole);
414 
415         for (double phi = center.getPhi() - alpha + 0.1; phi < center.getPhi() + alpha - 0.1; phi += 0.07) {
416             Location l = hexaWithHole.checkPoint(new S2Point(FastMath.PI / 4, phi));
417             if (phi < FastMath.PI / 6 || phi > FastMath.PI / 3) {
418                 assertEquals(Location.INSIDE,  l);
419             } else {
420                 assertEquals(Location.OUTSIDE, l);
421             }
422         }
423 
424         // there should be two separate boundary loops
425         assertEquals(2, hexaWithHole.getBoundaryLoops().size());
426 
427         assertEquals(hexa.getBoundarySize() + hole.getBoundarySize(), hexaWithHole.getBoundarySize(), 1.0e-10);
428         assertEquals(hexa.getSize() - hole.getSize(), hexaWithHole.getSize(), 1.0e-10);
429 
430         assertEquals(Location.INSIDE, hexaWithHole.checkPoint(hexaWithHole.getInteriorPoint()));
431 
432     }
433 
434     @Test
435     void testConcentricSubParts() {
436         double tol = 0.001;
437         Vector3D center = new Vector3D(1, 1, 1);
438         SphericalPolygonsSet hexaOut   = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.9,  6, tol);
439         SphericalPolygonsSet hexaIn    = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.8,  6, tol);
440         SphericalPolygonsSet pentaOut  = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.7,  5, tol);
441         SphericalPolygonsSet pentaIn   = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.6,  5, tol);
442         SphericalPolygonsSet quadriOut = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.5,  4, tol);
443         SphericalPolygonsSet quadriIn  = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.4,  4, tol);
444         SphericalPolygonsSet triOut    = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.25, 3, tol);
445         SphericalPolygonsSet triIn     = new SphericalPolygonsSet(center, Vector3D.PLUS_K, 0.15, 3, tol);
446 
447         RegionFactory<Sphere2D, S2Point, Circle, SubCircle> factory = new RegionFactory<>();
448         SphericalPolygonsSet hexa   = (SphericalPolygonsSet) factory.difference(hexaOut,   hexaIn);
449         SphericalPolygonsSet penta  = (SphericalPolygonsSet) factory.difference(pentaOut,  pentaIn);
450         SphericalPolygonsSet quadri = (SphericalPolygonsSet) factory.difference(quadriOut, quadriIn);
451         SphericalPolygonsSet tri    = (SphericalPolygonsSet) factory.difference(triOut,    triIn);
452         SphericalPolygonsSet concentric =
453                 (SphericalPolygonsSet) factory.union(factory.union(hexa, penta), factory.union(quadri, tri));
454 
455         // there should be two separate boundary loops
456         assertEquals(8, concentric.getBoundaryLoops().size());
457 
458         assertEquals(hexaOut.getBoundarySize()   + hexaIn.getBoundarySize()   +
459                             pentaOut.getBoundarySize()  + pentaIn.getBoundarySize()  +
460                             quadriOut.getBoundarySize() + quadriIn.getBoundarySize() +
461                             triOut.getBoundarySize()    + triIn.getBoundarySize(),
462                             concentric.getBoundarySize(), 1.0e-10);
463         assertEquals(hexaOut.getSize()   - hexaIn.getSize()   +
464                             pentaOut.getSize()  - pentaIn.getSize()  +
465                             quadriOut.getSize() - quadriIn.getSize() +
466                             triOut.getSize()    - triIn.getSize(),
467                             concentric.getSize(), 1.0e-10);
468 
469         // we expect lots of sign changes as we traverse all concentric rings
470         double phi = new S2Point(center).getPhi();
471         assertEquals(+0.207, concentric.projectToBoundary(new S2Point(-0.60,  phi)).getOffset(), 0.01);
472         assertEquals(-0.048, concentric.projectToBoundary(new S2Point(-0.21,  phi)).getOffset(), 0.01);
473         assertEquals(+0.027, concentric.projectToBoundary(new S2Point(-0.10,  phi)).getOffset(), 0.01);
474         assertEquals(-0.041, concentric.projectToBoundary(new S2Point( 0.01,  phi)).getOffset(), 0.01);
475         assertEquals(+0.049, concentric.projectToBoundary(new S2Point( 0.16,  phi)).getOffset(), 0.01);
476         assertEquals(-0.038, concentric.projectToBoundary(new S2Point( 0.29,  phi)).getOffset(), 0.01);
477         assertEquals(+0.097, concentric.projectToBoundary(new S2Point( 0.48,  phi)).getOffset(), 0.01);
478         assertEquals(-0.022, concentric.projectToBoundary(new S2Point( 0.64,  phi)).getOffset(), 0.01);
479         assertEquals(+0.072, concentric.projectToBoundary(new S2Point( 0.79,  phi)).getOffset(), 0.01);
480         assertEquals(-0.022, concentric.projectToBoundary(new S2Point( 0.93,  phi)).getOffset(), 0.01);
481         assertEquals(+0.091, concentric.projectToBoundary(new S2Point( 1.08,  phi)).getOffset(), 0.01);
482         assertEquals(-0.037, concentric.projectToBoundary(new S2Point( 1.28,  phi)).getOffset(), 0.01);
483         assertEquals(+0.051, concentric.projectToBoundary(new S2Point( 1.40,  phi)).getOffset(), 0.01);
484         assertEquals(-0.041, concentric.projectToBoundary(new S2Point( 1.55,  phi)).getOffset(), 0.01);
485         assertEquals(+0.027, concentric.projectToBoundary(new S2Point( 1.67,  phi)).getOffset(), 0.01);
486         assertEquals(-0.044, concentric.projectToBoundary(new S2Point( 1.79,  phi)).getOffset(), 0.01);
487         assertEquals(+0.201, concentric.projectToBoundary(new S2Point( 2.16,  phi)).getOffset(), 0.01);
488 
489         assertEquals(Location.INSIDE, concentric.checkPoint(concentric.getInteriorPoint()));
490 
491     }
492 
493     @Test
494     void testGeographicalMap() {
495 
496         SphericalPolygonsSet continental = buildSimpleZone(new double[][] {
497           { 51.14850,  2.51357 }, { 50.94660,  1.63900 }, { 50.12717,  1.33876 }, { 49.34737, -0.98946 },
498           { 49.77634, -1.93349 }, { 48.64442, -1.61651 }, { 48.90169, -3.29581 }, { 48.68416, -4.59234 },
499           { 47.95495, -4.49155 }, { 47.57032, -2.96327 }, { 46.01491, -1.19379 }, { 44.02261, -1.38422 },
500           { 43.42280, -1.90135 }, { 43.03401, -1.50277 }, { 42.34338,  1.82679 }, { 42.47301,  2.98599 },
501           { 43.07520,  3.10041 }, { 43.39965,  4.55696 }, { 43.12889,  6.52924 }, { 43.69384,  7.43518 },
502           { 44.12790,  7.54959 }, { 45.02851,  6.74995 }, { 45.33309,  7.09665 }, { 46.42967,  6.50009 },
503           { 46.27298,  6.02260 }, { 46.72577,  6.03738 }, { 47.62058,  7.46675 }, { 49.01778,  8.09927 },
504           { 49.20195,  6.65822 }, { 49.44266,  5.89775 }, { 49.98537,  4.79922 }
505         });
506         SphericalPolygonsSet corsica = buildSimpleZone(new double[][] {
507           { 42.15249,  9.56001 }, { 43.00998,  9.39000 }, { 42.62812,  8.74600 }, { 42.25651,  8.54421 },
508           { 41.58361,  8.77572 }, { 41.38000,  9.22975 }
509         });
510         RegionFactory<Sphere2D, S2Point, Circle, SubCircle> factory = new RegionFactory<>();
511         SphericalPolygonsSet zone = (SphericalPolygonsSet) factory.union(continental, corsica);
512         EnclosingBall<Sphere2D, S2Point> enclosing = zone.getEnclosingCap();
513         Vector3D enclosingCenter = enclosing.getCenter().getVector();
514 
515         double step = FastMath.toRadians(0.1);
516         for (Vertex loopStart : zone.getBoundaryLoops()) {
517             int count = 0;
518             for (Vertex v = loopStart; count == 0 || v != loopStart; v = v.getOutgoing().getEnd()) {
519                 ++count;
520                 for (int i = 0; i < FastMath.ceil(v.getOutgoing().getLength() / step); ++i) {
521                     Vector3D p = v.getOutgoing().getPointAt(i * step);
522                     assertTrue(Vector3D.angle(p, enclosingCenter) <= enclosing.getRadius());
523                 }
524             }
525         }
526 
527         S2Point supportPointA = s2Point(48.68416, -4.59234);
528         S2Point supportPointB = s2Point(41.38000,  9.22975);
529         assertEquals(enclosing.getRadius(), supportPointA.distance(enclosing.getCenter()), 1.0e-10);
530         assertEquals(enclosing.getRadius(), supportPointB.distance(enclosing.getCenter()), 1.0e-10);
531         assertEquals(0.5 * supportPointA.distance(supportPointB), enclosing.getRadius(), 1.0e-10);
532         assertEquals(2, enclosing.getSupportSize());
533 
534         EnclosingBall<Sphere2D, S2Point> continentalInscribed =
535                 ((SphericalPolygonsSet) factory.getComplement(continental)).getEnclosingCap();
536         Vector3D continentalCenter = continentalInscribed.getCenter().getVector();
537         assertEquals(2.2, FastMath.toDegrees(FastMath.PI - continentalInscribed.getRadius()), 0.1);
538         for (Vertex loopStart : continental.getBoundaryLoops()) {
539             int count = 0;
540             for (Vertex v = loopStart; count == 0 || v != loopStart; v = v.getOutgoing().getEnd()) {
541                 ++count;
542                 for (int i = 0; i < FastMath.ceil(v.getOutgoing().getLength() / step); ++i) {
543                     Vector3D p = v.getOutgoing().getPointAt(i * step);
544                     assertTrue(Vector3D.angle(p, continentalCenter) <= continentalInscribed.getRadius());
545                 }
546             }
547         }
548 
549         EnclosingBall<Sphere2D, S2Point> corsicaInscribed =
550                 ((SphericalPolygonsSet) factory.getComplement(corsica)).getEnclosingCap();
551         Vector3D corsicaCenter = corsicaInscribed.getCenter().getVector();
552         assertEquals(0.34, FastMath.toDegrees(FastMath.PI - corsicaInscribed.getRadius()), 0.01);
553         for (Vertex loopStart : corsica.getBoundaryLoops()) {
554             int count = 0;
555             for (Vertex v = loopStart; count == 0 || v != loopStart; v = v.getOutgoing().getEnd()) {
556                 ++count;
557                 for (int i = 0; i < FastMath.ceil(v.getOutgoing().getLength() / step); ++i) {
558                     Vector3D p = v.getOutgoing().getPointAt(i * step);
559                     assertTrue(Vector3D.angle(p, corsicaCenter) <= corsicaInscribed.getRadius());
560                 }
561             }
562         }
563 
564         assertEquals(Location.INSIDE, zone.checkPoint(zone.getInteriorPoint()));
565 
566         // TODO: this triggers an exception as we cannot find an interior point
567 //        checkInterior(FastMath.toRadians(-2), FastMath.toRadians(5), FastMath.toRadians(0.4),
568 //                      FastMath.toRadians(90 - 55), FastMath.toRadians(90 - 40), FastMath.toRadians(0.4),
569 //                      zone, 1.0e-8);
570 
571     }
572 
573     @Test
574     void testZigZagBoundary() {
575         SphericalPolygonsSet zone = new SphericalPolygonsSet(1.0e-6,
576                                                              new S2Point(-0.12630940610562444, 0.8998192093789258),
577                                                              new S2Point(-0.12731320182988207, 0.8963735568774486),
578                                                              new S2Point(-0.1351107624622557,  0.8978258663483273),
579                                                              new S2Point(-0.13545331405131725, 0.8966781238246179),
580                                                              new S2Point(-0.14324883017454967, 0.8981309629283796),
581                                                              new S2Point(-0.14359875625524995, 0.896983965573036),
582                                                              new S2Point(-0.14749650541159384, 0.8977109994666864),
583                                                              new S2Point(-0.14785037758231825, 0.8965644005442432),
584                                                              new S2Point(-0.15369807257448784, 0.8976550608135502),
585                                                              new S2Point(-0.1526225554339386,  0.9010934265410458),
586                                                              new S2Point(-0.14679028466684121, 0.9000043396997698),
587                                                              new S2Point(-0.14643807494172612, 0.9011511073761742),
588                                                              new S2Point(-0.1386609051963748,  0.8996991539048602),
589                                                              new S2Point(-0.13831601655974668, 0.9008466623902937),
590                                                              new S2Point(-0.1305365419828323,  0.8993961857946309),
591                                                              new S2Point(-0.1301989630405964,  0.9005444294061787));
592         assertEquals(Region.Location.INSIDE, zone.checkPoint(new S2Point(-0.145, 0.898)));
593         assertEquals(6.463e-5, zone.getSize(),         1.0e-7);
594         assertEquals(5.487e-2, zone.getBoundarySize(), 1.0e-4);
595     }
596 
597     @Test
598     void testGitHubIssue41() {
599         RegionFactory<Sphere2D, S2Point, Circle, SubCircle> regionFactory = new RegionFactory<>();
600         S2Point[] s2pA = new S2Point[]{
601                 new S2Point(new Vector3D(0.2122954606, -0.629606302,  0.7473463333)),
602                 new S2Point(new Vector3D(0.2120220248, -0.6296445493, 0.747391733)),
603                 new S2Point(new Vector3D(0.2119838016, -0.6298173178, 0.7472569934)),
604                 new S2Point(new Vector3D(0.2122571927, -0.6297790738, 0.7472116182))};
605 
606         S2Point[] s2pB = new S2Point[]{
607                 new S2Point(new Vector3D(0.2120291561, -0.629952069,  0.7471305292)),
608                 new S2Point(new Vector3D(0.2123026002, -0.6299138005, 0.7470851423)),
609                 new S2Point(new Vector3D(0.2123408927, -0.6297410403, 0.7472198923)),
610                 new S2Point(new Vector3D(0.2120674039, -0.6297793122, 0.7472653037))};
611 
612         final double tol = 0.0001;
613         final SphericalPolygonsSet spsA = new SphericalPolygonsSet(tol, s2pA);
614         final SphericalPolygonsSet spsB = new SphericalPolygonsSet(tol, s2pB);
615         assertEquals(0.61254e-7, spsA.getSize(),          1.0e-12);
616         assertEquals(1.00437e-3, spsA.getBoundarySize(),  1.0e-08);
617         assertEquals(0.61269e-7, spsB.getSize(),          1.0e-12);
618         assertEquals(1.00452e-3, spsB.getBoundarySize(),  1.0e-08);
619         SphericalPolygonsSet union = (SphericalPolygonsSet) regionFactory.union(spsA, spsB);
620 
621         // as the tolerance is very large with respect to polygons,
622         // the union is not computed properly, which is EXPECTED
623         // so the thresholds for the tests are large.
624         // the reference values have been computed with a much lower tolerance
625         assertEquals(1.15628e-7, union.getSize(),         4.0e-9);
626         assertEquals(1.53824e-3, union.getBoundarySize(), 3.0e-4);
627 
628     }
629 
630     @Test
631     void testGitHubIssue42A() {
632         // if building it was allowed (i.e. if the check for tolerance was removed)
633         // the BSP tree would be wrong, it would include a large extra chunk that contains
634         // a point that should really be outside
635         try {
636             doTestGitHubIssue42(1.0e-100);
637         } catch (MathIllegalArgumentException miae) {
638             assertEquals(LocalizedGeometryFormats.TOO_SMALL_TOLERANCE, miae.getSpecifier());
639             assertEquals(1.0e-100, (Double) miae.getParts()[0], 1.0e-110);
640             assertEquals("Sphere2D.SMALLEST_TOLERANCE", miae.getParts()[1]);
641             assertEquals(Sphere2D.SMALLEST_TOLERANCE, (Double) miae.getParts()[2], 1.0e-20);
642         }
643     }
644 
645     @Test
646     void testGitHubIssue42B() {
647         // the BSP tree is right, but size cannot be computed
648         try {
649             // success of this call is dependent on numerical noise.
650             // If it fails it should fail predictably.
651             doTestGitHubIssue42(9.0e-16);
652         } catch (MathIllegalStateException e) {
653             assertEquals(
654                 LocalizedGeometryFormats.OUTLINE_BOUNDARY_LOOP_OPEN, e.getSpecifier());
655         }
656         // Computations in Edge.split use angles up to 4 PI
657         double tol = FastMath.ulp(4 * FastMath.PI);
658         // works when tol >= ulp(largest angle used)
659         doTestGitHubIssue42(tol);
660     }
661 
662     private void doTestGitHubIssue42(double tolerance) throws MathIllegalArgumentException {
663         S2Point[] s2pA = new S2Point[]{
664             new S2Point(new Vector3D(0.1504230736114679,  -0.6603084987333554, 0.7357754993377947)),
665             new S2Point(new Vector3D(0.15011191112224423, -0.6603400871954631, 0.7358106980616113)),
666             new S2Point(new Vector3D(0.15008035620222715, -0.6605195692153062, 0.7356560238085725)),
667             new S2Point(new Vector3D(0.1503914563063968,  -0.6604879854490165, 0.7356208472763267))
668         };
669         S2Point outsidePoint = new S2Point(new Vector3D( 2, s2pA[0].getVector(),
670                                                         -1, s2pA[1].getVector(),
671                                                         -1, s2pA[2].getVector(),
672                                                          2, s2pA[3].getVector()).normalize());
673         // test all permutations because order matters when tolerance is small
674         for (int i = 0; i < s2pA.length; i++) {
675             S2Point[] points = new S2Point[s2pA.length];
676             for (int j = 0; j < s2pA.length; j++) {
677                 points[j] = s2pA[(i + j) % s2pA.length];
678             }
679             final SphericalPolygonsSet spsA = new SphericalPolygonsSet(tolerance, points);
680             assertEquals(Location.OUTSIDE, spsA.checkPoint(outsidePoint));
681             assertEquals(7.4547e-8, spsA.getSize(), 1.0e-12);
682         }
683     }
684 
685     @Test
686     void testZigZagBoundaryOversampledIssue46() {
687         final double tol = 1.0e-4;
688         // sample region, non-convex, not too big, not too small
689         final S2Point[] vertices = {
690                 new S2Point(-0.12630940610562444e1, (0.8998192093789258 - 0.89) * 100),
691                 new S2Point(-0.12731320182988207e1, (0.8963735568774486 - 0.89) * 100),
692                 new S2Point(-0.1351107624622557e1, (0.8978258663483273 - 0.89) * 100),
693                 new S2Point(-0.13545331405131725e1, (0.8966781238246179 - 0.89) * 100),
694                 new S2Point(-0.14324883017454967e1, (0.8981309629283796 - 0.89) * 100),
695                 new S2Point(-0.14359875625524995e1, (0.896983965573036 - 0.89) * 100),
696                 new S2Point(-0.14749650541159384e1, (0.8977109994666864 - 0.89) * 100),
697                 new S2Point(-0.14785037758231825e1, (0.8965644005442432 - 0.89) * 100),
698                 new S2Point(-0.15369807257448784e1, (0.8976550608135502 - 0.89) * 100),
699                 new S2Point(-0.1526225554339386e1, (0.9010934265410458 - 0.89) * 100),
700                 new S2Point(-0.14679028466684121e1, (0.9000043396997698 - 0.89) * 100),
701                 new S2Point(-0.14643807494172612e1, (0.9011511073761742 - 0.89) * 100),
702                 new S2Point(-0.1386609051963748e1, (0.8996991539048602 - 0.89) * 100),
703                 new S2Point(-0.13831601655974668e1, (0.9008466623902937 - 0.89) * 100),
704                 new S2Point(-0.1305365419828323e1, (0.8993961857946309 - 0.89) * 100),
705                 new S2Point(-0.1301989630405964e1, (0.9005444294061787 - 0.89) * 100)};
706         SphericalPolygonsSet zone = new SphericalPolygonsSet(tol, vertices);
707         // sample high resolution boundary
708         List<S2Point> points = new ArrayList<>();
709         final Vertex start = zone.getBoundaryLoops().get(0);
710         Vertex v = start;
711         double step = tol / 10;
712         do {
713             Edge outgoing = v.getOutgoing();
714             final double length = outgoing.getLength();
715             int n = (int) (length / step);
716             for (int i = 0; i < n; i++) {
717                 points.add(new S2Point(outgoing.getPointAt(i * step)));
718             }
719             v = outgoing.getEnd();
720         } while (v != start);
721         // create zone from high resolution boundary
722         zone = new SphericalPolygonsSet(tol, points.toArray(new S2Point[0]));
723         EnclosingBall<Sphere2D, S2Point> cap = zone.getEnclosingCap();
724         // check cap size is reasonable. The region is ~0.5 accross, could be < 0.25
725         assertTrue(cap.getRadius() < 0.5);
726         assertEquals(Location.INSIDE, zone.checkPoint(zone.getBarycenter()));
727         assertEquals(Location.INSIDE, zone.checkPoint(cap.getCenter()));
728         // extra tolerance at corners due to SPS tolerance being a hyperplaneThickness
729         // a factor of 3.1 corresponds to a edge intersection angle of ~19 degrees
730         final double cornerTol = 3.1 * tol;
731         for (S2Point vertex : vertices) {
732             // check original points are on the boundary
733             assertEquals(Location.BOUNDARY, zone.checkPoint(vertex), "" + vertex);
734             double offset = FastMath.abs(zone.projectToBoundary(vertex).getOffset());
735             assertEquals(0, offset, cornerTol, vertex + " offset: " + offset);
736             // check original points are within the cap
737             assertTrue(
738                     cap.contains(vertex, tol),
739                     "vertex: " + vertex + " distance: " + (vertex.distance(cap.getCenter()) - cap.getRadius()));
740         }
741     }
742 
743     @Test
744     void testPositiveOctantByVerticesDetailIssue46() {
745         double tol = 0.01;
746         double sinTol = FastMath.sin(tol);
747         Circle x = new Circle(Vector3D.PLUS_I, tol);
748         Circle z = new Circle(Vector3D.PLUS_K, tol);
749         double length = FastMath.PI / 2;
750         double step = tol / 10;
751         // sample high resolution boundary
752         int n = (int) (length / step);
753         List<S2Point> points = new ArrayList<>();
754         for (int i = 0; i < n; i++) {
755             double t = i * step;
756             points.add(new S2Point(z.getPointAt(z.getPhase(Vector3D.PLUS_I) + t)));
757         }
758         for (int i = 0; i < n; i++) {
759             double t = i * step;
760             points.add(new S2Point(x.getPointAt(x.getPhase(Vector3D.PLUS_J) + t)));
761         }
762         points.add(S2Point.PLUS_K);
763 
764         SphericalPolygonsSet octant = new SphericalPolygonsSet(tol, points.toArray(new S2Point[0]));
765         UnitSphereRandomVectorGenerator random =
766                 new UnitSphereRandomVectorGenerator(3, new Well1024a(0xb8fc5acc91044308L));
767         /* Where exactly the boundaries fall depends on which points are kept from
768          * decimation, which can vary by up to tol. So a point up to 2*tol away from a
769          * input point may be on the boundary. All input points are guaranteed to be on
770          * the boundary, just not the center of the boundary.
771          */
772         for (int i = 0; i < 1000; ++i) {
773             Vector3D v = new Vector3D(random.nextVector());
774             final Location actual = octant.checkPoint(new S2Point(v));
775             if ((v.getX() > sinTol) && (v.getY() > sinTol) && (v.getZ() > sinTol)) {
776                 if ((v.getX() > 2*sinTol) && (v.getY() > 2*sinTol) && (v.getZ() > 2*sinTol)) {
777                     // certainly inside
778                     assertEquals(Location.INSIDE, actual, "" + v);
779                 } else {
780                     // may be inside or boundary
781                     assertNotEquals(Location.OUTSIDE, actual, "" + v);
782                 }
783             } else if ((v.getX() < 0) || (v.getY() < 0) || (v.getZ() < 0)) {
784                 if ((v.getX() < -sinTol) || (v.getY() < -sinTol) || (v.getZ() < -sinTol)) {
785                     // certainly outside
786                     assertEquals(Location.OUTSIDE, actual);
787                 } else {
788                     // may be outside or boundary
789                     assertNotEquals(Location.INSIDE, actual);
790                 }
791             } else {
792                 // certainly boundary
793                 assertEquals(Location.BOUNDARY, actual);
794             }
795         }
796         // all input points are on the boundary
797         for (S2Point point : points) {
798             assertEquals(Location.BOUNDARY, octant.checkPoint(point), "" + point);
799         }
800     }
801 
802     /** Check that constructing a region from 0, 1, 2, 3 points along a circle works. */
803     @Test
804     void testConstructingFromFewPointsIssue46() {
805         double tol  = 1e-9;
806         List<S2Point> points = new ArrayList<>();
807         Circle circle = new Circle(Vector3D.PLUS_K, tol);
808 
809         // no points
810         SphericalPolygonsSet sps = new SphericalPolygonsSet(tol, new S2Point[0]);
811         assertEquals(4*FastMath.PI, sps.getSize(), tol);
812 
813         // one point, does not define a valid boundary
814         points.add(new S2Point(circle.getPointAt(0)));
815         try {
816             new SphericalPolygonsSet(tol, points.toArray(new S2Point[0]));
817             fail("expcected exception");
818         } catch (MathRuntimeException e) {
819             // expected
820         }
821 
822         // two points, defines hemisphere but not orientation
823         points.add(new S2Point(circle.getPointAt(FastMath.PI / 2)));
824         sps = new SphericalPolygonsSet(tol, points.toArray(new S2Point[0]));
825         assertEquals(2*FastMath.PI, sps.getSize(), tol);
826 
827         // three points
828         points.add(0, new S2Point(circle.getPointAt(FastMath.PI)));
829         sps = new SphericalPolygonsSet(tol, points.toArray(new S2Point[0]));
830         assertEquals(2*FastMath.PI, sps.getSize(), tol);
831         assertEquals(0, sps.getBarycenter().distance(new S2Point(Vector3D.PLUS_K)), tol);
832 
833         // four points
834         points.add(1, new S2Point(circle.getPointAt(3 * FastMath.PI / 2)));
835         sps = new SphericalPolygonsSet(tol, points.toArray(new S2Point[0]));
836         assertEquals(2 * FastMath.PI, sps.getSize(), tol);
837         assertEquals(0, sps.getBarycenter().distance(new S2Point(Vector3D.PLUS_K)), tol);
838 
839         // many points in semi-circle
840         sps = new SphericalPolygonsSet(tol,
841                 new S2Point(circle.getPointAt(0)),
842                 new S2Point(circle.getPointAt(-0.3)),
843                 new S2Point(circle.getPointAt(-0.2)),
844                 new S2Point(circle.getPointAt(-0.1)));
845         assertEquals(2 * FastMath.PI, sps.getSize(), tol);
846         assertEquals(0, sps.getBarycenter().distance(new S2Point(Vector3D.PLUS_K)), tol);
847     }
848 
849     @Test
850     void testDefensiveProgrammingCheck() {
851         // this tests defensive programming code that seems almost unreachable otherwise
852         try {
853             Method searchHelper = SphericalPolygonsSet.class.getDeclaredMethod("searchHelper",
854                                                                                IntPredicate.class,
855                                                                                Integer.TYPE, Integer.TYPE);
856             searchHelper.setAccessible(true);
857             searchHelper.invoke(null, (IntPredicate) (n -> true), 1, 0);
858             fail("an exception should have been thrown");
859         } catch (InvocationTargetException ite) {
860             MathIllegalArgumentException miae = (MathIllegalArgumentException) ite.getCause();
861             assertEquals(LocalizedCoreFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT, miae.getSpecifier());
862         } catch (NoSuchMethodException | IllegalAccessException | IllegalArgumentException e) {
863             fail(e.getLocalizedMessage());
864         }
865     }
866 
867     /**
868      * Tests the Hipparchus {@link RegionFactory#intersection(Region, Region)}
869      * method.
870      */
871     @Test
872     void TestIntersectionOrder() {
873 
874         final S2Point[] vertices1 = {
875             new S2Point(0.0193428339344826, 1.5537209444301618),
876             new S2Point(0.0178197572212936, 1.553415699912148),
877             new S2Point(0.01628496406053076, 1.5531081515279537),
878             new S2Point(0.016284670226196844, 1.5531096373947835),
879             new S2Point(0.019342540199680208, 1.5537224293848613)
880         };
881 
882         final S2Point[] vertices2 = {
883             new S2Point(0.016, 1.555),
884             new S2Point(0.017453292519943295, 1.555),
885             new S2Point(0.017453292519943295, 1.5533430342749532),
886             new S2Point(0.016, 1.5533430342749532)
887         };
888 
889         final RegionFactory<Sphere2D, S2Point, Circle, SubCircle> regionFactory = new RegionFactory<>();
890 
891         // thickness is small enough for proper computation of very small intersection
892         double thickness1 = 4.96740426e-11;
893         final SphericalPolygonsSet sps1 = new SphericalPolygonsSet(thickness1, vertices1);
894         final SphericalPolygonsSet sps2 = new SphericalPolygonsSet(thickness1, vertices2);
895         assertEquals(1.4886e-12, regionFactory.intersection(sps1, sps2).getSize(), 1.0e-15);
896         assertEquals(1.4881e-12, regionFactory.intersection(sps2, sps1).getSize(), 1.0e-15);
897 
898         // thickness is too large, very small intersection is not computed properly in one case
899         double thickness2 = 4.96740427e-11;
900         final SphericalPolygonsSet sps3 = new SphericalPolygonsSet(thickness2, vertices1);
901         final SphericalPolygonsSet sps4 = new SphericalPolygonsSet(thickness2, vertices2);
902         assertEquals(1.4886e-12, regionFactory.intersection(sps3, sps4).getSize(), 1.0e-15);
903         assertEquals(0.0,        regionFactory.intersection(sps4, sps3).getSize(), 1.0e-15);
904 
905     }
906 
907     @Test
908     public void testIssueOrekit1388A() {
909         doTestIssueOrekit1388(true);
910     }
911 
912     @Test
913     public void testIssueOrekit1388B() {
914         doTestIssueOrekit1388(false);
915     }
916 
917     private void doTestIssueOrekit1388(final boolean order) {
918         final double[][] coordinates1 = new double[][] {
919             { 18.52684751402596,  -76.97880893719434 },
920             { 18.451108862175584, -76.99484778988442 },
921             { 18.375369256045143, -77.01087679277504 },
922             { 18.299628701801268, -77.02689599675749 },
923             { 18.223887203723567, -77.04290545732495 },
924             { 18.148144771769385, -77.05890521550272 },
925             { 18.072401410187016, -77.0748953260324  },
926             { 17.99665712514784,  -77.09087584010511 },
927             { 17.92091192468999,  -77.10684680260262 },
928             { 17.84516581113023,  -77.12280827309398 },
929             { 17.769418792522664, -77.13876029654094 },
930             { 17.747659863099422, -77.14334087084347 },
931             { 17.67571798336192,  -77.15846791369165 },
932             { 17.624293265977183, -77.16928381433733 },
933             { 17.5485398681768,   -77.18520934447962 },
934             { 17.526779103104783, -77.1897823275402  },
935             { 17.49650619905315,  -77.0342031192472  },
936             { 17.588661518962343, -77.01473903648854 },
937             { 17.728574326965138, -76.98517769352242 },
938             { 17.80416324015021,  -76.96919708557023 },
939             { 17.87969526622326,  -76.95321858415689 },
940             { 17.955280973332677, -76.93721874766547 },
941             { 18.030855567607098, -76.92121123297645 },
942             { 18.106414929680927, -76.90519686376611 },
943             { 18.182031502555215, -76.88916022728444 },
944             { 18.257597934434987, -76.87312403715188 },
945             { 18.3331742522667,   -76.85707550881591 },
946             { 18.408750874895002, -76.84101662269072 },
947             { 18.57249082100609,  -76.80620195239239 },
948             { 18.602585205425896, -76.96276018365842 }
949         };
950 
951         final double[][] coordinates2 = new double[][] {
952             { 18.338614038907608, -78.37885677406668 },
953             { 18.195574802144037, -78.24425107003432 },
954             { 18.20775293886321,  -78.0711865934217  },
955             { 18.07679345301507,  -77.95901517339438 },
956             { 18.006705181057598, -77.85325354879791 },
957             { 17.857293838883137, -77.73787723105598 },
958             { 17.854243316622103, -77.57442744758828 },
959             { 17.875595873376014, -77.38213358468467 },
960             { 17.72607423578937,  -77.23470828979222 },
961             { 17.71386286451302,  -77.12253686976543 },
962             { 17.790170276013725, -77.14817605148616 },
963             { 17.869495404611797, -77.14497115377101 },
964             { 17.854243309397717, -76.9302429967729  },
965             { 17.954882874700132, -76.84371075846688 },
966             { 17.94268718313505,  -76.6898756681441  },
967             { 17.869495397388064, -76.54886016868198 },
968             { 17.863394719203555, -76.35015651034861 },
969             { 17.93049065091843,  -76.23478019260665 },
970             { 18.155989976776553, -76.32451732862788 },
971             { 18.22601854027039,  -76.63218750927341 },
972             { 18.33861403170316,  -76.85653034932697 },
973             { 18.405527980074993, -76.97831646249921 },
974             { 18.4541763474828,   -77.28598664314421 },
975             { 18.496732365966466, -77.705828243816   },
976             { 18.451136227912485, -78.00708862903122 },
977             { 18.405527980074993, -78.25707065080552 }
978         };
979 
980         final double[][] expectedIn = new double[][] {
981                 { 18.408, -77.003 },
982                 { 18.338, -76.857 },
983                 { 17.869, -77.117 },
984                 { 17.857, -76.959 },
985                 { 17.761, -77.139 },
986                 { 17.715, -77.125 },
987                 { 17.935, -77.055 }
988         };
989 
990         final double[][] expectedOut = new double[][] {
991                 { 17.794, -77.145 },
992                 { 17.736, -76.981 },
993                 { 17.715, -77.138 },
994                 { 18.153, -77.059 },
995                 { 18.232, -76.877 },
996                 { 18.373, -76.917 },
997                 { 17.871, -77.261 } // this is the point that was wrongly inside the intersection
998         };
999 
1000         SphericalPolygonsSet shape1 = buildSimpleZone(coordinates1);
1001         SphericalPolygonsSet shape2 = buildSimpleZone(coordinates2);
1002         Region<Sphere2D, S2Point, Circle, SubCircle> intersection =
1003                 order ?
1004                 new RegionFactory<Sphere2D, S2Point, Circle, SubCircle>().intersection(shape1.copySelf(), shape2.copySelf()) :
1005                 new RegionFactory<Sphere2D, S2Point, Circle, SubCircle>().intersection(shape2.copySelf(), shape1.copySelf());
1006 
1007         for (final double[] doubles : expectedIn) {
1008             Assertions.assertEquals(Location.INSIDE, intersection.checkPoint(s2Point(doubles[0], doubles[1])));
1009         }
1010 
1011         for (final double[] doubles : expectedOut) {
1012             Assertions.assertEquals(Location.OUTSIDE, intersection.checkPoint(s2Point(doubles[0], doubles[1])));
1013         }
1014 
1015     }
1016 
1017     private SubCircle create(Vector3D pole, Vector3D x, Vector3D y,
1018                              double tolerance, double ... limits) {
1019         RegionFactory<Sphere1D, S1Point, LimitAngle, SubLimitAngle> factory = new RegionFactory<>();
1020         Circle                           circle  = new Circle(pole, tolerance);
1021         Circle phased = Circle.getTransform(new Rotation(circle.getXAxis(), circle.getYAxis(), x, y)).apply(circle);
1022         ArcsSet set = (ArcsSet) factory.getComplement(new ArcsSet(tolerance));
1023         for (int i = 0; i < limits.length; i += 2) {
1024             set = (ArcsSet) factory.union(set, new ArcsSet(limits[i], limits[i + 1], tolerance));
1025         }
1026         return new SubCircle(phased, set);
1027     }
1028 
1029     private SphericalPolygonsSet buildSimpleZone(double[][] points) {
1030         final S2Point[] vertices = new S2Point[points.length];
1031         for (int i = 0; i < points.length; ++i) {
1032             vertices[i] = s2Point(points[i][0], points[i][1]);
1033         }
1034         return new SphericalPolygonsSet(1.0e-10, vertices);
1035     }
1036 
1037     private S2Point s2Point(double latitude, double longitude) {
1038         return new S2Point(FastMath.toRadians(longitude), FastMath.toRadians(90.0 - latitude));
1039     }
1040 
1041 }