1 /* Copyright 2002-2024 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.Hyperplane;
33 import org.hipparchus.geometry.partitioning.RegionFactory;
34 import org.hipparchus.geometry.partitioning.SubHyperplane;
35 import org.hipparchus.geometry.spherical.oned.ArcsSet;
36 import org.hipparchus.geometry.spherical.twod.Circle;
37 import org.hipparchus.geometry.spherical.twod.S2Point;
38 import org.hipparchus.geometry.spherical.twod.Sphere2D;
39 import org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet;
40 import org.hipparchus.geometry.spherical.twod.SubCircle;
41 import org.hipparchus.util.FastMath;
42 import org.hipparchus.util.MathUtils;
43 import org.orekit.bodies.GeodeticPoint;
44 import org.orekit.bodies.OneAxisEllipsoid;
45 import org.orekit.errors.OrekitException;
46 import org.orekit.errors.OrekitInternalError;
47
48 /** Class used to tessellate an interest zone on an ellipsoid in either
49 * {@link Tile tiles} or grids of {@link GeodeticPoint geodetic points}.
50 * <p>
51 * This class is typically used for Earth Observation missions, in order to
52 * create tiles or grids that may be used as the basis of visibility event
53 * detectors. Tiles are used when surface-related elements are needed, the
54 * tiles created completely cover the zone of interest. Grids are used when
55 * point-related elements are needed, the points created lie entirely within
56 * the zone of interest.
57 * </p>
58 * <p>
59 * One should note that as tessellation essentially creates a 2 dimensional
60 * almost Cartesian map, it can never perfectly fulfill geometrical dimensions
61 * because neither sphere nor ellipsoid are developable surfaces. This implies
62 * that the tesselation will always be distorted, and distortion increases as
63 * the size of the zone to be tessellated increases.
64 * </p>
65 * @author Luc Maisonobe
66 * @since 7.1
67 */
68 public class EllipsoidTessellator {
69
70 /** Safety limit to avoid infinite loops during tesselation due to numerical noise.
71 * @since 10.3.1
72 */
73 private static final int MAX_ITER = 1000;
74
75 /** Number of segments tiles sides are split into for tiles fine positioning. */
76 private final int quantization;
77
78 /** Aiming used for orienting tiles. */
79 private final TileAiming aiming;
80
81 /** Underlying ellipsoid. */
82 private final OneAxisEllipsoid ellipsoid;
83
84 /** Simple constructor.
85 * <p>
86 * The {@code quantization} parameter is used internally to adjust points positioning.
87 * For example when quantization is set to 4, a complete tile that has 4 corner points
88 * separated by the tile lengths will really be computed on a grid containing 25 points
89 * (5 rows of 5 points, as each side will be split in 4 segments, hence will have 5
90 * points). This quantization allows rough adjustment to balance margins around the
91 * zone of interest and improves geometric accuracy as the along and across directions
92 * are readjusted at each points.
93 * </p>
94 * <p>
95 * It is recommended to use at least 2 as the quantization parameter for tiling. The
96 * rationale is that using only 1 for quantization would imply all points used are tiles
97 * vertices, and hence would lead small zones to generate 4 tiles with a shared vertex
98 * inside the zone and the 4 tiles covering the four quadrants at North-West, North-East,
99 * South-East and South-West. A quantization value of at least 2 allows to shift the
100 * tiles so the center point is an inside point rather than a tile vertex, hence allowing
101 * a single tile to cover the small zone. A value even greater like 4 or 8 would allow even
102 * finer positioning to balance the tiles margins around the zone.
103 * </p>
104 * @param ellipsoid underlying ellipsoid
105 * @param aiming aiming used for orienting tiles
106 * @param quantization number of segments tiles sides are split into for tiles fine positioning
107 */
108 public EllipsoidTessellator(final OneAxisEllipsoid ellipsoid, final TileAiming aiming,
109 final int quantization) {
110 this.ellipsoid = ellipsoid;
111 this.aiming = aiming;
112 this.quantization = quantization;
113 }
114
115 /** Tessellate a zone of interest into tiles.
116 * <p>
117 * The created tiles will completely cover the zone of interest.
118 * </p>
119 * <p>
120 * The distance between a vertex at a tile corner and the vertex at the same corner
121 * in the next vertex are computed by subtracting the overlap width (resp. overlap length)
122 * from the full width (resp. full length). If for example the full width is specified to
123 * be 55 km and the overlap in width is specified to be +5 km, successive tiles would span
124 * as follows:
125 * </p>
126 * <ul>
127 * <li>tile 1 covering from 0 km to 55 km</li>
128 * <li>tile 2 covering from 50 km to 105 km</li>
129 * <li>tile 3 covering from 100 km to 155 km</li>
130 * <li>...</li>
131 * </ul>
132 * <p>
133 * In order to achieve the same 50 km step but using a 5 km gap instead of an overlap, one would
134 * need to specify the full width to be 45 km and the overlap to be -5 km. With these settings,
135 * successive tiles would span as follows:
136 * </p>
137 * <ul>
138 * <li>tile 1 covering from 0 km to 45 km</li>
139 * <li>tile 2 covering from 50 km to 95 km</li>
140 * <li>tile 3 covering from 100 km to 155 km</li>
141 * <li>...</li>
142 * </ul>
143 * @param zone zone of interest to tessellate
144 * @param fullWidth full tiles width as a distance on surface, including overlap (in meters)
145 * @param fullLength full tiles length as a distance on surface, including overlap (in meters)
146 * @param widthOverlap overlap between adjacent tiles (in meters), if negative the tiles
147 * will have a gap between each other instead of an overlap
148 * @param lengthOverlap overlap between adjacent tiles (in meters), if negative the tiles
149 * will have a gap between each other instead of an overlap
150 * @param truncateLastWidth if true, the first tiles strip will be started as close as
151 * possible to the zone of interest, and the last tiles strip will have its width reduced
152 * to also remain close to the zone of interest; if false all tiles strip will have the
153 * same {@code fullWidth} and they will be balanced around zone of interest
154 * @param truncateLastLength if true, the first tile in each strip will be started as close as
155 * possible to the zone of interest, and the last tile in each strip will have its length reduced
156 * to also remain close to the zone of interest; if false all tiles in each strip will have the
157 * same {@code fullLength} and they will be balanced around zone of interest
158 * @return a list of lists of tiles covering the zone of interest,
159 * each sub-list corresponding to a part not connected to the other
160 * parts (for example for islands)
161 */
162 public List<List<Tile>> tessellate(final SphericalPolygonsSet zone,
163 final double fullWidth, final double fullLength,
164 final double widthOverlap, final double lengthOverlap,
165 final boolean truncateLastWidth, final boolean truncateLastLength) {
166
167 final double splitWidth = (fullWidth - widthOverlap) / quantization;
168 final double splitLength = (fullLength - lengthOverlap) / quantization;
169 final Map<Mesh, List<Tile>> map = new IdentityHashMap<>();
170 final RegionFactory<Sphere2D> factory = new RegionFactory<>();
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 double splitWidth = width / quantization;
245 final double splitLength = length / quantization;
246 final Map<Mesh, List<GeodeticPoint>> map = new IdentityHashMap<>();
247 final RegionFactory<Sphere2D> factory = new RegionFactory<>();
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>().intersection(zone.copySelf(), quadrilateral).isEmpty()) {
431
432 // the tile does cover part of the zone, it contributes to the tessellation
433 tiles.add(new Tile(toGeodetic(s2p0), toGeodetic(s2p1), toGeodetic(s2p2), toGeodetic(s2p3)));
434 rangePairs.add(new RangePair(acrossPair, alongPair));
435
436 }
437
438 }
439 }
440
441 // ensure the taxicab boundary follows the built tile sides
442 // this is done outside of the previous loop in order
443 // to avoid one tile changing the min/max indices of the
444 // neighboring tile as they share some nodes that will be enabled here
445 for (final RangePair rangePair : rangePairs) {
446 for (int c = rangePair.across.lower; c < rangePair.across.upper; ++c) {
447 mesh.addNode(rangePair.along.lower, c + 1).setEnabled();
448 mesh.addNode(rangePair.along.upper, c).setEnabled();
449 }
450 for (int l = rangePair.along.lower; l < rangePair.along.upper; ++l) {
451 mesh.addNode(l, rangePair.across.lower).setEnabled();
452 mesh.addNode(l + 1, rangePair.across.upper).setEnabled();
453 }
454 }
455
456 return tiles;
457
458 }
459
460 /**
461 * Extract a sample of points from a mesh.
462 *
463 * @param mesh mesh from which grid should be extracted
464 * @return extracted grid
465 */
466 private List<GeodeticPoint> extractSample(final Mesh mesh) {
467
468 // find how to select sample points taking quantization into account
469 // to have the largest possible number of points while still
470 // being inside the zone of interest
471 int selectedAcrossModulus = -1;
472 int selectedAlongModulus = -1;
473 int selectedCount = -1;
474 for (int acrossModulus = 0; acrossModulus < quantization; ++acrossModulus) {
475 for (int alongModulus = 0; alongModulus < quantization; ++alongModulus) {
476
477 // count how many points would be selected for the current modulus
478 int count = 0;
479 for (int across = mesh.getMinAcrossIndex() + acrossModulus;
480 across <= mesh.getMaxAcrossIndex();
481 across += quantization) {
482 for (int along = mesh.getMinAlongIndex() + alongModulus;
483 along <= mesh.getMaxAlongIndex();
484 along += quantization) {
485 final Mesh.Node node = mesh.getNode(along, across);
486 if (node != null && node.isInside()) {
487 ++count;
488 }
489 }
490 }
491
492 if (count > selectedCount) {
493 // current modulus are better than the selected ones
494 selectedAcrossModulus = acrossModulus;
495 selectedAlongModulus = alongModulus;
496 selectedCount = count;
497 }
498 }
499 }
500
501 // extract the sample points
502 final List<GeodeticPoint> sample = new ArrayList<>(selectedCount);
503 for (int across = mesh.getMinAcrossIndex() + selectedAcrossModulus;
504 across <= mesh.getMaxAcrossIndex();
505 across += quantization) {
506 for (int along = mesh.getMinAlongIndex() + selectedAlongModulus;
507 along <= mesh.getMaxAlongIndex();
508 along += quantization) {
509 final Mesh.Node node = mesh.getNode(along, across);
510 if (node != null && node.isInside()) {
511 sample.add(toGeodetic(node.getS2P()));
512 }
513 }
514 }
515
516 return sample;
517
518 }
519
520 /** Merge two meshes together.
521 * @param mesh1 first mesh
522 * @param mesh2 second mesh
523 * @param mergingSeeds collection where to put the nodes created during the merge
524 * @return merged mesh (really one of the instances)
525 */
526 private Mesh mergeMeshes(final Mesh mesh1, final Mesh mesh2,
527 final Collection<Mesh.Node> mergingSeeds) {
528
529 // select the way merge will be performed
530 final Mesh larger;
531 final Mesh smaller;
532 if (mesh1.getNumberOfNodes() >= mesh2.getNumberOfNodes()) {
533 // the larger new mesh should absorb the smaller existing mesh
534 larger = mesh1;
535 smaller = mesh2;
536 } else {
537 // the larger existing mesh should absorb the smaller new mesh
538 larger = mesh2;
539 smaller = mesh1;
540 }
541
542 // prepare seed nodes for next iteration
543 for (final Mesh.Node insideNode : smaller.getInsideNodes()) {
544
545 // beware we cannot reuse the node itself as the two meshes are not aligned!
546 // we have to create new nodes around the previous location
547 Mesh.Node node = larger.getClosestExistingNode(insideNode.getV());
548
549 while (estimateAlongMotion(node, insideNode.getV()) > +mesh1.getAlongGap()) {
550 // the node is before desired index in the along direction
551 // we need to create intermediates nodes up to the desired index
552 node = larger.addNode(node.getAlongIndex() + 1, node.getAcrossIndex());
553 }
554
555 while (estimateAlongMotion(node, insideNode.getV()) < -mesh1.getAlongGap()) {
556 // the node is after desired index in the along direction
557 // we need to create intermediates nodes up to the desired index
558 node = larger.addNode(node.getAlongIndex() - 1, node.getAcrossIndex());
559 }
560
561 while (estimateAcrossMotion(node, insideNode.getV()) > +mesh1.getAcrossGap()) {
562 // the node is before desired index in the across direction
563 // we need to create intermediates nodes up to the desired index
564 node = larger.addNode(node.getAlongIndex(), node.getAcrossIndex() + 1);
565 }
566
567 while (estimateAcrossMotion(node, insideNode.getV()) < -mesh1.getAcrossGap()) {
568 // the node is after desired index in the across direction
569 // we need to create intermediates nodes up to the desired index
570 node = larger.addNode(node.getAlongIndex(), node.getAcrossIndex() - 1);
571 }
572
573 // now we are close to the inside node,
574 // make sure the four surrounding nodes are available
575 final int otherAlong = (estimateAlongMotion(node, insideNode.getV()) < 0.0) ?
576 node.getAlongIndex() - 1 : node.getAlongIndex() + 1;
577 final int otherAcross = (estimateAcrossMotion(node, insideNode.getV()) < 0.0) ?
578 node.getAcrossIndex() - 1 : node.getAcrossIndex() + 1;
579 addNode(node.getAlongIndex(), node.getAcrossIndex(), larger, mergingSeeds);
580 addNode(node.getAlongIndex(), otherAcross, larger, mergingSeeds);
581 addNode(otherAlong, node.getAcrossIndex(), larger, mergingSeeds);
582 addNode(otherAlong, otherAcross, larger, mergingSeeds);
583
584 }
585
586 return larger;
587
588 }
589
590 /** Ensure all 8 neighbors of a node are in the mesh.
591 * @param base base node
592 * @param mesh complete mesh containing nodes
593 * @param newNodes queue where new node must be put
594 */
595 private void addAllNeighborsIfNeeded(final Mesh.Node base, final Mesh mesh,
596 final Collection<Mesh.Node> newNodes) {
597 addNode(base.getAlongIndex() - 1, base.getAcrossIndex() - 1, mesh, newNodes);
598 addNode(base.getAlongIndex() - 1, base.getAcrossIndex(), mesh, newNodes);
599 addNode(base.getAlongIndex() - 1, base.getAcrossIndex() + 1, mesh, newNodes);
600 addNode(base.getAlongIndex(), base.getAcrossIndex() - 1, mesh, newNodes);
601 addNode(base.getAlongIndex(), base.getAcrossIndex() + 1, mesh, newNodes);
602 addNode(base.getAlongIndex() + 1, base.getAcrossIndex() - 1, mesh, newNodes);
603 addNode(base.getAlongIndex() + 1, base.getAcrossIndex(), mesh, newNodes);
604 addNode(base.getAlongIndex() + 1, base.getAcrossIndex() + 1, mesh, newNodes);
605 }
606
607 /** Add a node to a mesh if not already present.
608 * @param alongIndex index in the along direction
609 * @param acrossIndex index in the across direction
610 * @param mesh complete mesh containing nodes
611 * @param newNodes queue where new node must be put
612 */
613 private void addNode(final int alongIndex, final int acrossIndex,
614 final Mesh mesh, final Collection<Mesh.Node> newNodes) {
615
616 final Mesh.Node node = mesh.addNode(alongIndex, acrossIndex);
617
618 if (!node.isEnabled()) {
619 // enable the node
620 node.setEnabled();
621 newNodes.add(node);
622 }
623
624 }
625
626 /** Convert a point on the unit 2-sphere to geodetic coordinates.
627 * @param point point on the unit 2-sphere
628 * @return geodetic point (arbitrarily set at altitude 0)
629 */
630 protected GeodeticPoint toGeodetic(final S2Point point) {
631 return new GeodeticPoint(0.5 * FastMath.PI - point.getPhi(), point.getTheta(), 0.0);
632 }
633
634 /** Build a simple zone (connected zone without holes).
635 * <p>
636 * In order to build more complex zones (not connected or with
637 * holes), the user should directly call Hipparchus
638 * {@link SphericalPolygonsSet} constructors and
639 * {@link RegionFactory region factory} if set operations
640 * are needed (union, intersection, difference ...).
641 * </p>
642 * <p>
643 * Take care that the vertices boundary points must be given <em>counterclockwise</em>.
644 * Using the wrong order defines the complementary of the real zone,
645 * and will often result in tessellation failure as the zone is too
646 * wide.
647 * </p>
648 * @param tolerance angular separation below which points are considered
649 * equal (typically 1.0e-10)
650 * @param points vertices of the boundary, in <em>counterclockwise</em>
651 * order, each point being a two-elements arrays with latitude at index 0
652 * and longitude at index 1
653 * @return a zone defined on the unit 2-sphere
654 */
655 public static SphericalPolygonsSet buildSimpleZone(final double tolerance,
656 final double[]... points) {
657 final S2Point[] vertices = new S2Point[points.length];
658 for (int i = 0; i < points.length; ++i) {
659 vertices[i] = new S2Point(points[i][1], 0.5 * FastMath.PI - points[i][0]);
660 }
661 return new SphericalPolygonsSet(tolerance, vertices);
662 }
663
664 /** Build a simple zone (connected zone without holes).
665 * <p>
666 * In order to build more complex zones (not connected or with
667 * holes), the user should directly call Hipparchus
668 * {@link SphericalPolygonsSet} constructors and
669 * {@link RegionFactory region factory} if set operations
670 * are needed (union, intersection, difference ...).
671 * </p>
672 * <p>
673 * Take care that the vertices boundary points must be given <em>counterclockwise</em>.
674 * Using the wrong order defines the complementary of the real zone,
675 * and will often result in tessellation failure as the zone is too
676 * wide.
677 * </p>
678 * @param tolerance angular separation below which points are considered
679 * equal (typically 1.0e-10)
680 * @param points vertices of the boundary, in <em>counterclockwise</em>
681 * order
682 * @return a zone defined on the unit 2-sphere
683 */
684 public static SphericalPolygonsSet buildSimpleZone(final double tolerance,
685 final GeodeticPoint... points) {
686 final S2Point[] vertices = new S2Point[points.length];
687 for (int i = 0; i < points.length; ++i) {
688 vertices[i] = new S2Point(points[i].getLongitude(),
689 0.5 * FastMath.PI - points[i].getLatitude());
690 }
691 return new SphericalPolygonsSet(tolerance, vertices);
692 }
693
694 /** Estimate an approximate motion in the along direction.
695 * @param start node at start of motion
696 * @param end desired point at end of motion
697 * @return approximate motion in the along direction
698 */
699 private double estimateAlongMotion(final Mesh.Node start, final Vector3D end) {
700 return Vector3D.dotProduct(start.getAlong(), end.subtract(start.getV()));
701 }
702
703 /** Estimate an approximate motion in the across direction.
704 * @param start node at start of motion
705 * @param end desired point at end of motion
706 * @return approximate motion in the across direction
707 */
708 private double estimateAcrossMotion(final Mesh.Node start, final Vector3D end) {
709 return Vector3D.dotProduct(start.getAcross(), end.subtract(start.getV()));
710 }
711
712 /** Check if an arc meets the inside of a zone.
713 * @param s1 first point
714 * @param s2 second point
715 * @param zone zone to check arc against
716 * @return true if the arc meets the inside of the zone
717 */
718 private boolean meetInside(final S2Point s1, final S2Point s2,
719 final SphericalPolygonsSet zone) {
720 final Circle circle = new Circle(s1, s2, zone.getTolerance());
721 final double alpha1 = circle.toSubSpace(s1).getAlpha();
722 final double alpha2 = MathUtils.normalizeAngle(circle.toSubSpace(s2).getAlpha(),
723 alpha1 + FastMath.PI);
724 final SubCircle sub = new SubCircle(circle,
725 new ArcsSet(alpha1, alpha2, zone.getTolerance()));
726 return recurseMeetInside(zone.getTree(false), sub);
727
728 }
729
730 /** Check if an arc meets the inside of a zone.
731 * <p>
732 * This method is heavily based on the Characterization class from
733 * Hipparchus library, also distributed under the terms
734 * of the Apache Software License V2.
735 * </p>
736 * @param node spherical zone node
737 * @param sub arc to characterize
738 * @return true if the arc meets the inside of the zone
739 */
740 private boolean recurseMeetInside(final BSPTree<Sphere2D> node, final SubHyperplane<Sphere2D> sub) {
741
742 if (node.getCut() == null) {
743 // we have reached a leaf node
744 if (sub.isEmpty()) {
745 return false;
746 } else {
747 return (Boolean) node.getAttribute();
748 }
749 } else {
750 final Hyperplane<Sphere2D> hyperplane = node.getCut().getHyperplane();
751 final SubHyperplane.SplitSubHyperplane<Sphere2D> split = sub.split(hyperplane);
752 switch (split.getSide()) {
753 case PLUS:
754 return recurseMeetInside(node.getPlus(), sub);
755 case MINUS:
756 return recurseMeetInside(node.getMinus(), sub);
757 case BOTH:
758 if (recurseMeetInside(node.getPlus(), split.getPlus())) {
759 return true;
760 } else {
761 return recurseMeetInside(node.getMinus(), split.getMinus());
762 }
763 default:
764 // this should not happen
765 throw new OrekitInternalError(null);
766 }
767 }
768 }
769
770 /** Get an iterator over mesh nodes indices.
771 * @param minIndex minimum node index
772 * @param maxIndex maximum node index
773 * @param truncateLast true if we can reduce last tile
774 * @return iterator over mesh nodes indices
775 */
776 private Iterable<Range> nodesIndices(final int minIndex, final int maxIndex, final boolean truncateLast) {
777
778 final int first;
779 if (truncateLast) {
780
781 // truncate last tile rather than balance tiles around the zone of interest
782 first = minIndex;
783
784 } else {
785
786 // balance tiles around the zone of interest rather than truncate last tile
787
788 // number of tiles needed to cover the full indices range
789 final int range = maxIndex - minIndex;
790 final int nbTiles = (range + quantization - 1) / quantization;
791
792 // extra nodes that must be added to complete the tiles
793 final int extraNodes = nbTiles * quantization - range;
794
795 // balance the extra nodes before min index and after maxIndex
796 final int extraBefore = (extraNodes + 1) / 2;
797
798 first = minIndex - extraBefore;
799
800 }
801
802 return new Iterable<Range>() {
803
804 /** {@inheritDoc} */
805 @Override
806 public Iterator<Range> iterator() {
807 return new Iterator<Range>() {
808
809 private int nextLower = first;
810
811 /** {@inheritDoc} */
812 @Override
813 public boolean hasNext() {
814 return nextLower < maxIndex;
815 }
816
817 /** {@inheritDoc} */
818 @Override
819 public Range next() {
820
821 if (nextLower >= maxIndex) {
822 throw new NoSuchElementException();
823 }
824 final int lower = nextLower;
825
826 nextLower += quantization;
827 if (truncateLast && nextLower > maxIndex) {
828 // truncate last tile
829 nextLower = maxIndex;
830 }
831
832 return new Range(lower, nextLower);
833
834 }
835
836 /** {@inheritDoc} */
837 @Override
838 public void remove() {
839 throw new UnsupportedOperationException();
840 }
841
842 };
843 }
844 };
845
846 }
847
848 /** Local class for a range of indices to be used for building a tile. */
849 private static class Range {
850
851 /** Lower index. */
852 private final int lower;
853
854 /** Upper index. */
855 private final int upper;
856
857 /** Simple constructor.
858 * @param lower lower index
859 * @param upper upper index
860 */
861 Range(final int lower, final int upper) {
862 this.lower = lower;
863 this.upper = upper;
864 }
865
866 }
867
868 /** Local class for a pair of ranges of indices to be used for building a tile. */
869 private static class RangePair {
870
871 /** Across range. */
872 private final Range across;
873
874 /** Along range. */
875 private final Range along;
876
877 /** Simple constructor.
878 * @param across across range
879 * @param along along range
880 */
881 RangePair(final Range across, final Range along) {
882 this.across = across;
883 this.along = along;
884 }
885
886 }
887
888 }