1 /* Copyright 2002-2025 CS GROUP
2 * Licensed to CS GROUP (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS 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 * http://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 package org.orekit.models.earth.tessellation;
18
19 import java.util.ArrayList;
20 import java.util.Collection;
21 import java.util.IdentityHashMap;
22 import java.util.Iterator;
23 import java.util.LinkedList;
24 import java.util.List;
25 import java.util.Map;
26 import java.util.NoSuchElementException;
27 import java.util.Queue;
28
29 import org.hipparchus.exception.LocalizedCoreFormats;
30 import org.hipparchus.geometry.euclidean.threed.Vector3D;
31 import org.hipparchus.geometry.partitioning.BSPTree;
32 import org.hipparchus.geometry.partitioning.RegionFactory;
33 import org.hipparchus.geometry.partitioning.SubHyperplane;
34 import org.hipparchus.geometry.spherical.oned.ArcsSet;
35 import org.hipparchus.geometry.spherical.twod.Circle;
36 import org.hipparchus.geometry.spherical.twod.S2Point;
37 import org.hipparchus.geometry.spherical.twod.Sphere2D;
38 import org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet;
39 import org.hipparchus.geometry.spherical.twod.SubCircle;
40 import org.hipparchus.util.FastMath;
41 import org.hipparchus.util.MathUtils;
42 import org.orekit.bodies.GeodeticPoint;
43 import org.orekit.bodies.OneAxisEllipsoid;
44 import org.orekit.errors.OrekitException;
45 import org.orekit.errors.OrekitInternalError;
46
47 /** Class used to tessellate an interest zone on an ellipsoid in either
48 * {@link Tile tiles} or grids of {@link GeodeticPoint geodetic points}.
49 * <p>
50 * This class is typically used for Earth Observation missions, in order to
51 * create tiles or grids that may be used as the basis of visibility event
52 * detectors. Tiles are used when surface-related elements are needed, the
53 * tiles created completely cover the zone of interest. Grids are used when
54 * point-related elements are needed, the points created lie entirely within
55 * the zone of interest.
56 * </p>
57 * <p>
58 * One should note that as tessellation essentially creates a 2 dimensional
59 * almost Cartesian map, it can never perfectly fulfill geometrical dimensions
60 * because neither sphere nor ellipsoid are developable surfaces. This implies
61 * that the tesselation will always be distorted, and distortion increases as
62 * the size of the zone to be tessellated increases.
63 * </p>
64 * @author Luc Maisonobe
65 * @since 7.1
66 */
67 public class EllipsoidTessellator {
68
69 /** Safety limit to avoid infinite loops during tesselation due to numerical noise.
70 * @since 10.3.1
71 */
72 private static final int MAX_ITER = 1000;
73
74 /** Number of segments tiles sides are split into for tiles fine positioning. */
75 private final int quantization;
76
77 /** Aiming used for orienting tiles. */
78 private final TileAiming aiming;
79
80 /** Underlying ellipsoid. */
81 private final OneAxisEllipsoid ellipsoid;
82
83 /** Simple constructor.
84 * <p>
85 * The {@code quantization} parameter is used internally to adjust points positioning.
86 * For example when quantization is set to 4, a complete tile that has 4 corner points
87 * separated by the tile lengths will really be computed on a grid containing 25 points
88 * (5 rows of 5 points, as each side will be split in 4 segments, hence will have 5
89 * points). This quantization allows rough adjustment to balance margins around the
90 * zone of interest and improves geometric accuracy as the along and across directions
91 * are readjusted at each points.
92 * </p>
93 * <p>
94 * It is recommended to use at least 2 as the quantization parameter for tiling. The
95 * rationale is that using only 1 for quantization would imply all points used are tiles
96 * vertices, and hence would lead small zones to generate 4 tiles with a shared vertex
97 * inside the zone and the 4 tiles covering the four quadrants at North-West, North-East,
98 * South-East and South-West. A quantization value of at least 2 allows to shift the
99 * tiles so the center point is an inside point rather than a tile vertex, hence allowing
100 * a single tile to cover the small zone. A value even greater like 4 or 8 would allow even
101 * finer positioning to balance the tiles margins around the zone.
102 * </p>
103 * @param ellipsoid underlying ellipsoid
104 * @param aiming aiming used for orienting tiles
105 * @param quantization number of segments tiles sides are split into for tiles fine positioning
106 */
107 public EllipsoidTessellator(final OneAxisEllipsoid ellipsoid, final TileAiming aiming,
108 final int quantization) {
109 this.ellipsoid = ellipsoid;
110 this.aiming = aiming;
111 this.quantization = quantization;
112 }
113
114 /** Tessellate a zone of interest into tiles.
115 * <p>
116 * The created tiles will completely cover the zone of interest.
117 * </p>
118 * <p>
119 * The distance between a vertex at a tile corner and the vertex at the same corner
120 * in the next vertex are computed by subtracting the overlap width (resp. overlap length)
121 * from the full width (resp. full length). If for example the full width is specified to
122 * be 55 km and the overlap in width is specified to be +5 km, successive tiles would span
123 * as follows:
124 * </p>
125 * <ul>
126 * <li>tile 1 covering from 0 km to 55 km</li>
127 * <li>tile 2 covering from 50 km to 105 km</li>
128 * <li>tile 3 covering from 100 km to 155 km</li>
129 * <li>...</li>
130 * </ul>
131 * <p>
132 * In order to achieve the same 50 km step but using a 5 km gap instead of an overlap, one would
133 * need to specify the full width to be 45 km and the overlap to be -5 km. With these settings,
134 * successive tiles would span as follows:
135 * </p>
136 * <ul>
137 * <li>tile 1 covering from 0 km to 45 km</li>
138 * <li>tile 2 covering from 50 km to 95 km</li>
139 * <li>tile 3 covering from 100 km to 155 km</li>
140 * <li>...</li>
141 * </ul>
142 * @param zone zone of interest to tessellate
143 * @param fullWidth full tiles width as a distance on surface, including overlap (in meters)
144 * @param fullLength full tiles length as a distance on surface, including overlap (in meters)
145 * @param widthOverlap overlap between adjacent tiles (in meters), if negative the tiles
146 * will have a gap between each other instead of an overlap
147 * @param lengthOverlap overlap between adjacent tiles (in meters), if negative the tiles
148 * will have a gap between each other instead of an overlap
149 * @param truncateLastWidth if true, the first tiles strip will be started as close as
150 * possible to the zone of interest, and the last tiles strip will have its width reduced
151 * to also remain close to the zone of interest; if false all tiles strip will have the
152 * same {@code fullWidth} and they will be balanced around zone of interest
153 * @param truncateLastLength if true, the first tile in each strip will be started as close as
154 * possible to the zone of interest, and the last tile in each strip will have its length reduced
155 * to also remain close to the zone of interest; if false all tiles in each strip will have the
156 * same {@code fullLength} and they will be balanced around zone of interest
157 * @return a list of lists of tiles covering the zone of interest,
158 * each sub-list corresponding to a part not connected to the other
159 * parts (for example for islands)
160 */
161 public List<List<Tile>> tessellate(final SphericalPolygonsSet zone,
162 final double fullWidth, final double fullLength,
163 final double widthOverlap, final double lengthOverlap,
164 final boolean truncateLastWidth, final boolean truncateLastLength) {
165
166 final RegionFactory<Sphere2D, S2Point, Circle, SubCircle> factory = new RegionFactory<>();
167
168 final double splitWidth = (fullWidth - widthOverlap) / quantization;
169 final double splitLength = (fullLength - lengthOverlap) / quantization;
170 final Map<Mesh, List<Tile>> map = new IdentityHashMap<>();
171 SphericalPolygonsSet remaining = (SphericalPolygonsSet) zone.copySelf();
172 S2Point inside = getInsidePoint(remaining);
173
174 int count = 0;
175 while (inside != null) {
176
177 if (++count > MAX_ITER) {
178 throw new OrekitException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, MAX_ITER);
179 }
180
181 // find a mesh covering at least one connected part of the zone
182 final List<Mesh.Node> mergingSeeds = new ArrayList<>();
183 Mesh mesh = new Mesh(ellipsoid, zone, aiming, splitLength, splitWidth, inside);
184 mergingSeeds.add(mesh.getNode(0, 0));
185 List<Tile> tiles = null;
186 while (!mergingSeeds.isEmpty()) {
187
188 // expand the mesh around the seed
189 neighborExpandMesh(mesh, mergingSeeds, zone);
190
191 // extract the tiles from the mesh
192 // this further expands the mesh so tiles dimensions are multiples of quantization,
193 // hence it must be performed here before checking meshes independence
194 tiles = extractTiles(mesh, zone, lengthOverlap, widthOverlap, truncateLastWidth, truncateLastLength);
195
196 // check the mesh is independent from existing meshes
197 mergingSeeds.clear();
198 for (final Map.Entry<Mesh, List<Tile>> entry : map.entrySet()) {
199 if (!factory.intersection(mesh.getCoverage(), entry.getKey().getCoverage()).isEmpty()) {
200 // the meshes are not independent, they intersect each other!
201
202 // merge the two meshes together
203 mesh = mergeMeshes(mesh, entry.getKey(), mergingSeeds);
204 map.remove(entry.getKey());
205 break;
206
207 }
208 }
209
210 }
211
212 // remove the part of the zone covered by the mesh
213 remaining = (SphericalPolygonsSet) factory.difference(remaining, mesh.getCoverage());
214 inside = getInsidePoint(remaining);
215
216 map.put(mesh, tiles);
217
218 }
219
220 // concatenate the lists from the independent meshes
221 final List<List<Tile>> tilesLists = new ArrayList<>(map.size());
222 for (final Map.Entry<Mesh, List<Tile>> entry : map.entrySet()) {
223 tilesLists.add(entry.getValue());
224 }
225
226 return tilesLists;
227
228 }
229
230 /** Sample a zone of interest into a grid sample of {@link GeodeticPoint geodetic points}.
231 * <p>
232 * The created points will be entirely within the zone of interest.
233 * </p>
234 * @param zone zone of interest to sample
235 * @param width grid sample cells width as a distance on surface (in meters)
236 * @param length grid sample cells length as a distance on surface (in meters)
237 * @return a list of lists of points sampling the zone of interest,
238 * each sub-list corresponding to a part not connected to the other
239 * parts (for example for islands)
240 */
241 public List<List<GeodeticPoint>> sample(final SphericalPolygonsSet zone,
242 final double width, final double length) {
243
244 final RegionFactory<Sphere2D, S2Point, Circle, SubCircle> factory = new RegionFactory<>();
245 final double splitWidth = width / quantization;
246 final double splitLength = length / quantization;
247 final Map<Mesh, List<GeodeticPoint>> map = new IdentityHashMap<>();
248 SphericalPolygonsSet remaining = (SphericalPolygonsSet) zone.copySelf();
249 S2Point inside = getInsidePoint(remaining);
250
251 int count = 0;
252 while (inside != null) {
253
254 if (++count > MAX_ITER) {
255 throw new OrekitException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, MAX_ITER);
256 }
257
258 // find a mesh covering at least one connected part of the zone
259 final List<Mesh.Node> mergingSeeds = new ArrayList<>();
260 Mesh mesh = new Mesh(ellipsoid, zone, aiming, splitLength, splitWidth, inside);
261 mergingSeeds.add(mesh.getNode(0, 0));
262 List<GeodeticPoint> sample = null;
263 while (!mergingSeeds.isEmpty()) {
264
265 // expand the mesh around the seed
266 neighborExpandMesh(mesh, mergingSeeds, zone);
267
268 // extract the sample from the mesh
269 // this further expands the mesh so sample cells dimensions are multiples of quantization,
270 // hence it must be performed here before checking meshes independence
271 sample = extractSample(mesh);
272
273 // check the mesh is independent from existing meshes
274 mergingSeeds.clear();
275 for (final Map.Entry<Mesh, List<GeodeticPoint>> entry : map.entrySet()) {
276 if (!factory.intersection(mesh.getCoverage(), entry.getKey().getCoverage()).isEmpty()) {
277 // the meshes are not independent, they intersect each other!
278
279 // merge the two meshes together
280 mesh = mergeMeshes(mesh, entry.getKey(), mergingSeeds);
281 map.remove(entry.getKey());
282 break;
283
284 }
285 }
286
287 }
288
289 // remove the part of the zone covered by the mesh
290 remaining = (SphericalPolygonsSet) factory.difference(remaining, mesh.getCoverage());
291 inside = getInsidePoint(remaining);
292
293 map.put(mesh, sample);
294
295 }
296
297 // concatenate the lists from the independent meshes
298 final List<List<GeodeticPoint>> sampleLists = new ArrayList<>(map.size());
299 for (final Map.Entry<Mesh, List<GeodeticPoint>> entry : map.entrySet()) {
300 sampleLists.add(entry.getValue());
301 }
302
303 return sampleLists;
304
305 }
306
307 /** Get an inside point from a zone of interest.
308 * @param zone zone to mesh
309 * @return a point inside the zone or null if zone is empty or too thin
310 */
311 private S2Point getInsidePoint(final SphericalPolygonsSet zone) {
312
313 final InsidePointFinder finder = new InsidePointFinder(zone);
314 zone.getTree(false).visit(finder);
315 return finder.getInsidePoint();
316
317 }
318
319 /** Expand a mesh so it surrounds at least one connected part of a zone.
320 * <p>
321 * This part of mesh expansion is neighbors based. It includes the seed
322 * node neighbors, and their neighbors, and the neighbors of their
323 * neighbors until the path-connected sub-parts of the zone these nodes
324 * belong to are completely surrounded by the mesh taxicab boundary.
325 * </p>
326 * @param mesh mesh to expand
327 * @param seeds seed nodes (already in the mesh) from which to start expansion
328 * @param zone zone to mesh
329 */
330 private void neighborExpandMesh(final Mesh mesh, final Collection<Mesh.Node> seeds,
331 final SphericalPolygonsSet zone) {
332
333 // mesh expansion loop
334 boolean expanding = true;
335 final Queue<Mesh.Node> newNodes = new LinkedList<>(seeds);
336 int count = 0;
337 while (expanding) {
338
339 if (++count > MAX_ITER) {
340 throw new OrekitException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, MAX_ITER);
341 }
342
343 // first expansion step: set up the mesh so that all its
344 // inside nodes are completely surrounded by at least
345 // one layer of outside nodes
346 while (!newNodes.isEmpty()) {
347
348 // retrieve an active node
349 final Mesh.Node node = newNodes.remove();
350
351 if (node.isInside()) {
352 // the node is inside the zone, the mesh must contain its 8 neighbors
353 addAllNeighborsIfNeeded(node, mesh, newNodes);
354 }
355
356 }
357
358 // second expansion step: check if the loop of outside nodes
359 // completely surrounds the zone, i.e. there are no peaks
360 // pointing out of the loop between two nodes
361 expanding = false;
362 final List<Mesh.Node> boundary = mesh.getTaxicabBoundary(false);
363 if (boundary.size() > 1) {
364 Mesh.Node previous = boundary.get(boundary.size() - 1);
365 for (final Mesh.Node node : boundary) {
366 if (meetInside(previous.getS2P(), node.getS2P(), zone)) {
367 // part of the mesh boundary is still inside the zone!
368 // the mesh must be expanded again
369 addAllNeighborsIfNeeded(previous, mesh, newNodes);
370 addAllNeighborsIfNeeded(node, mesh, newNodes);
371 expanding = true;
372 }
373 previous = node;
374 }
375 }
376
377 }
378
379 }
380
381 /** Extract tiles from a mesh.
382 * @param mesh mesh from which tiles should be extracted
383 * @param zone zone covered by the mesh
384 * @param lengthOverlap overlap between adjacent tiles
385 * @param widthOverlap overlap between adjacent tiles
386 * @param truncateLastWidth true if we can reduce last tile width
387 * @param truncateLastLength true if we can reduce last tile length
388 * @return extracted tiles
389 */
390 private List<Tile> extractTiles(final Mesh mesh, final SphericalPolygonsSet zone,
391 final double lengthOverlap, final double widthOverlap,
392 final boolean truncateLastWidth, final boolean truncateLastLength) {
393
394 final List<Tile> tiles = new ArrayList<>();
395 final List<RangePair> rangePairs = new ArrayList<>();
396
397 final int minAcross = mesh.getMinAcrossIndex();
398 final int maxAcross = mesh.getMaxAcrossIndex();
399 for (Range acrossPair : nodesIndices(minAcross, maxAcross, truncateLastWidth)) {
400
401 int minAlong = mesh.getMaxAlongIndex() + 1;
402 int maxAlong = mesh.getMinAlongIndex() - 1;
403 for (int c = acrossPair.lower; c <= acrossPair.upper; ++c) {
404 minAlong = FastMath.min(minAlong, mesh.getMinAlongIndex(c));
405 maxAlong = FastMath.max(maxAlong, mesh.getMaxAlongIndex(c));
406 }
407
408 for (Range alongPair : nodesIndices(minAlong, maxAlong, truncateLastLength)) {
409
410 // get the base vertex nodes
411 final Mesh.Node node0 = mesh.addNode(alongPair.lower, acrossPair.lower);
412 final Mesh.Node node1 = mesh.addNode(alongPair.upper, acrossPair.lower);
413 final Mesh.Node node2 = mesh.addNode(alongPair.upper, acrossPair.upper);
414 final Mesh.Node node3 = mesh.addNode(alongPair.lower, acrossPair.upper);
415
416 // apply tile overlap
417 final S2Point s2p0 = node0.move(new Vector3D(-0.5 * lengthOverlap, node0.getAlong(),
418 -0.5 * widthOverlap, node0.getAcross()));
419 final S2Point s2p1 = node1.move(new Vector3D(+0.5 * lengthOverlap, node1.getAlong(),
420 -0.5 * widthOverlap, node1.getAcross()));
421 final S2Point s2p2 = node2.move(new Vector3D(+0.5 * lengthOverlap, node2.getAlong(),
422 +0.5 * widthOverlap, node2.getAcross()));
423 final S2Point s2p3 = node3.move(new Vector3D(-0.5 * lengthOverlap, node2.getAlong(),
424 +0.5 * widthOverlap, node2.getAcross()));
425
426 // create a quadrilateral region corresponding to the candidate tile
427 final SphericalPolygonsSet quadrilateral =
428 new SphericalPolygonsSet(zone.getTolerance(), s2p0, s2p1, s2p2, s2p3);
429
430 if (!new RegionFactory<Sphere2D, S2Point, Circle, SubCircle>().
431 intersection(zone.copySelf(), quadrilateral).isEmpty()) {
432
433 // the tile does cover part of the zone, it contributes to the tessellation
434 tiles.add(new Tile(toGeodetic(s2p0), toGeodetic(s2p1), toGeodetic(s2p2), toGeodetic(s2p3)));
435 rangePairs.add(new RangePair(acrossPair, alongPair));
436
437 }
438
439 }
440 }
441
442 // ensure the taxicab boundary follows the built tile sides
443 // this is done outside of the previous loop in order
444 // to avoid one tile changing the min/max indices of the
445 // neighboring tile as they share some nodes that will be enabled here
446 for (final RangePair rangePair : rangePairs) {
447 for (int c = rangePair.across.lower; c < rangePair.across.upper; ++c) {
448 mesh.addNode(rangePair.along.lower, c + 1).setEnabled();
449 mesh.addNode(rangePair.along.upper, c).setEnabled();
450 }
451 for (int l = rangePair.along.lower; l < rangePair.along.upper; ++l) {
452 mesh.addNode(l, rangePair.across.lower).setEnabled();
453 mesh.addNode(l + 1, rangePair.across.upper).setEnabled();
454 }
455 }
456
457 return tiles;
458
459 }
460
461 /**
462 * Extract a sample of points from a mesh.
463 *
464 * @param mesh mesh from which grid should be extracted
465 * @return extracted grid
466 */
467 private List<GeodeticPoint> extractSample(final Mesh mesh) {
468
469 // find how to select sample points taking quantization into account
470 // to have the largest possible number of points while still
471 // being inside the zone of interest
472 int selectedAcrossModulus = -1;
473 int selectedAlongModulus = -1;
474 int selectedCount = -1;
475 for (int acrossModulus = 0; acrossModulus < quantization; ++acrossModulus) {
476 for (int alongModulus = 0; alongModulus < quantization; ++alongModulus) {
477
478 // count how many points would be selected for the current modulus
479 int count = 0;
480 for (int across = mesh.getMinAcrossIndex() + acrossModulus;
481 across <= mesh.getMaxAcrossIndex();
482 across += quantization) {
483 for (int along = mesh.getMinAlongIndex() + alongModulus;
484 along <= mesh.getMaxAlongIndex();
485 along += quantization) {
486 final Mesh.Node node = mesh.getNode(along, across);
487 if (node != null && node.isInside()) {
488 ++count;
489 }
490 }
491 }
492
493 if (count > selectedCount) {
494 // current modulus are better than the selected ones
495 selectedAcrossModulus = acrossModulus;
496 selectedAlongModulus = alongModulus;
497 selectedCount = count;
498 }
499 }
500 }
501
502 // extract the sample points
503 final List<GeodeticPoint> sample = new ArrayList<>(selectedCount);
504 for (int across = mesh.getMinAcrossIndex() + selectedAcrossModulus;
505 across <= mesh.getMaxAcrossIndex();
506 across += quantization) {
507 for (int along = mesh.getMinAlongIndex() + selectedAlongModulus;
508 along <= mesh.getMaxAlongIndex();
509 along += quantization) {
510 final Mesh.Node node = mesh.getNode(along, across);
511 if (node != null && node.isInside()) {
512 sample.add(toGeodetic(node.getS2P()));
513 }
514 }
515 }
516
517 return sample;
518
519 }
520
521 /** Merge two meshes together.
522 * @param mesh1 first mesh
523 * @param mesh2 second mesh
524 * @param mergingSeeds collection where to put the nodes created during the merge
525 * @return merged mesh (really one of the instances)
526 */
527 private Mesh mergeMeshes(final Mesh mesh1, final Mesh mesh2,
528 final Collection<Mesh.Node> mergingSeeds) {
529
530 // select the way merge will be performed
531 final Mesh larger;
532 final Mesh smaller;
533 if (mesh1.getNumberOfNodes() >= mesh2.getNumberOfNodes()) {
534 // the larger new mesh should absorb the smaller existing mesh
535 larger = mesh1;
536 smaller = mesh2;
537 } else {
538 // the larger existing mesh should absorb the smaller new mesh
539 larger = mesh2;
540 smaller = mesh1;
541 }
542
543 // prepare seed nodes for next iteration
544 for (final Mesh.Node insideNode : smaller.getInsideNodes()) {
545
546 // beware we cannot reuse the node itself as the two meshes are not aligned!
547 // we have to create new nodes around the previous location
548 Mesh.Node node = larger.getClosestExistingNode(insideNode.getV());
549
550 while (estimateAlongMotion(node, insideNode.getV()) > +mesh1.getAlongGap()) {
551 // the node is before desired index in the along direction
552 // we need to create intermediates nodes up to the desired index
553 node = larger.addNode(node.getAlongIndex() + 1, node.getAcrossIndex());
554 }
555
556 while (estimateAlongMotion(node, insideNode.getV()) < -mesh1.getAlongGap()) {
557 // the node is after desired index in the along direction
558 // we need to create intermediates nodes up to the desired index
559 node = larger.addNode(node.getAlongIndex() - 1, node.getAcrossIndex());
560 }
561
562 while (estimateAcrossMotion(node, insideNode.getV()) > +mesh1.getAcrossGap()) {
563 // the node is before desired index in the across direction
564 // we need to create intermediates nodes up to the desired index
565 node = larger.addNode(node.getAlongIndex(), node.getAcrossIndex() + 1);
566 }
567
568 while (estimateAcrossMotion(node, insideNode.getV()) < -mesh1.getAcrossGap()) {
569 // the node is after desired index in the across direction
570 // we need to create intermediates nodes up to the desired index
571 node = larger.addNode(node.getAlongIndex(), node.getAcrossIndex() - 1);
572 }
573
574 // now we are close to the inside node,
575 // make sure the four surrounding nodes are available
576 final int otherAlong = (estimateAlongMotion(node, insideNode.getV()) < 0.0) ?
577 node.getAlongIndex() - 1 : node.getAlongIndex() + 1;
578 final int otherAcross = (estimateAcrossMotion(node, insideNode.getV()) < 0.0) ?
579 node.getAcrossIndex() - 1 : node.getAcrossIndex() + 1;
580 addNode(node.getAlongIndex(), node.getAcrossIndex(), larger, mergingSeeds);
581 addNode(node.getAlongIndex(), otherAcross, larger, mergingSeeds);
582 addNode(otherAlong, node.getAcrossIndex(), larger, mergingSeeds);
583 addNode(otherAlong, otherAcross, larger, mergingSeeds);
584
585 }
586
587 return larger;
588
589 }
590
591 /** Ensure all 8 neighbors of a node are in the mesh.
592 * @param base base node
593 * @param mesh complete mesh containing nodes
594 * @param newNodes queue where new node must be put
595 */
596 private void addAllNeighborsIfNeeded(final Mesh.Node base, final Mesh mesh,
597 final Collection<Mesh.Node> newNodes) {
598 addNode(base.getAlongIndex() - 1, base.getAcrossIndex() - 1, mesh, newNodes);
599 addNode(base.getAlongIndex() - 1, base.getAcrossIndex(), mesh, newNodes);
600 addNode(base.getAlongIndex() - 1, base.getAcrossIndex() + 1, mesh, newNodes);
601 addNode(base.getAlongIndex(), base.getAcrossIndex() - 1, mesh, newNodes);
602 addNode(base.getAlongIndex(), base.getAcrossIndex() + 1, mesh, newNodes);
603 addNode(base.getAlongIndex() + 1, base.getAcrossIndex() - 1, mesh, newNodes);
604 addNode(base.getAlongIndex() + 1, base.getAcrossIndex(), mesh, newNodes);
605 addNode(base.getAlongIndex() + 1, base.getAcrossIndex() + 1, mesh, newNodes);
606 }
607
608 /** Add a node to a mesh if not already present.
609 * @param alongIndex index in the along direction
610 * @param acrossIndex index in the across direction
611 * @param mesh complete mesh containing nodes
612 * @param newNodes queue where new node must be put
613 */
614 private void addNode(final int alongIndex, final int acrossIndex,
615 final Mesh mesh, final Collection<Mesh.Node> newNodes) {
616
617 final Mesh.Node node = mesh.addNode(alongIndex, acrossIndex);
618
619 if (!node.isEnabled()) {
620 // enable the node
621 node.setEnabled();
622 newNodes.add(node);
623 }
624
625 }
626
627 /** Convert a point on the unit 2-sphere to geodetic coordinates.
628 * @param point point on the unit 2-sphere
629 * @return geodetic point (arbitrarily set at altitude 0)
630 */
631 protected GeodeticPoint toGeodetic(final S2Point point) {
632 return new GeodeticPoint(0.5 * FastMath.PI - point.getPhi(), point.getTheta(), 0.0);
633 }
634
635 /** Build a simple zone (connected zone without holes).
636 * <p>
637 * In order to build more complex zones (not connected or with
638 * holes), the user should directly call Hipparchus
639 * {@link SphericalPolygonsSet} constructors and
640 * {@link RegionFactory region factory} if set operations
641 * are needed (union, intersection, difference ...).
642 * </p>
643 * <p>
644 * Take care that the vertices boundary points must be given <em>counterclockwise</em>.
645 * Using the wrong order defines the complementary of the real zone,
646 * and will often result in tessellation failure as the zone is too
647 * wide.
648 * </p>
649 * @param tolerance angular separation below which points are considered
650 * equal (typically 1.0e-10)
651 * @param points vertices of the boundary, in <em>counterclockwise</em>
652 * order, each point being a two-elements arrays with latitude at index 0
653 * and longitude at index 1
654 * @return a zone defined on the unit 2-sphere
655 */
656 public static SphericalPolygonsSet buildSimpleZone(final double tolerance,
657 final double[]... points) {
658 final S2Point[] vertices = new S2Point[points.length];
659 for (int i = 0; i < points.length; ++i) {
660 vertices[i] = new S2Point(points[i][1], 0.5 * FastMath.PI - points[i][0]);
661 }
662 return new SphericalPolygonsSet(tolerance, vertices);
663 }
664
665 /** Build a simple zone (connected zone without holes).
666 * <p>
667 * In order to build more complex zones (not connected or with
668 * holes), the user should directly call Hipparchus
669 * {@link SphericalPolygonsSet} constructors and
670 * {@link RegionFactory region factory} if set operations
671 * are needed (union, intersection, difference ...).
672 * </p>
673 * <p>
674 * Take care that the vertices boundary points must be given <em>counterclockwise</em>.
675 * Using the wrong order defines the complementary of the real zone,
676 * and will often result in tessellation failure as the zone is too
677 * wide.
678 * </p>
679 * @param tolerance angular separation below which points are considered
680 * equal (typically 1.0e-10)
681 * @param points vertices of the boundary, in <em>counterclockwise</em>
682 * order
683 * @return a zone defined on the unit 2-sphere
684 */
685 public static SphericalPolygonsSet buildSimpleZone(final double tolerance,
686 final GeodeticPoint... points) {
687 final S2Point[] vertices = new S2Point[points.length];
688 for (int i = 0; i < points.length; ++i) {
689 vertices[i] = new S2Point(points[i].getLongitude(),
690 0.5 * FastMath.PI - points[i].getLatitude());
691 }
692 return new SphericalPolygonsSet(tolerance, vertices);
693 }
694
695 /** Estimate an approximate motion in the along direction.
696 * @param start node at start of motion
697 * @param end desired point at end of motion
698 * @return approximate motion in the along direction
699 */
700 private double estimateAlongMotion(final Mesh.Node start, final Vector3D end) {
701 return Vector3D.dotProduct(start.getAlong(), end.subtract(start.getV()));
702 }
703
704 /** Estimate an approximate motion in the across direction.
705 * @param start node at start of motion
706 * @param end desired point at end of motion
707 * @return approximate motion in the across direction
708 */
709 private double estimateAcrossMotion(final Mesh.Node start, final Vector3D end) {
710 return Vector3D.dotProduct(start.getAcross(), end.subtract(start.getV()));
711 }
712
713 /** Check if an arc meets the inside of a zone.
714 * @param s1 first point
715 * @param s2 second point
716 * @param zone zone to check arc against
717 * @return true if the arc meets the inside of the zone
718 */
719 private boolean meetInside(final S2Point s1, final S2Point s2,
720 final SphericalPolygonsSet zone) {
721 final Circle circle = new Circle(s1, s2, zone.getTolerance());
722 final double alpha1 = circle.toSubSpace(s1).getAlpha();
723 final double alpha2 = MathUtils.normalizeAngle(circle.toSubSpace(s2).getAlpha(),
724 alpha1 + FastMath.PI);
725 final SubCircle sub = new SubCircle(circle,
726 new ArcsSet(alpha1, alpha2, zone.getTolerance()));
727 return recurseMeetInside(zone.getTree(false), sub);
728
729 }
730
731 /** Check if an arc meets the inside of a zone.
732 * <p>
733 * This method is heavily based on the Characterization class from
734 * Hipparchus library, also distributed under the terms
735 * of the Apache Software License V2.
736 * </p>
737 * @param node spherical zone node
738 * @param sub arc to characterize
739 * @return true if the arc meets the inside of the zone
740 */
741 private boolean recurseMeetInside(final BSPTree<Sphere2D, S2Point, Circle, SubCircle> node,
742 final SubHyperplane<Sphere2D, S2Point, Circle, SubCircle> sub) {
743
744 if (node.getCut() == null) {
745 // we have reached a leaf node
746 if (sub.isEmpty()) {
747 return false;
748 } else {
749 return (Boolean) node.getAttribute();
750 }
751 } else {
752 final Circle hyperplane = node.getCut().getHyperplane();
753 final SubHyperplane.SplitSubHyperplane<Sphere2D, S2Point, Circle, SubCircle> split = sub.split(hyperplane);
754 switch (split.getSide()) {
755 case PLUS:
756 return recurseMeetInside(node.getPlus(), sub);
757 case MINUS:
758 return recurseMeetInside(node.getMinus(), sub);
759 case BOTH:
760 if (recurseMeetInside(node.getPlus(), split.getPlus())) {
761 return true;
762 } else {
763 return recurseMeetInside(node.getMinus(), split.getMinus());
764 }
765 default:
766 // this should not happen
767 throw new OrekitInternalError(null);
768 }
769 }
770 }
771
772 /** Get an iterator over mesh nodes indices.
773 * @param minIndex minimum node index
774 * @param maxIndex maximum node index
775 * @param truncateLast true if we can reduce last tile
776 * @return iterator over mesh nodes indices
777 */
778 private Iterable<Range> nodesIndices(final int minIndex, final int maxIndex, final boolean truncateLast) {
779
780 final int first;
781 if (truncateLast) {
782
783 // truncate last tile rather than balance tiles around the zone of interest
784 first = minIndex;
785
786 } else {
787
788 // balance tiles around the zone of interest rather than truncate last tile
789
790 // number of tiles needed to cover the full indices range
791 final int range = maxIndex - minIndex;
792 final int nbTiles = (range + quantization - 1) / quantization;
793
794 // extra nodes that must be added to complete the tiles
795 final int extraNodes = nbTiles * quantization - range;
796
797 // balance the extra nodes before min index and after maxIndex
798 final int extraBefore = (extraNodes + 1) / 2;
799
800 first = minIndex - extraBefore;
801
802 }
803
804 return new Iterable<Range>() {
805
806 /** {@inheritDoc} */
807 @Override
808 public Iterator<Range> iterator() {
809 return new Iterator<Range>() {
810
811 private int nextLower = first;
812
813 /** {@inheritDoc} */
814 @Override
815 public boolean hasNext() {
816 return nextLower < maxIndex;
817 }
818
819 /** {@inheritDoc} */
820 @Override
821 public Range next() {
822
823 if (nextLower >= maxIndex) {
824 throw new NoSuchElementException();
825 }
826 final int lower = nextLower;
827
828 nextLower += quantization;
829 if (truncateLast && nextLower > maxIndex) {
830 // truncate last tile
831 nextLower = maxIndex;
832 }
833
834 return new Range(lower, nextLower);
835
836 }
837
838 /** {@inheritDoc} */
839 @Override
840 public void remove() {
841 throw new UnsupportedOperationException();
842 }
843
844 };
845 }
846 };
847
848 }
849
850 /** Local class for a range of indices to be used for building a tile. */
851 private static class Range {
852
853 /** Lower index. */
854 private final int lower;
855
856 /** Upper index. */
857 private final int upper;
858
859 /** Simple constructor.
860 * @param lower lower index
861 * @param upper upper index
862 */
863 Range(final int lower, final int upper) {
864 this.lower = lower;
865 this.upper = upper;
866 }
867
868 }
869
870 /** Local class for a pair of ranges of indices to be used for building a tile. */
871 private static class RangePair {
872
873 /** Across range. */
874 private final Range across;
875
876 /** Along range. */
877 private final Range along;
878
879 /** Simple constructor.
880 * @param across across range
881 * @param along along range
882 */
883 RangePair(final Range across, final Range along) {
884 this.across = across;
885 this.along = along;
886 }
887
888 }
889
890 }