1 /* Copyright 2002-2020 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.HashMap;
21 import java.util.List;
22 import java.util.Map;
23
24 import org.hipparchus.geometry.euclidean.threed.Vector3D;
25 import org.hipparchus.geometry.euclidean.twod.Vector2D;
26 import org.hipparchus.geometry.partitioning.Region.Location;
27 import org.hipparchus.geometry.spherical.twod.S2Point;
28 import org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet;
29 import org.hipparchus.util.FastMath;
30 import org.hipparchus.util.Precision;
31 import org.orekit.bodies.Ellipse;
32 import org.orekit.bodies.GeodeticPoint;
33 import org.orekit.bodies.OneAxisEllipsoid;
34 import org.orekit.errors.OrekitException;
35 import org.orekit.errors.OrekitMessages;
36
37 /** Class creating a mesh for spherical zones tessellation.
38 * @author Luc Maisonobe
39 */
40 class Mesh {
41
42 /** Underlying ellipsoid. */
43 private final OneAxisEllipsoid ellipsoid;
44
45 /** Zone of interest to tessellate. */
46 private final SphericalPolygonsSet zone;
47
48 /** Zone covered by the mesh. */
49 private SphericalPolygonsSet coverage;
50
51 /** Aiming used for orienting tiles. */
52 private final TileAiming aiming;
53
54 /** Distance between nodes in the along direction. */
55 private final double alongGap;
56
57 /** Distance between nodes in the across direction. */
58 private final double acrossGap;
59
60 /** Map containing nodes. */
61 private final Map<Long, Node> nodes;
62
63 /** Minimum along tile index. */
64 private int minAlongIndex;
65
66 /** Maximum along tile index. */
67 private int maxAlongIndex;
68
69 /** Minimum across tile index. */
70 private int minAcrossIndex;
71
72 /** Maximum across tile index. */
73 private int maxAcrossIndex;
74
75 /** Simple constructor.
76 * @param ellipsoid underlying ellipsoid
77 * @param zone zone of interest to tessellate
78 * @param aiming aiming used for orienting tiles
79 * @param alongGap distance between nodes in the along direction
80 * @param acrossGap distance between nodes in the across direction
81 * @param start location of the first node.
82 */
83 Mesh(final OneAxisEllipsoid ellipsoid, final SphericalPolygonsSet zone,
84 final TileAiming aiming, final double alongGap, final double acrossGap,
85 final S2Point start) {
86 this.ellipsoid = ellipsoid;
87 this.zone = zone;
88 this.coverage = null;
89 this.aiming = aiming;
90 this.alongGap = alongGap;
91 this.acrossGap = acrossGap;
92 this.nodes = new HashMap<Long, Node>();
93 this.minAlongIndex = 0;
94 this.maxAlongIndex = 0;
95 this.minAcrossIndex = 0;
96 this.maxAcrossIndex = 0;
97
98 // safety check for singular points (typically poles)
99 for (final GeodeticPoint singular : aiming.getSingularPoints()) {
100 final S2Point s2p = new S2Point(singular.getLongitude(), 0.5 * FastMath.PI - singular.getLatitude());
101 if (zone.checkPoint(s2p) != Location.OUTSIDE) {
102 throw new OrekitException(OrekitMessages.CANNOT_COMPUTE_AIMING_AT_SINGULAR_POINT,
103 FastMath.toDegrees(singular.getLatitude()),
104 FastMath.toDegrees(singular.getLongitude()));
105 }
106 }
107
108 // create an enabled first node at origin
109 final Node origin = new Node(start, 0, 0);
110 origin.setEnabled();
111
112 // force the first node to be considered inside
113 // It may appear outside if the zone is very thin and
114 // BSPTree.checkPoint selects a very close but wrong
115 // tree leaf tree for the point. Even in this case,
116 // we want the mesh to be properly defined and surround
117 // the area
118 origin.forceInside();
119
120 store(origin);
121
122 }
123
124 /** Get the minimum along tile index for enabled nodes.
125 * @return minimum along tile index for enabled nodes
126 */
127 public int getMinAlongIndex() {
128 return minAlongIndex;
129 }
130
131 /** Get the maximum along tile index for enabled nodes.
132 * @return maximum along tile index for enabled nodes
133 */
134 public int getMaxAlongIndex() {
135 return maxAlongIndex;
136 }
137
138 /** Get the minimum along tile index for enabled nodes for a specific across index.
139 * @param acrossIndex across index to use
140 * @return minimum along tile index for enabled nodes for a specific across index
141 * or {@link #getMaxAlongIndex() getMaxAlongIndex() + 1} if there
142 * are no nodes with the specified acrossIndex.
143 */
144 public int getMinAlongIndex(final int acrossIndex) {
145 for (int alongIndex = minAlongIndex; alongIndex <= maxAlongIndex; ++alongIndex) {
146 final Node node = getNode(alongIndex, acrossIndex);
147 if (node != null && node.isEnabled()) {
148 return alongIndex;
149 }
150 }
151 return maxAlongIndex + 1;
152 }
153
154 /** Get the maximum along tile index for enabled nodes for a specific across index.
155 * @param acrossIndex across index to use
156 * @return maximum along tile index for enabled nodes for a specific across index
157 * or {@link #getMinAlongIndex() getMinAlongIndex() - 1} if there
158 * are no nodes with the specified acrossIndex.
159 */
160 public int getMaxAlongIndex(final int acrossIndex) {
161 for (int alongIndex = maxAlongIndex; alongIndex >= minAlongIndex; --alongIndex) {
162 final Node node = getNode(alongIndex, acrossIndex);
163 if (node != null && node.isEnabled()) {
164 return alongIndex;
165 }
166 }
167 return minAlongIndex - 1;
168 }
169
170 /** Get the minimum across tile index.
171 * @return minimum across tile index
172 */
173 public int getMinAcrossIndex() {
174 return minAcrossIndex;
175 }
176
177 /** Get the maximum across tile index.
178 * @return maximum across tile index
179 */
180 public int getMaxAcrossIndex() {
181 return maxAcrossIndex;
182 }
183
184 /** Get the number of nodes.
185 * @return number of nodes
186 */
187 public int getNumberOfNodes() {
188 return nodes.size();
189 }
190
191 /** Get the distance between nodes in the along direction.
192 * @return distance between nodes in the along direction
193 */
194 public double getAlongGap() {
195 return alongGap;
196 }
197
198 /** Get the distance between nodes in the across direction.
199 * @return distance between nodes in the across direction
200 */
201 public double getAcrossGap() {
202 return acrossGap;
203 }
204
205 /** Retrieve a node from its indices.
206 * @param alongIndex index in the along direction
207 * @param acrossIndex index in the across direction
208 * @return node or null if no nodes are available at these indices
209 * @see #addNode(int, int)
210 */
211 public Node getNode(final int alongIndex, final int acrossIndex) {
212 return nodes.get(key(alongIndex, acrossIndex));
213 }
214
215 /** Add a node.
216 * <p>
217 * This method is similar to {@link #getNode(int, int) getNode} except
218 * it creates the node if it doesn't alreay exists. All created nodes
219 * have a status set to {@code disabled}.
220 * </p>
221 * @param alongIndex index in the along direction
222 * @param acrossIndex index in the across direction
223 * @return node at specified indices, guaranteed to be non-null
224 * @see #getNode(int, int)
225 */
226 public Node addNode(final int alongIndex, final int acrossIndex) {
227
228 // create intermediate (disabled) nodes, up to specified indices
229 Node node = getExistingAncestor(alongIndex, acrossIndex);
230 while (node.getAlongIndex() != alongIndex || node.getAcrossIndex() != acrossIndex) {
231 final Direction direction;
232 if (node.getAlongIndex() < alongIndex) {
233 direction = Direction.PLUS_ALONG;
234 } else if (node.getAlongIndex() > alongIndex) {
235 direction = Direction.MINUS_ALONG;
236 } else if (node.getAcrossIndex() < acrossIndex) {
237 direction = Direction.PLUS_ACROSS;
238 } else {
239 direction = Direction.MINUS_ACROSS;
240 }
241 final S2Point s2p = node.move(direction.motion(node, alongGap, acrossGap));
242 node = new Node(s2p, direction.neighborAlongIndex(node), direction.neighborAcrossIndex(node));
243 store(node);
244 }
245
246 return node;
247
248 }
249
250 /** Find the closest existing ancestor of a node.
251 * <p>
252 * The path from origin to any node is first in the along direction,
253 * and then in the across direction. Using always the same path pattern
254 * ensures consistent distortion of the mesh.
255 * </p>
256 * @param alongIndex index in the along direction
257 * @param acrossIndex index in the across direction
258 * @return an existing node in the path from origin to specified indices
259 */
260 private Node getExistingAncestor(final int alongIndex, final int acrossIndex) {
261
262 // start from the desired node indices
263 int l = alongIndex;
264 int c = acrossIndex;
265 Node node = getNode(l, c);
266
267 // rewind the path backward, up to origin,
268 // that is first in the across direction, then in the along direction
269 // the loop WILL end as there is at least one node at (0, 0)
270 while (node == null) {
271 if (c != 0) {
272 c += c > 0 ? -1 : +1;
273 } else {
274 l += l > 0 ? -1 : +1;
275 }
276 node = getNode(l, c);
277 }
278
279 // we have found an existing ancestor
280 return node;
281
282 }
283
284 /** Get the nodes that lie inside the interest zone.
285 * @return nodes that lie inside the interest zone
286 */
287 public List<Node> getInsideNodes() {
288 final List<Node> insideNodes = new ArrayList<Node>();
289 for (final Map.Entry<Long, Node> entry : nodes.entrySet()) {
290 if (entry.getValue().isInside()) {
291 insideNodes.add(entry.getValue());
292 }
293 }
294 return insideNodes;
295 }
296
297 /** Find the existing node closest to a location.
298 * @param alongIndex index in the along direction
299 * @param acrossIndex index in the across direction
300 * @return node or null if no node is available at these indices
301 */
302 public Node getClosestExistingNode(final int alongIndex, final int acrossIndex) {
303
304 // we know at least the first node (at 0, 0) exists, we use its
305 // taxicab distance to the desired location as a search limit
306 final int maxD = FastMath.max(FastMath.abs(alongIndex), FastMath.abs(acrossIndex));
307
308 // search for an existing point in increasing taxicab distances
309 for (int d = 0; d < maxD; ++d) {
310 for (int deltaAcross = 0; deltaAcross <= d; ++deltaAcross) {
311 final int deltaAlong = d - deltaAcross;
312 Node node = getNode(alongIndex - deltaAlong, acrossIndex - deltaAcross);
313 if (node != null) {
314 return node;
315 }
316 if (deltaAcross != 0) {
317 node = getNode(alongIndex - deltaAlong, acrossIndex + deltaAcross);
318 if (node != null) {
319 return node;
320 }
321 }
322 if (deltaAlong != 0) {
323 node = getNode(alongIndex + deltaAlong, acrossIndex - deltaAcross);
324 if (node != null) {
325 return node;
326 }
327 if (deltaAcross != 0) {
328 node = getNode(alongIndex + deltaAlong, acrossIndex + deltaAcross);
329 if (node != null) {
330 return node;
331 }
332 }
333 }
334 }
335 }
336
337 // at least the first node always exists
338 return getNode(0, 0);
339
340 }
341
342 /** Find the existing node closest to a location.
343 * @param location reference location in Cartesian coordinates
344 * @return node or null if no node is available at these indices
345 */
346 public Node getClosestExistingNode(final Vector3D location) {
347 Node selected = null;
348 double min = Double.POSITIVE_INFINITY;
349 for (final Map.Entry<Long, Node> entry : nodes.entrySet()) {
350 final double distance = Vector3D.distance(location, entry.getValue().getV());
351 if (distance < min) {
352 selected = entry.getValue();
353 min = distance;
354 }
355 }
356 return selected;
357 }
358
359 /** Get the oriented list of <em>enabled</em> nodes at mesh boundary, in taxicab geometry.
360 * @param simplified if true, don't include intermediate points along almost straight lines
361 * @return list of nodes
362 */
363 public List<Node> getTaxicabBoundary(final boolean simplified) {
364
365 final List<Node> boundary = new ArrayList<Node>();
366 if (nodes.size() < 2) {
367 boundary.add(getNode(0, 0));
368 } else {
369
370 // search for the lower left corner
371 Node lowerLeft = null;
372 for (int i = minAlongIndex; lowerLeft == null && i <= maxAlongIndex; ++i) {
373 for (int j = minAcrossIndex; lowerLeft == null && j <= maxAcrossIndex; ++j) {
374 lowerLeft = getNode(i, j);
375 if (lowerLeft != null && !lowerLeft.isEnabled()) {
376 lowerLeft = null;
377 }
378 }
379 }
380
381 // loop counterclockwise around the mesh
382 Direction direction = Direction.MINUS_ACROSS;
383 Node node = lowerLeft;
384 do {
385 boundary.add(node);
386 Node neighbor = null;
387 do {
388 direction = direction.next();
389 neighbor = getNode(direction.neighborAlongIndex(node),
390 direction.neighborAcrossIndex(node));
391 } while (neighbor == null || !neighbor.isEnabled());
392 direction = direction.next().next();
393 node = neighbor;
394 } while (node != lowerLeft);
395 }
396
397 // filter out infinitely thin parts corresponding to spikes
398 // joining outliers points to the main mesh
399 boolean changed = true;
400 while (changed && boundary.size() > 1) {
401 changed = false;
402 final int n = boundary.size();
403 for (int i = 0; i < n; ++i) {
404 final int previousIndex = (i + n - 1) % n;
405 final int nextIndex = (i + 1) % n;
406 if (boundary.get(previousIndex) == boundary.get(nextIndex)) {
407 // the current point is an infinitely thin spike, remove it
408 boundary.remove(FastMath.max(i, nextIndex));
409 boundary.remove(FastMath.min(i, nextIndex));
410 changed = true;
411 break;
412 }
413 }
414 }
415
416 if (simplified) {
417 for (int i = 0; i < boundary.size(); ++i) {
418 final int n = boundary.size();
419 final Node previous = boundary.get((i + n - 1) % n);
420 final int pl = previous.getAlongIndex();
421 final int pc = previous.getAcrossIndex();
422 final Node current = boundary.get(i);
423 final int cl = current.getAlongIndex();
424 final int cc = current.getAcrossIndex();
425 final Node next = boundary.get((i + 1) % n);
426 final int nl = next.getAlongIndex();
427 final int nc = next.getAcrossIndex();
428 if ((pl == cl && cl == nl) || (pc == cc && cc == nc)) {
429 // the current point is a spurious intermediate in a straight line, remove it
430 boundary.remove(i--);
431 }
432 }
433 }
434
435 return boundary;
436
437 }
438
439 /** Get the zone covered by the mesh.
440 * @return mesh coverage
441 */
442 public SphericalPolygonsSet getCoverage() {
443
444 if (coverage == null) {
445
446 // lazy build of mesh coverage
447 final List<Mesh.Node> boundary = getTaxicabBoundary(true);
448 final S2Point[] vertices = new S2Point[boundary.size()];
449 for (int i = 0; i < vertices.length; ++i) {
450 vertices[i] = boundary.get(i).getS2P();
451 }
452 coverage = new SphericalPolygonsSet(zone.getTolerance(), vertices);
453 }
454
455 // as caller may modify the BSP tree, we must provide a copy of our safe instance
456 return (SphericalPolygonsSet) coverage.copySelf();
457
458 }
459
460 /** Store a node.
461 * @param node to add
462 */
463 private void store(final Node node) {
464
465 // the new node invalidates current estimation of the coverage
466 coverage = null;
467
468 nodes.put(key(node.alongIndex, node.acrossIndex), node);
469
470 }
471
472 /** Convert along and across indices to map key.
473 * @param alongIndex index in the along direction
474 * @param acrossIndex index in the across direction
475 * @return key map key
476 */
477 private long key(final int alongIndex, final int acrossIndex) {
478 return ((long) alongIndex) << 32 | (((long) acrossIndex) & 0xFFFFl);
479 }
480
481 /** Container for mesh nodes. */
482 public class Node {
483
484 /** Node position in spherical coordinates. */
485 private final S2Point s2p;
486
487 /** Node position in Cartesian coordinates. */
488 private final Vector3D v;
489
490 /** Normalized along tile direction. */
491 private final Vector3D along;
492
493 /** Normalized across tile direction. */
494 private final Vector3D across;
495
496 /** Indicator for node location with respect to interest zone. */
497 private boolean insideZone;
498
499 /** Index in the along direction. */
500 private final int alongIndex;
501
502 /** Index in the across direction. */
503 private final int acrossIndex;
504
505 /** Indicator for construction nodes used only as intermediate points (disabled) vs. real nodes (enabled). */
506 private boolean enabled;
507
508 /** Create a node.
509 * @param s2p position in spherical coordinates (my be null)
510 * @param alongIndex index in the along direction
511 * @param acrossIndex index in the across direction
512 */
513 private Node(final S2Point s2p, final int alongIndex, final int acrossIndex) {
514 final GeodeticPointint.html#GeodeticPoint">GeodeticPoint gp = new GeodeticPoint(0.5 * FastMath.PI - s2p.getPhi(), s2p.getTheta(), 0.0);
515 this.v = ellipsoid.transform(gp);
516 this.s2p = s2p;
517 this.along = aiming.alongTileDirection(v, gp);
518 this.across = Vector3D.crossProduct(v, along).normalize();
519 this.insideZone = zone.checkPoint(s2p) != Location.OUTSIDE;
520 this.alongIndex = alongIndex;
521 this.acrossIndex = acrossIndex;
522 this.enabled = false;
523 }
524
525 /** Set the enabled property.
526 */
527 public void setEnabled() {
528
529 // store status
530 this.enabled = true;
531
532 // update min/max indices
533 minAlongIndex = FastMath.min(minAlongIndex, alongIndex);
534 maxAlongIndex = FastMath.max(maxAlongIndex, alongIndex);
535 minAcrossIndex = FastMath.min(minAcrossIndex, acrossIndex);
536 maxAcrossIndex = FastMath.max(maxAcrossIndex, acrossIndex);
537
538 }
539
540 /** Check if a node is enabled.
541 * @return true if the node is enabled
542 */
543 public boolean isEnabled() {
544 return enabled;
545 }
546
547 /** Get the node position in spherical coordinates.
548 * @return vode position in spherical coordinates
549 */
550 public S2Point getS2P() {
551 return s2p;
552 }
553
554 /** Get the node position in Cartesian coordinates.
555 * @return vode position in Cartesian coordinates
556 */
557 public Vector3D getV() {
558 return v;
559 }
560
561 /** Get the normalized along tile direction.
562 * @return normalized along tile direction
563 */
564 public Vector3D getAlong() {
565 return along;
566 }
567
568 /** Get the normalized across tile direction.
569 * @return normalized across tile direction
570 */
571 public Vector3D getAcross() {
572 return across;
573 }
574
575 /** Force the node location to be considered inside interest zone.
576 */
577 private void forceInside() {
578 insideZone = true;
579 }
580
581 /** Check if the node location is inside interest zone.
582 * @return true if the node location is inside interest zone
583 */
584 public boolean isInside() {
585 return insideZone;
586 }
587
588 /** Get the index in the along direction.
589 * @return index in the along direction
590 */
591 public int getAlongIndex() {
592 return alongIndex;
593 }
594
595 /** Get the index in the across direction.
596 * @return index in the across direction
597 */
598 public int getAcrossIndex() {
599 return acrossIndex;
600 }
601
602 /** Move to a nearby point along surface.
603 * <p>
604 * The motion will be approximated, assuming the body surface has constant
605 * curvature along the motion direction. The approximation will be accurate
606 * for short distances, and error will increase as distance increases.
607 * </p>
608 * @param motion straight motion, which must be curved back to surface
609 * @return arrival point, approximately at specified distance from node
610 */
611 public S2Point move(final Vector3D motion) {
612
613 // safety check for too small motion
614 if (motion.getNorm() < Precision.EPSILON * v.getNorm()) {
615 return s2p;
616 }
617
618 // find elliptic plane section
619 final Vector3D normal = Vector3D.crossProduct(v, motion);
620 final Ellipse planeSection = ellipsoid.getPlaneSection(v, normal);
621
622 // find the center of curvature (point on the evolute) below start point
623 final Vector2D omega2D = planeSection.getCenterOfCurvature(planeSection.toPlane(v));
624 final Vector3D omega3D = planeSection.toSpace(omega2D);
625
626 // compute approximated arrival point, assuming constant radius of curvature
627 final Vector3D delta = v.subtract(omega3D);
628 final double theta = motion.getNorm() / delta.getNorm();
629 final Vector3D approximated = new Vector3D(1, omega3D,
630 FastMath.cos(theta), delta,
631 FastMath.sin(theta) / theta, motion);
632
633 // convert to spherical coordinates
634 final GeodeticPoint approximatedGP = ellipsoid.transform(approximated, ellipsoid.getBodyFrame(), null);
635 return new S2Point(approximatedGP.getLongitude(), 0.5 * FastMath.PI - approximatedGP.getLatitude());
636
637 }
638
639 }
640
641 }