EllipsoidTessellator.java

  1. /* Copyright 2002-2022 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. import java.util.ArrayList;
  19. import java.util.Collection;
  20. import java.util.IdentityHashMap;
  21. import java.util.Iterator;
  22. import java.util.LinkedList;
  23. import java.util.List;
  24. import java.util.Map;
  25. import java.util.NoSuchElementException;
  26. import java.util.Queue;

  27. import org.hipparchus.exception.LocalizedCoreFormats;
  28. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  29. import org.hipparchus.geometry.partitioning.BSPTree;
  30. import org.hipparchus.geometry.partitioning.Hyperplane;
  31. import org.hipparchus.geometry.partitioning.RegionFactory;
  32. import org.hipparchus.geometry.partitioning.SubHyperplane;
  33. import org.hipparchus.geometry.spherical.oned.ArcsSet;
  34. import org.hipparchus.geometry.spherical.twod.Circle;
  35. import org.hipparchus.geometry.spherical.twod.S2Point;
  36. import org.hipparchus.geometry.spherical.twod.Sphere2D;
  37. import org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet;
  38. import org.hipparchus.geometry.spherical.twod.SubCircle;
  39. import org.hipparchus.util.FastMath;
  40. import org.hipparchus.util.MathUtils;
  41. import org.orekit.bodies.GeodeticPoint;
  42. import org.orekit.bodies.OneAxisEllipsoid;
  43. import org.orekit.errors.OrekitException;
  44. import org.orekit.errors.OrekitInternalError;

  45. /** Class used to tessellate an interest zone on an ellipsoid in either
  46.  * {@link Tile tiles} or grids of {@link GeodeticPoint geodetic points}.
  47.  * <p>
  48.  * This class is typically used for Earth Observation missions, in order to
  49.  * create tiles or grids that may be used as the basis of visibility event
  50.  * detectors. Tiles are used when surface-related elements are needed, the
  51.  * tiles created completely cover the zone of interest. Grids are used when
  52.  * point-related elements are needed, the points created lie entirely within
  53.  * the zone of interest.
  54.  * </p>
  55.  * <p>
  56.  * One should note that as tessellation essentially creates a 2 dimensional
  57.  * almost Cartesian map, it can never perfectly fulfill geometrical dimensions
  58.  * because neither sphere nor ellipsoid are developable surfaces. This implies
  59.  * that the tesselation will always be distorted, and distortion increases as
  60.  * the size of the zone to be tessellated increases.
  61.  * </p>
  62.  * @author Luc Maisonobe
  63.  * @since 7.1
  64.  */
  65. public class EllipsoidTessellator {

  66.     /** Safety limit to avoid infinite loops during tesselation due to numerical noise.
  67.      * @since 10.3.1
  68.      */
  69.     private static final int MAX_ITER = 1000;

  70.     /** Number of segments tiles sides are split into for tiles fine positioning. */
  71.     private final int quantization;

  72.     /** Aiming used for orienting tiles. */
  73.     private final TileAiming aiming;

  74.     /** Underlying ellipsoid. */
  75.     private final OneAxisEllipsoid ellipsoid;

  76.     /** Simple constructor.
  77.      * <p>
  78.      * The {@code quantization} parameter is used internally to adjust points positioning.
  79.      * For example when quantization is set to 4, a complete tile that has 4 corner points
  80.      * separated by the tile lengths will really be computed on a grid containing 25 points
  81.      * (5 rows of 5 points, as each side will be split in 4 segments, hence will have 5
  82.      * points). This quantization allows rough adjustment to balance margins around the
  83.      * zone of interest and improves geometric accuracy as the along and across directions
  84.      * are readjusted at each points.
  85.      * </p>
  86.      * <p>
  87.      * It is recommended to use at least 2 as the quantization parameter for tiling. The
  88.      * rationale is that using only 1 for quantization would imply all points used are tiles
  89.      * vertices, and hence would lead small zones to generate 4 tiles with a shared vertex
  90.      * inside the zone and the 4 tiles covering the four quadrants at North-West, North-East,
  91.      * South-East and South-West. A quantization value of at least 2 allows to shift the
  92.      * tiles so the center point is an inside point rather than a tile vertex, hence allowing
  93.      * a single tile to cover the small zone. A value even greater like 4 or 8 would allow even
  94.      * finer positioning to balance the tiles margins around the zone.
  95.      * </p>
  96.      * @param ellipsoid underlying ellipsoid
  97.      * @param aiming aiming used for orienting tiles
  98.      * @param quantization number of segments tiles sides are split into for tiles fine positioning
  99.      */
  100.     public EllipsoidTessellator(final OneAxisEllipsoid ellipsoid, final TileAiming aiming,
  101.                                 final int quantization) {
  102.         this.ellipsoid    = ellipsoid;
  103.         this.aiming       = aiming;
  104.         this.quantization = quantization;
  105.     }

  106.     /** Tessellate a zone of interest into tiles.
  107.      * <p>
  108.      * The created tiles will completely cover the zone of interest.
  109.      * </p>
  110.      * <p>
  111.      * The distance between a vertex at a tile corner and the vertex at the same corner
  112.      * in the next vertex are computed by subtracting the overlap width (resp. overlap length)
  113.      * from the full width (resp. full length). If for example the full width is specified to
  114.      * be 55 km and the overlap in width is specified to be +5 km, successive tiles would span
  115.      * as follows:
  116.      * </p>
  117.      * <ul>
  118.      *   <li>tile 1 covering from   0 km to  55 km</li>
  119.      *   <li>tile 2 covering from  50 km to 105 km</li>
  120.      *   <li>tile 3 covering from 100 km to 155 km</li>
  121.      *   <li>...</li>
  122.      * </ul>
  123.      * <p>
  124.      * In order to achieve the same 50 km step but using a 5 km gap instead of an overlap, one would
  125.      * need to specify the full width to be 45 km and the overlap to be -5 km. With these settings,
  126.      * successive tiles would span as follows:
  127.      * </p>
  128.      * <ul>
  129.      *   <li>tile 1 covering from   0 km to  45 km</li>
  130.      *   <li>tile 2 covering from  50 km to  95 km</li>
  131.      *   <li>tile 3 covering from 100 km to 155 km</li>
  132.      *   <li>...</li>
  133.      * </ul>
  134.      * @param zone zone of interest to tessellate
  135.      * @param fullWidth full tiles width as a distance on surface, including overlap (in meters)
  136.      * @param fullLength full tiles length as a distance on surface, including overlap (in meters)
  137.      * @param widthOverlap overlap between adjacent tiles (in meters), if negative the tiles
  138.      * will have a gap between each other instead of an overlap
  139.      * @param lengthOverlap overlap between adjacent tiles (in meters), if negative the tiles
  140.      * will have a gap between each other instead of an overlap
  141.      * @param truncateLastWidth if true, the first tiles strip will be started as close as
  142.      * possible to the zone of interest, and the last tiles strip will have its width reduced
  143.      * to also remain close to the zone of interest; if false all tiles strip will have the
  144.      * same {@code fullWidth} and they will be balanced around zone of interest
  145.      * @param truncateLastLength if true, the first tile in each strip will be started as close as
  146.      * possible to the zone of interest, and the last tile in each strip will have its length reduced
  147.      * to also remain close to the zone of interest; if false all tiles in each strip will have the
  148.      * same {@code fullLength} and they will be balanced around zone of interest
  149.      * @return a list of lists of tiles covering the zone of interest,
  150.      * each sub-list corresponding to a part not connected to the other
  151.      * parts (for example for islands)
  152.      */
  153.     public List<List<Tile>> tessellate(final SphericalPolygonsSet zone,
  154.                                        final double fullWidth, final double fullLength,
  155.                                        final double widthOverlap, final double lengthOverlap,
  156.                                        final boolean truncateLastWidth, final boolean truncateLastLength) {

  157.         final double                  splitWidth  = (fullWidth  - widthOverlap)  / quantization;
  158.         final double                  splitLength = (fullLength - lengthOverlap) / quantization;
  159.         final Map<Mesh, List<Tile>>   map         = new IdentityHashMap<Mesh, List<Tile>>();
  160.         final RegionFactory<Sphere2D> factory     = new RegionFactory<Sphere2D>();
  161.         SphericalPolygonsSet          remaining   = (SphericalPolygonsSet) zone.copySelf();
  162.         S2Point                       inside      = getInsidePoint(remaining);

  163.         int count = 0;
  164.         while (inside != null) {

  165.             if (++count > MAX_ITER) {
  166.                 throw new OrekitException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, MAX_ITER);
  167.             }

  168.             // find a mesh covering at least one connected part of the zone
  169.             final List<Mesh.Node> mergingSeeds = new ArrayList<Mesh.Node>();
  170.             Mesh mesh = new Mesh(ellipsoid, zone, aiming, splitLength, splitWidth, inside);
  171.             mergingSeeds.add(mesh.getNode(0, 0));
  172.             List<Tile> tiles = null;
  173.             while (!mergingSeeds.isEmpty()) {

  174.                 // expand the mesh around the seed
  175.                 neighborExpandMesh(mesh, mergingSeeds, zone);

  176.                 // extract the tiles from the mesh
  177.                 // this further expands the mesh so tiles dimensions are multiples of quantization,
  178.                 // hence it must be performed here before checking meshes independence
  179.                 tiles = extractTiles(mesh, zone, lengthOverlap, widthOverlap, truncateLastWidth, truncateLastLength);

  180.                 // check the mesh is independent from existing meshes
  181.                 mergingSeeds.clear();
  182.                 for (final Map.Entry<Mesh, List<Tile>> entry : map.entrySet()) {
  183.                     if (!factory.intersection(mesh.getCoverage(), entry.getKey().getCoverage()).isEmpty()) {
  184.                         // the meshes are not independent, they intersect each other!

  185.                         // merge the two meshes together
  186.                         mesh = mergeMeshes(mesh, entry.getKey(), mergingSeeds);
  187.                         map.remove(entry.getKey());
  188.                         break;

  189.                     }
  190.                 }

  191.             }

  192.             // remove the part of the zone covered by the mesh
  193.             remaining = (SphericalPolygonsSet) factory.difference(remaining, mesh.getCoverage());
  194.             inside    = getInsidePoint(remaining);

  195.             map.put(mesh, tiles);

  196.         }

  197.         // concatenate the lists from the independent meshes
  198.         final List<List<Tile>> tilesLists = new ArrayList<List<Tile>>(map.size());
  199.         for (final Map.Entry<Mesh, List<Tile>> entry : map.entrySet()) {
  200.             tilesLists.add(entry.getValue());
  201.         }

  202.         return tilesLists;

  203.     }

  204.     /** Sample a zone of interest into a grid sample of {@link GeodeticPoint geodetic points}.
  205.      * <p>
  206.      * The created points will be entirely within the zone of interest.
  207.      * </p>
  208.      * @param zone zone of interest to sample
  209.      * @param width grid sample cells width as a distance on surface (in meters)
  210.      * @param length grid sample cells length as a distance on surface (in meters)
  211.      * @return a list of lists of points sampling the zone of interest,
  212.      * each sub-list corresponding to a part not connected to the other
  213.      * parts (for example for islands)
  214.      */
  215.     public List<List<GeodeticPoint>> sample(final SphericalPolygonsSet zone,
  216.                                             final double width, final double length) {

  217.         final double                         splitWidth  = width  / quantization;
  218.         final double                         splitLength = length / quantization;
  219.         final Map<Mesh, List<GeodeticPoint>> map         = new IdentityHashMap<Mesh, List<GeodeticPoint>>();
  220.         final RegionFactory<Sphere2D>        factory     = new RegionFactory<Sphere2D>();
  221.         SphericalPolygonsSet                 remaining   = (SphericalPolygonsSet) zone.copySelf();
  222.         S2Point                              inside      = getInsidePoint(remaining);

  223.         int count = 0;
  224.         while (inside != null) {

  225.             if (++count > MAX_ITER) {
  226.                 throw new OrekitException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, MAX_ITER);
  227.             }

  228.             // find a mesh covering at least one connected part of the zone
  229.             final List<Mesh.Node> mergingSeeds = new ArrayList<Mesh.Node>();
  230.             Mesh mesh = new Mesh(ellipsoid, zone, aiming, splitLength, splitWidth, inside);
  231.             mergingSeeds.add(mesh.getNode(0, 0));
  232.             List<GeodeticPoint> sample = null;
  233.             while (!mergingSeeds.isEmpty()) {

  234.                 // expand the mesh around the seed
  235.                 neighborExpandMesh(mesh, mergingSeeds, zone);

  236.                 // extract the sample from the mesh
  237.                 // this further expands the mesh so sample cells dimensions are multiples of quantization,
  238.                 // hence it must be performed here before checking meshes independence
  239.                 sample = extractSample(mesh, zone);

  240.                 // check the mesh is independent from existing meshes
  241.                 mergingSeeds.clear();
  242.                 for (final Map.Entry<Mesh, List<GeodeticPoint>> entry : map.entrySet()) {
  243.                     if (!factory.intersection(mesh.getCoverage(), entry.getKey().getCoverage()).isEmpty()) {
  244.                         // the meshes are not independent, they intersect each other!

  245.                         // merge the two meshes together
  246.                         mesh = mergeMeshes(mesh, entry.getKey(), mergingSeeds);
  247.                         map.remove(entry.getKey());
  248.                         break;

  249.                     }
  250.                 }

  251.             }

  252.             // remove the part of the zone covered by the mesh
  253.             remaining = (SphericalPolygonsSet) factory.difference(remaining, mesh.getCoverage());
  254.             inside    = getInsidePoint(remaining);

  255.             map.put(mesh, sample);

  256.         }

  257.         // concatenate the lists from the independent meshes
  258.         final List<List<GeodeticPoint>> sampleLists = new ArrayList<List<GeodeticPoint>>(map.size());
  259.         for (final Map.Entry<Mesh, List<GeodeticPoint>> entry : map.entrySet()) {
  260.             sampleLists.add(entry.getValue());
  261.         }

  262.         return sampleLists;

  263.     }

  264.     /** Get an inside point from a zone of interest.
  265.      * @param zone zone to mesh
  266.      * @return a point inside the zone or null if zone is empty or too thin
  267.      */
  268.     private S2Point getInsidePoint(final SphericalPolygonsSet zone) {

  269.         final InsideFinder finder = new InsideFinder(zone);
  270.         zone.getTree(false).visit(finder);
  271.         return finder.getInsidePoint();

  272.     }

  273.     /** Expand a mesh so it surrounds at least one connected part of a zone.
  274.      * <p>
  275.      * This part of mesh expansion is neighbors based. It includes the seed
  276.      * node neighbors, and their neighbors, and the neighbors of their
  277.      * neighbors until the path-connected sub-parts of the zone these nodes
  278.      * belong to are completely surrounded by the mesh taxicab boundary.
  279.      * </p>
  280.      * @param mesh mesh to expand
  281.      * @param seeds seed nodes (already in the mesh) from which to start expansion
  282.      * @param zone zone to mesh
  283.      */
  284.     private void neighborExpandMesh(final Mesh mesh, final Collection<Mesh.Node> seeds,
  285.                                     final SphericalPolygonsSet zone) {

  286.         // mesh expansion loop
  287.         boolean expanding = true;
  288.         final Queue<Mesh.Node> newNodes = new LinkedList<Mesh.Node>();
  289.         newNodes.addAll(seeds);
  290.         int count = 0;
  291.         while (expanding) {

  292.             if (++count > MAX_ITER) {
  293.                 throw new OrekitException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, MAX_ITER);
  294.             }

  295.             // first expansion step: set up the mesh so that all its
  296.             // inside nodes are completely surrounded by at least
  297.             // one layer of outside nodes
  298.             while (!newNodes.isEmpty()) {

  299.                 // retrieve an active node
  300.                 final Mesh.Node node = newNodes.remove();

  301.                 if (node.isInside()) {
  302.                     // the node is inside the zone, the mesh must contain its 8 neighbors
  303.                     addAllNeighborsIfNeeded(node, mesh, newNodes);
  304.                 }

  305.             }

  306.             // second expansion step: check if the loop of outside nodes
  307.             // completely surrounds the zone, i.e. there are no peaks
  308.             // pointing out of the loop between two nodes
  309.             expanding = false;
  310.             final List<Mesh.Node> boundary = mesh.getTaxicabBoundary(false);
  311.             if (boundary.size() > 1) {
  312.                 Mesh.Node previous = boundary.get(boundary.size() - 1);
  313.                 for (final Mesh.Node node : boundary) {
  314.                     if (meetInside(previous.getS2P(), node.getS2P(), zone)) {
  315.                         // part of the mesh boundary is still inside the zone!
  316.                         // the mesh must be expanded again
  317.                         addAllNeighborsIfNeeded(previous, mesh, newNodes);
  318.                         addAllNeighborsIfNeeded(node,     mesh, newNodes);
  319.                         expanding = true;
  320.                     }
  321.                     previous = node;
  322.                 }
  323.             }

  324.         }

  325.     }

  326.     /** Extract tiles from a mesh.
  327.      * @param mesh mesh from which tiles should be extracted
  328.      * @param zone zone covered by the mesh
  329.      * @param lengthOverlap overlap between adjacent tiles
  330.      * @param widthOverlap overlap between adjacent tiles
  331.      * @param truncateLastWidth true if we can reduce last tile width
  332.      * @param truncateLastLength true if we can reduce last tile length
  333.      * @return extracted tiles
  334.      */
  335.     private List<Tile> extractTiles(final Mesh mesh, final SphericalPolygonsSet zone,
  336.                                     final double lengthOverlap, final double widthOverlap,
  337.                                     final boolean truncateLastWidth, final boolean truncateLastLength) {

  338.         final List<Tile>      tiles = new ArrayList<Tile>();
  339.         final List<RangePair> rangePairs = new ArrayList<RangePair>();

  340.         final int minAcross = mesh.getMinAcrossIndex();
  341.         final int maxAcross = mesh.getMaxAcrossIndex();
  342.         for (Range acrossPair : nodesIndices(minAcross, maxAcross, truncateLastWidth)) {

  343.             int minAlong = mesh.getMaxAlongIndex() + 1;
  344.             int maxAlong = mesh.getMinAlongIndex() - 1;
  345.             for (int c = acrossPair.lower; c <= acrossPair.upper; ++c) {
  346.                 minAlong = FastMath.min(minAlong, mesh.getMinAlongIndex(c));
  347.                 maxAlong = FastMath.max(maxAlong, mesh.getMaxAlongIndex(c));
  348.             }

  349.             for (Range alongPair : nodesIndices(minAlong, maxAlong, truncateLastLength)) {

  350.                 // get the base vertex nodes
  351.                 final Mesh.Node node0 = mesh.addNode(alongPair.lower, acrossPair.lower);
  352.                 final Mesh.Node node1 = mesh.addNode(alongPair.upper, acrossPair.lower);
  353.                 final Mesh.Node node2 = mesh.addNode(alongPair.upper, acrossPair.upper);
  354.                 final Mesh.Node node3 = mesh.addNode(alongPair.lower, acrossPair.upper);

  355.                 // apply tile overlap
  356.                 final S2Point s2p0 = node0.move(new Vector3D(-0.5 * lengthOverlap, node0.getAlong(),
  357.                                                              -0.5 * widthOverlap,  node0.getAcross()));
  358.                 final S2Point s2p1 = node1.move(new Vector3D(+0.5 * lengthOverlap, node1.getAlong(),
  359.                                                              -0.5 * widthOverlap,  node1.getAcross()));
  360.                 final S2Point s2p2 = node2.move(new Vector3D(+0.5 * lengthOverlap, node2.getAlong(),
  361.                                                              +0.5 * widthOverlap,  node2.getAcross()));
  362.                 final S2Point s2p3 = node3.move(new Vector3D(-0.5 * lengthOverlap, node2.getAlong(),
  363.                                                              +0.5 * widthOverlap,  node2.getAcross()));

  364.                 // create a quadrilateral region corresponding to the candidate tile
  365.                 final SphericalPolygonsSet quadrilateral =
  366.                         new SphericalPolygonsSet(zone.getTolerance(), s2p0, s2p1, s2p2, s2p3);

  367.                 if (!new RegionFactory<Sphere2D>().intersection(zone.copySelf(), quadrilateral).isEmpty()) {

  368.                     // the tile does cover part of the zone, it contributes to the tessellation
  369.                     tiles.add(new Tile(toGeodetic(s2p0), toGeodetic(s2p1), toGeodetic(s2p2), toGeodetic(s2p3)));
  370.                     rangePairs.add(new RangePair(acrossPair, alongPair));

  371.                 }

  372.             }
  373.         }

  374.         // ensure the taxicab boundary follows the built tile sides
  375.         // this is done outside of the previous loop in order
  376.         // to avoid one tile changing the min/max indices of the
  377.         // neighboring tile as they share some nodes that will be enabled here
  378.         for (final RangePair rangePair : rangePairs) {
  379.             for (int c = rangePair.across.lower; c < rangePair.across.upper; ++c) {
  380.                 mesh.addNode(rangePair.along.lower, c + 1).setEnabled();
  381.                 mesh.addNode(rangePair.along.upper, c).setEnabled();
  382.             }
  383.             for (int l = rangePair.along.lower; l < rangePair.along.upper; ++l) {
  384.                 mesh.addNode(l,     rangePair.across.lower).setEnabled();
  385.                 mesh.addNode(l + 1, rangePair.across.upper).setEnabled();
  386.             }
  387.         }

  388.         return tiles;

  389.     }

  390.     /** Extract a sample of points from a mesh.
  391.      * @param mesh mesh from which grid should be extracted
  392.      * @param zone zone covered by the mesh
  393.      * @return extracted grid
  394.      */
  395.     private List<GeodeticPoint> extractSample(final Mesh mesh, final SphericalPolygonsSet zone) {

  396.         // find how to select sample points taking quantization into account
  397.         // to have the largest possible number of points while still
  398.         // being inside the zone of interest
  399.         int selectedAcrossModulus = -1;
  400.         int selectedAlongModulus  = -1;
  401.         int selectedCount         = -1;
  402.         for (int acrossModulus = 0; acrossModulus < quantization; ++acrossModulus) {
  403.             for (int alongModulus = 0; alongModulus < quantization; ++alongModulus) {

  404.                 // count how many points would be selected for the current modulus
  405.                 int count = 0;
  406.                 for (int across = mesh.getMinAcrossIndex() + acrossModulus;
  407.                         across <= mesh.getMaxAcrossIndex();
  408.                         across += quantization) {
  409.                     for (int along = mesh.getMinAlongIndex() + alongModulus;
  410.                             along <= mesh.getMaxAlongIndex();
  411.                             along += quantization) {
  412.                         final Mesh.Node  node = mesh.getNode(along, across);
  413.                         if (node != null && node.isInside()) {
  414.                             ++count;
  415.                         }
  416.                     }
  417.                 }

  418.                 if (count > selectedCount) {
  419.                     // current modulus are better than the selected ones
  420.                     selectedAcrossModulus = acrossModulus;
  421.                     selectedAlongModulus  = alongModulus;
  422.                     selectedCount         = count;
  423.                 }
  424.             }
  425.         }

  426.         // extract the sample points
  427.         final List<GeodeticPoint> sample = new ArrayList<GeodeticPoint>(selectedCount);
  428.         for (int across = mesh.getMinAcrossIndex() + selectedAcrossModulus;
  429.                 across <= mesh.getMaxAcrossIndex();
  430.                 across += quantization) {
  431.             for (int along = mesh.getMinAlongIndex() + selectedAlongModulus;
  432.                     along <= mesh.getMaxAlongIndex();
  433.                     along += quantization) {
  434.                 final Mesh.Node  node = mesh.getNode(along, across);
  435.                 if (node != null && node.isInside()) {
  436.                     sample.add(toGeodetic(node.getS2P()));
  437.                 }
  438.             }
  439.         }

  440.         return sample;

  441.     }

  442.     /** Merge two meshes together.
  443.      * @param mesh1 first mesh
  444.      * @param mesh2 second mesh
  445.      * @param mergingSeeds collection where to put the nodes created during the merge
  446.      * @return merged mesh (really one of the instances)
  447.      */
  448.     private Mesh mergeMeshes(final Mesh mesh1, final Mesh mesh2,
  449.                              final Collection<Mesh.Node> mergingSeeds) {

  450.         // select the way merge will be performed
  451.         final Mesh larger;
  452.         final Mesh smaller;
  453.         if (mesh1.getNumberOfNodes() >= mesh2.getNumberOfNodes()) {
  454.             // the larger new mesh should absorb the smaller existing mesh
  455.             larger  = mesh1;
  456.             smaller = mesh2;
  457.         } else {
  458.             // the larger existing mesh should absorb the smaller new mesh
  459.             larger  = mesh2;
  460.             smaller = mesh1;
  461.         }

  462.         // prepare seed nodes for next iteration
  463.         for (final Mesh.Node insideNode : smaller.getInsideNodes()) {

  464.             // beware we cannot reuse the node itself as the two meshes are not aligned!
  465.             // we have to create new nodes around the previous location
  466.             Mesh.Node node = larger.getClosestExistingNode(insideNode.getV());

  467.             while (estimateAlongMotion(node, insideNode.getV()) > +mesh1.getAlongGap()) {
  468.                 // the node is before desired index in the along direction
  469.                 // we need to create intermediates nodes up to the desired index
  470.                 node = larger.addNode(node.getAlongIndex() + 1, node.getAcrossIndex());
  471.             }

  472.             while (estimateAlongMotion(node, insideNode.getV()) < -mesh1.getAlongGap()) {
  473.                 // the node is after desired index in the along direction
  474.                 // we need to create intermediates nodes up to the desired index
  475.                 node = larger.addNode(node.getAlongIndex() - 1, node.getAcrossIndex());
  476.             }

  477.             while (estimateAcrossMotion(node, insideNode.getV()) > +mesh1.getAcrossGap()) {
  478.                 // the node is before desired index in the across direction
  479.                 // we need to create intermediates nodes up to the desired index
  480.                 node = larger.addNode(node.getAlongIndex(), node.getAcrossIndex() + 1);
  481.             }

  482.             while (estimateAcrossMotion(node, insideNode.getV()) < -mesh1.getAcrossGap()) {
  483.                 // the node is after desired index in the across direction
  484.                 // we need to create intermediates nodes up to the desired index
  485.                 node = larger.addNode(node.getAlongIndex(), node.getAcrossIndex() - 1);
  486.             }

  487.             // now we are close to the inside node,
  488.             // make sure the four surrounding nodes are available
  489.             final int otherAlong  = (estimateAlongMotion(node, insideNode.getV()) < 0.0) ?
  490.                                     node.getAlongIndex()  - 1 : node.getAlongIndex() + 1;
  491.             final int otherAcross = (estimateAcrossMotion(node, insideNode.getV()) < 0.0) ?
  492.                                     node.getAcrossIndex()  - 1 : node.getAcrossIndex() + 1;
  493.             addNode(node.getAlongIndex(), node.getAcrossIndex(), larger, mergingSeeds);
  494.             addNode(node.getAlongIndex(), otherAcross,           larger, mergingSeeds);
  495.             addNode(otherAlong,           node.getAcrossIndex(), larger, mergingSeeds);
  496.             addNode(otherAlong,           otherAcross,           larger, mergingSeeds);

  497.         }

  498.         return larger;

  499.     }

  500.     /** Ensure all 8 neighbors of a node are in the mesh.
  501.      * @param base base node
  502.      * @param mesh complete mesh containing nodes
  503.      * @param newNodes queue where new node must be put
  504.      */
  505.     private void addAllNeighborsIfNeeded(final Mesh.Node base, final Mesh mesh,
  506.                                          final Collection<Mesh.Node> newNodes) {
  507.         addNode(base.getAlongIndex() - 1, base.getAcrossIndex() - 1, mesh, newNodes);
  508.         addNode(base.getAlongIndex() - 1, base.getAcrossIndex(),     mesh, newNodes);
  509.         addNode(base.getAlongIndex() - 1, base.getAcrossIndex() + 1, mesh, newNodes);
  510.         addNode(base.getAlongIndex(),     base.getAcrossIndex() - 1, mesh, newNodes);
  511.         addNode(base.getAlongIndex(),     base.getAcrossIndex() + 1, mesh, newNodes);
  512.         addNode(base.getAlongIndex() + 1, base.getAcrossIndex() - 1, mesh, newNodes);
  513.         addNode(base.getAlongIndex() + 1, base.getAcrossIndex(),     mesh, newNodes);
  514.         addNode(base.getAlongIndex() + 1, base.getAcrossIndex() + 1, mesh, newNodes);
  515.     }

  516.     /** Add a node to a mesh if not already present.
  517.      * @param alongIndex index in the along direction
  518.      * @param acrossIndex index in the across direction
  519.      * @param mesh complete mesh containing nodes
  520.      * @param newNodes queue where new node must be put
  521.      */
  522.     private void addNode(final int alongIndex, final int acrossIndex,
  523.                          final Mesh mesh, final Collection<Mesh.Node> newNodes) {

  524.         final Mesh.Node node = mesh.addNode(alongIndex, acrossIndex);

  525.         if (!node.isEnabled()) {
  526.             // enable the node
  527.             node.setEnabled();
  528.             newNodes.add(node);
  529.         }

  530.     }

  531.     /** Convert a point on the unit 2-sphere to geodetic coordinates.
  532.      * @param point point on the unit 2-sphere
  533.      * @return geodetic point (arbitrarily set at altitude 0)
  534.      */
  535.     protected GeodeticPoint toGeodetic(final S2Point point) {
  536.         return new GeodeticPoint(0.5 * FastMath.PI - point.getPhi(), point.getTheta(), 0.0);
  537.     }

  538.     /** Build a simple zone (connected zone without holes).
  539.      * <p>
  540.      * In order to build more complex zones (not connected or with
  541.      * holes), the user should directly call Hipparchus
  542.      * {@link SphericalPolygonsSet} constructors and
  543.      * {@link RegionFactory region factory} if set operations
  544.      * are needed (union, intersection, difference ...).
  545.      * </p>
  546.      * <p>
  547.      * Take care that the vertices boundary points must be given <em>counterclockwise</em>.
  548.      * Using the wrong order defines the complementary of the real zone,
  549.      * and will often result in tessellation failure as the zone is too
  550.      * wide.
  551.      * </p>
  552.      * @param tolerance angular separation below which points are considered
  553.      * equal (typically 1.0e-10)
  554.      * @param points vertices of the boundary, in <em>counterclockwise</em>
  555.      * order, each point being a two-elements arrays with latitude at index 0
  556.      * and longitude at index 1
  557.      * @return a zone defined on the unit 2-sphere
  558.      */
  559.     public static SphericalPolygonsSet buildSimpleZone(final double tolerance,
  560.                                                        final double[]... points) {
  561.         final S2Point[] vertices = new S2Point[points.length];
  562.         for (int i = 0; i < points.length; ++i) {
  563.             vertices[i] = new S2Point(points[i][1], 0.5 * FastMath.PI - points[i][0]);
  564.         }
  565.         return new SphericalPolygonsSet(tolerance, vertices);
  566.     }

  567.     /** Build a simple zone (connected zone without holes).
  568.      * <p>
  569.      * In order to build more complex zones (not connected or with
  570.      * holes), the user should directly call Hipparchus
  571.      * {@link SphericalPolygonsSet} constructors and
  572.      * {@link RegionFactory region factory} if set operations
  573.      * are needed (union, intersection, difference ...).
  574.      * </p>
  575.      * <p>
  576.      * Take care that the vertices boundary points must be given <em>counterclockwise</em>.
  577.      * Using the wrong order defines the complementary of the real zone,
  578.      * and will often result in tessellation failure as the zone is too
  579.      * wide.
  580.      * </p>
  581.      * @param tolerance angular separation below which points are considered
  582.      * equal (typically 1.0e-10)
  583.      * @param points vertices of the boundary, in <em>counterclockwise</em>
  584.      * order
  585.      * @return a zone defined on the unit 2-sphere
  586.      */
  587.     public static SphericalPolygonsSet buildSimpleZone(final double tolerance,
  588.                                                        final GeodeticPoint... points) {
  589.         final S2Point[] vertices = new S2Point[points.length];
  590.         for (int i = 0; i < points.length; ++i) {
  591.             vertices[i] = new S2Point(points[i].getLongitude(),
  592.                                       0.5 * FastMath.PI - points[i].getLatitude());
  593.         }
  594.         return new SphericalPolygonsSet(tolerance, vertices);
  595.     }

  596.     /** Estimate an approximate motion in the along direction.
  597.      * @param start node at start of motion
  598.      * @param end desired point at end of motion
  599.      * @return approximate motion in the along direction
  600.      */
  601.     private double estimateAlongMotion(final Mesh.Node start, final Vector3D end) {
  602.         return Vector3D.dotProduct(start.getAlong(), end.subtract(start.getV()));
  603.     }

  604.     /** Estimate an approximate motion in the across direction.
  605.      * @param start node at start of motion
  606.      * @param end desired point at end of motion
  607.      * @return approximate motion in the across direction
  608.      */
  609.     private double estimateAcrossMotion(final Mesh.Node start, final Vector3D end) {
  610.         return Vector3D.dotProduct(start.getAcross(), end.subtract(start.getV()));
  611.     }

  612.     /** Check if an arc meets the inside of a zone.
  613.      * @param s1 first point
  614.      * @param s2 second point
  615.      * @param zone zone to check arc against
  616.      * @return true if the arc meets the inside of the zone
  617.      */
  618.     private boolean meetInside(final S2Point s1, final S2Point s2,
  619.                                final SphericalPolygonsSet zone) {
  620.         final Circle  circle = new Circle(s1, s2, zone.getTolerance());
  621.         final double alpha1  = circle.toSubSpace(s1).getAlpha();
  622.         final double alpha2  = MathUtils.normalizeAngle(circle.toSubSpace(s2).getAlpha(),
  623.                                                         alpha1 + FastMath.PI);
  624.         final SubCircle sub  = new SubCircle(circle,
  625.                                              new ArcsSet(alpha1, alpha2, zone.getTolerance()));
  626.         return recurseMeetInside(zone.getTree(false), sub);

  627.     }

  628.     /** Check if an arc meets the inside of a zone.
  629.      * <p>
  630.      * This method is heavily based on the Characterization class from
  631.      * Hipparchus library, also distributed under the terms
  632.      * of the Apache Software License V2.
  633.      * </p>
  634.      * @param node spherical zone node
  635.      * @param sub arc to characterize
  636.      * @return true if the arc meets the inside of the zone
  637.      */
  638.     private boolean recurseMeetInside(final BSPTree<Sphere2D> node, final SubHyperplane<Sphere2D> sub) {

  639.         if (node.getCut() == null) {
  640.             // we have reached a leaf node
  641.             if (sub.isEmpty()) {
  642.                 return false;
  643.             } else {
  644.                 return (Boolean) node.getAttribute();
  645.             }
  646.         } else {
  647.             final Hyperplane<Sphere2D> hyperplane = node.getCut().getHyperplane();
  648.             final SubHyperplane.SplitSubHyperplane<Sphere2D> split = sub.split(hyperplane);
  649.             switch (split.getSide()) {
  650.                 case PLUS:
  651.                     return recurseMeetInside(node.getPlus(),  sub);
  652.                 case MINUS:
  653.                     return recurseMeetInside(node.getMinus(), sub);
  654.                 case BOTH:
  655.                     if (recurseMeetInside(node.getPlus(), split.getPlus())) {
  656.                         return true;
  657.                     } else {
  658.                         return recurseMeetInside(node.getMinus(), split.getMinus());
  659.                     }
  660.                 default:
  661.                     // this should not happen
  662.                     throw new OrekitInternalError(null);
  663.             }
  664.         }
  665.     }

  666.     /** Get an iterator over mesh nodes indices.
  667.      * @param minIndex minimum node index
  668.      * @param maxIndex maximum node index
  669.      * @param truncateLast true if we can reduce last tile
  670.      * @return iterator over mesh nodes indices
  671.      */
  672.     private Iterable<Range> nodesIndices(final int minIndex, final int maxIndex, final boolean truncateLast) {

  673.         final int first;
  674.         if (truncateLast) {

  675.             // truncate last tile rather than balance tiles around the zone of interest
  676.             first = minIndex;

  677.         } else {

  678.             // balance tiles around the zone of interest rather than truncate last tile

  679.             // number of tiles needed to cover the full indices range
  680.             final int range      = maxIndex - minIndex;
  681.             final int nbTiles    = (range + quantization - 1) / quantization;

  682.             // extra nodes that must be added to complete the tiles
  683.             final int extraNodes = nbTiles * quantization  - range;

  684.             // balance the extra nodes before min index and after maxIndex
  685.             final int extraBefore = (extraNodes + 1) / 2;

  686.             first = minIndex - extraBefore;

  687.         }

  688.         return new Iterable<Range>() {

  689.             /** {@inheritDoc} */
  690.             @Override
  691.             public Iterator<Range> iterator() {
  692.                 return new Iterator<Range>() {

  693.                     private int nextLower = first;

  694.                     /** {@inheritDoc} */
  695.                     @Override
  696.                     public boolean hasNext() {
  697.                         return nextLower < maxIndex;
  698.                     }

  699.                     /** {@inheritDoc} */
  700.                     @Override
  701.                     public Range next() {

  702.                         if (nextLower >= maxIndex) {
  703.                             throw new NoSuchElementException();
  704.                         }
  705.                         final int lower = nextLower;

  706.                         nextLower += quantization;
  707.                         if (truncateLast && nextLower > maxIndex && lower < maxIndex) {
  708.                             // truncate last tile
  709.                             nextLower = maxIndex;
  710.                         }

  711.                         return new Range(lower, nextLower);

  712.                     }

  713.                     /** {@inheritDoc} */
  714.                     @Override
  715.                     public void remove() {
  716.                         throw new UnsupportedOperationException();
  717.                     }

  718.                 };
  719.             }
  720.         };

  721.     }

  722.     /** Local class for a range of indices to be used for building a tile. */
  723.     private static class Range {

  724.         /** Lower index. */
  725.         private final int lower;

  726.         /** Upper index. */
  727.         private final int upper;

  728.         /** Simple constructor.
  729.          * @param lower lower index
  730.          * @param upper upper index
  731.          */
  732.         Range(final int lower, final int upper) {
  733.             this.lower = lower;
  734.             this.upper = upper;
  735.         }

  736.     }

  737.     /** Local class for a pair of ranges of indices to be used for building a tile. */
  738.     private static class RangePair {

  739.         /** Across range. */
  740.         private final Range across;

  741.         /** Along range. */
  742.         private final Range along;

  743.         /** Simple constructor.
  744.          * @param across across range
  745.          * @param along along range
  746.          */
  747.         RangePair(final Range across, final Range along) {
  748.             this.across = across;
  749.             this.along  = along;
  750.         }

  751.     }

  752. }