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