EllipsoidTessellator.java

  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. 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.geometry.euclidean.threed.Vector3D;
  28. import org.hipparchus.geometry.partitioning.BSPTree;
  29. import org.hipparchus.geometry.partitioning.Hyperplane;
  30. import org.hipparchus.geometry.partitioning.RegionFactory;
  31. import org.hipparchus.geometry.partitioning.SubHyperplane;
  32. import org.hipparchus.geometry.spherical.oned.ArcsSet;
  33. import org.hipparchus.geometry.spherical.twod.Circle;
  34. import org.hipparchus.geometry.spherical.twod.S2Point;
  35. import org.hipparchus.geometry.spherical.twod.Sphere2D;
  36. import org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet;
  37. import org.hipparchus.geometry.spherical.twod.SubCircle;
  38. import org.hipparchus.util.FastMath;
  39. import org.hipparchus.util.MathUtils;
  40. import org.orekit.bodies.GeodeticPoint;
  41. import org.orekit.bodies.OneAxisEllipsoid;
  42. import org.orekit.errors.OrekitInternalError;

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

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

  66.     /** Aiming used for orienting tiles. */
  67.     private final TileAiming aiming;

  68.     /** Underlying ellipsoid. */
  69.     private final OneAxisEllipsoid ellipsoid;

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

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

  151.         final double                  splitWidth  = (fullWidth  - widthOverlap)  / quantization;
  152.         final double                  splitLength = (fullLength - lengthOverlap) / quantization;
  153.         final Map<Mesh, List<Tile>>   map         = new IdentityHashMap<Mesh, List<Tile>>();
  154.         final RegionFactory<Sphere2D> factory     = new RegionFactory<Sphere2D>();
  155.         SphericalPolygonsSet          remaining   = (SphericalPolygonsSet) zone.copySelf();
  156.         S2Point                       inside      = getInsidePoint(remaining);

  157.         while (inside != null) {

  158.             // find a mesh covering at least one connected part of the zone
  159.             final List<Mesh.Node> mergingSeeds = new ArrayList<Mesh.Node>();
  160.             Mesh mesh = new Mesh(ellipsoid, zone, aiming, splitLength, splitWidth, inside);
  161.             mergingSeeds.add(mesh.getNode(0, 0));
  162.             List<Tile> tiles = null;
  163.             while (!mergingSeeds.isEmpty()) {

  164.                 // expand the mesh around the seed
  165.                 neighborExpandMesh(mesh, mergingSeeds, zone);

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

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

  175.                         // merge the two meshes together
  176.                         mesh = mergeMeshes(mesh, entry.getKey(), mergingSeeds);
  177.                         map.remove(entry.getKey());
  178.                         break;

  179.                     }
  180.                 }

  181.             }

  182.             // remove the part of the zone covered by the mesh
  183.             remaining = (SphericalPolygonsSet) factory.difference(remaining, mesh.getCoverage());
  184.             inside    = getInsidePoint(remaining);

  185.             map.put(mesh, tiles);

  186.         }

  187.         // concatenate the lists from the independent meshes
  188.         final List<List<Tile>> tilesLists = new ArrayList<List<Tile>>(map.size());
  189.         for (final Map.Entry<Mesh, List<Tile>> entry : map.entrySet()) {
  190.             tilesLists.add(entry.getValue());
  191.         }

  192.         return tilesLists;

  193.     }

  194.     /** Sample a zone of interest into a grid sample of {@link GeodeticPoint geodetic points}.
  195.      * <p>
  196.      * The created points will be entirely within the zone of interest.
  197.      * </p>
  198.      * @param zone zone of interest to sample
  199.      * @param width grid sample cells width as a distance on surface (in meters)
  200.      * @param length grid sample cells length as a distance on surface (in meters)
  201.      * @return a list of lists of points sampling the zone of interest,
  202.      * each sub-list corresponding to a part not connected to the other
  203.      * parts (for example for islands)
  204.      */
  205.     public List<List<GeodeticPoint>> sample(final SphericalPolygonsSet zone,
  206.                                             final double width, final double length) {

  207.         final double                         splitWidth  = width  / quantization;
  208.         final double                         splitLength = length / quantization;
  209.         final Map<Mesh, List<GeodeticPoint>> map         = new IdentityHashMap<Mesh, List<GeodeticPoint>>();
  210.         final RegionFactory<Sphere2D>        factory     = new RegionFactory<Sphere2D>();
  211.         SphericalPolygonsSet                 remaining   = (SphericalPolygonsSet) zone.copySelf();
  212.         S2Point                              inside      = getInsidePoint(remaining);

  213.         while (inside != null) {

  214.             // find a mesh covering at least one connected part of the zone
  215.             final List<Mesh.Node> mergingSeeds = new ArrayList<Mesh.Node>();
  216.             Mesh mesh = new Mesh(ellipsoid, zone, aiming, splitLength, splitWidth, inside);
  217.             mergingSeeds.add(mesh.getNode(0, 0));
  218.             List<GeodeticPoint> sample = null;
  219.             while (!mergingSeeds.isEmpty()) {

  220.                 // expand the mesh around the seed
  221.                 neighborExpandMesh(mesh, mergingSeeds, zone);

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

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

  231.                         // merge the two meshes together
  232.                         mesh = mergeMeshes(mesh, entry.getKey(), mergingSeeds);
  233.                         map.remove(entry.getKey());
  234.                         break;

  235.                     }
  236.                 }

  237.             }

  238.             // remove the part of the zone covered by the mesh
  239.             remaining = (SphericalPolygonsSet) factory.difference(remaining, mesh.getCoverage());
  240.             inside    = getInsidePoint(remaining);

  241.             map.put(mesh, sample);

  242.         }

  243.         // concatenate the lists from the independent meshes
  244.         final List<List<GeodeticPoint>> sampleLists = new ArrayList<List<GeodeticPoint>>(map.size());
  245.         for (final Map.Entry<Mesh, List<GeodeticPoint>> entry : map.entrySet()) {
  246.             sampleLists.add(entry.getValue());
  247.         }

  248.         return sampleLists;

  249.     }

  250.     /** Get an inside point from a zone of interest.
  251.      * @param zone zone to mesh
  252.      * @return a point inside the zone or null if zone is empty or too thin
  253.      */
  254.     private S2Point getInsidePoint(final SphericalPolygonsSet zone) {

  255.         final InsideFinder finder = new InsideFinder(zone);
  256.         zone.getTree(false).visit(finder);
  257.         return finder.getInsidePoint();

  258.     }

  259.     /** Expand a mesh so it surrounds at least one connected part of a zone.
  260.      * <p>
  261.      * This part of mesh expansion is neighbors based. It includes the seed
  262.      * node neighbors, and their neighbors, and the neighbors of their
  263.      * neighbors until the path-connected sub-parts of the zone these nodes
  264.      * belong to are completely surrounded by the mesh taxicab boundary.
  265.      * </p>
  266.      * @param mesh mesh to expand
  267.      * @param seeds seed nodes (already in the mesh) from which to start expansion
  268.      * @param zone zone to mesh
  269.      */
  270.     private void neighborExpandMesh(final Mesh mesh, final Collection<Mesh.Node> seeds,
  271.                                     final SphericalPolygonsSet zone) {

  272.         // mesh expansion loop
  273.         boolean expanding = true;
  274.         final Queue<Mesh.Node> newNodes = new LinkedList<Mesh.Node>();
  275.         newNodes.addAll(seeds);
  276.         while (expanding) {

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

  281.                 // retrieve an active node
  282.                 final Mesh.Node node = newNodes.remove();

  283.                 if (node.isInside()) {
  284.                     // the node is inside the zone, the mesh must contain its 8 neighbors
  285.                     addAllNeighborsIfNeeded(node, mesh, newNodes);
  286.                 }

  287.             }

  288.             // second expansion step: check if the loop of outside nodes
  289.             // completely surrounds the zone, i.e. there are no peaks
  290.             // pointing out of the loop between two nodes
  291.             expanding = false;
  292.             final List<Mesh.Node> boundary = mesh.getTaxicabBoundary(false);
  293.             if (boundary.size() > 1) {
  294.                 Mesh.Node previous = boundary.get(boundary.size() - 1);
  295.                 for (final Mesh.Node node : boundary) {
  296.                     if (meetInside(previous.getS2P(), node.getS2P(), zone)) {
  297.                         // part of the mesh boundary is still inside the zone!
  298.                         // the mesh must be expanded again
  299.                         addAllNeighborsIfNeeded(previous, mesh, newNodes);
  300.                         addAllNeighborsIfNeeded(node,     mesh, newNodes);
  301.                         expanding = true;
  302.                     }
  303.                     previous = node;
  304.                 }
  305.             }

  306.         }

  307.     }

  308.     /** Extract tiles from a mesh.
  309.      * @param mesh mesh from which tiles should be extracted
  310.      * @param zone zone covered by the mesh
  311.      * @param lengthOverlap overlap between adjacent tiles
  312.      * @param widthOverlap overlap between adjacent tiles
  313.      * @param truncateLastWidth true if we can reduce last tile width
  314.      * @param truncateLastLength true if we can reduce last tile length
  315.      * @return extracted tiles
  316.      */
  317.     private List<Tile> extractTiles(final Mesh mesh, final SphericalPolygonsSet zone,
  318.                                     final double lengthOverlap, final double widthOverlap,
  319.                                     final boolean truncateLastWidth, final boolean truncateLastLength) {

  320.         final List<Tile>      tiles = new ArrayList<Tile>();
  321.         final List<RangePair> rangePairs = new ArrayList<RangePair>();

  322.         final int minAcross = mesh.getMinAcrossIndex();
  323.         final int maxAcross = mesh.getMaxAcrossIndex();
  324.         for (Range acrossPair : nodesIndices(minAcross, maxAcross, truncateLastWidth)) {

  325.             int minAlong = mesh.getMaxAlongIndex() + 1;
  326.             int maxAlong = mesh.getMinAlongIndex() - 1;
  327.             for (int c = acrossPair.lower; c <= acrossPair.upper; ++c) {
  328.                 minAlong = FastMath.min(minAlong, mesh.getMinAlongIndex(c));
  329.                 maxAlong = FastMath.max(maxAlong, mesh.getMaxAlongIndex(c));
  330.             }

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

  332.                 // get the base vertex nodes
  333.                 final Mesh.Node node0 = mesh.addNode(alongPair.lower, acrossPair.lower);
  334.                 final Mesh.Node node1 = mesh.addNode(alongPair.upper, acrossPair.lower);
  335.                 final Mesh.Node node2 = mesh.addNode(alongPair.upper, acrossPair.upper);
  336.                 final Mesh.Node node3 = mesh.addNode(alongPair.lower, acrossPair.upper);

  337.                 // apply tile overlap
  338.                 final S2Point s2p0 = node0.move(new Vector3D(-0.5 * lengthOverlap, node0.getAlong(),
  339.                                                              -0.5 * widthOverlap,  node0.getAcross()));
  340.                 final S2Point s2p1 = node1.move(new Vector3D(+0.5 * lengthOverlap, node1.getAlong(),
  341.                                                              -0.5 * widthOverlap,  node1.getAcross()));
  342.                 final S2Point s2p2 = node2.move(new Vector3D(+0.5 * lengthOverlap, node2.getAlong(),
  343.                                                              +0.5 * widthOverlap,  node2.getAcross()));
  344.                 final S2Point s2p3 = node3.move(new Vector3D(-0.5 * lengthOverlap, node2.getAlong(),
  345.                                                              +0.5 * widthOverlap,  node2.getAcross()));

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

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

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

  353.                 }

  354.             }
  355.         }

  356.         // ensure the taxicab boundary follows the built tile sides
  357.         // this is done outside of the previous loop in order
  358.         // to avoid one tile changing the min/max indices of the
  359.         // neighboring tile as they share some nodes that will be enabled here
  360.         for (final RangePair rangePair : rangePairs) {
  361.             for (int c = rangePair.across.lower; c < rangePair.across.upper; ++c) {
  362.                 mesh.addNode(rangePair.along.lower, c + 1).setEnabled();
  363.                 mesh.addNode(rangePair.along.upper, c).setEnabled();
  364.             }
  365.             for (int l = rangePair.along.lower; l < rangePair.along.upper; ++l) {
  366.                 mesh.addNode(l,     rangePair.across.lower).setEnabled();
  367.                 mesh.addNode(l + 1, rangePair.across.upper).setEnabled();
  368.             }
  369.         }

  370.         return tiles;

  371.     }

  372.     /** Extract a sample of points from a mesh.
  373.      * @param mesh mesh from which grid should be extracted
  374.      * @param zone zone covered by the mesh
  375.      * @return extracted grid
  376.      */
  377.     private List<GeodeticPoint> extractSample(final Mesh mesh, final SphericalPolygonsSet zone) {

  378.         // find how to select sample points taking quantization into account
  379.         // to have the largest possible number of points while still
  380.         // being inside the zone of interest
  381.         int selectedAcrossModulus = -1;
  382.         int selectedAlongModulus  = -1;
  383.         int selectedCount         = -1;
  384.         for (int acrossModulus = 0; acrossModulus < quantization; ++acrossModulus) {
  385.             for (int alongModulus = 0; alongModulus < quantization; ++alongModulus) {

  386.                 // count how many points would be selected for the current modulus
  387.                 int count = 0;
  388.                 for (int across = mesh.getMinAcrossIndex() + acrossModulus;
  389.                         across <= mesh.getMaxAcrossIndex();
  390.                         across += quantization) {
  391.                     for (int along = mesh.getMinAlongIndex() + alongModulus;
  392.                             along <= mesh.getMaxAlongIndex();
  393.                             along += quantization) {
  394.                         final Mesh.Node  node = mesh.getNode(along, across);
  395.                         if (node != null && node.isInside()) {
  396.                             ++count;
  397.                         }
  398.                     }
  399.                 }

  400.                 if (count > selectedCount) {
  401.                     // current modulus are better than the selected ones
  402.                     selectedAcrossModulus = acrossModulus;
  403.                     selectedAlongModulus  = alongModulus;
  404.                     selectedCount         = count;
  405.                 }
  406.             }
  407.         }

  408.         // extract the sample points
  409.         final List<GeodeticPoint> sample = new ArrayList<GeodeticPoint>(selectedCount);
  410.         for (int across = mesh.getMinAcrossIndex() + selectedAcrossModulus;
  411.                 across <= mesh.getMaxAcrossIndex();
  412.                 across += quantization) {
  413.             for (int along = mesh.getMinAlongIndex() + selectedAlongModulus;
  414.                     along <= mesh.getMaxAlongIndex();
  415.                     along += quantization) {
  416.                 final Mesh.Node  node = mesh.getNode(along, across);
  417.                 if (node != null && node.isInside()) {
  418.                     sample.add(toGeodetic(node.getS2P()));
  419.                 }
  420.             }
  421.         }

  422.         return sample;

  423.     }

  424.     /** Merge two meshes together.
  425.      * @param mesh1 first mesh
  426.      * @param mesh2 second mesh
  427.      * @param mergingSeeds collection where to put the nodes created during the merge
  428.      * @return merged mesh (really one of the instances)
  429.      */
  430.     private Mesh mergeMeshes(final Mesh mesh1, final Mesh mesh2,
  431.                              final Collection<Mesh.Node> mergingSeeds) {

  432.         // select the way merge will be performed
  433.         final Mesh larger;
  434.         final Mesh smaller;
  435.         if (mesh1.getNumberOfNodes() >= mesh2.getNumberOfNodes()) {
  436.             // the larger new mesh should absorb the smaller existing mesh
  437.             larger  = mesh1;
  438.             smaller = mesh2;
  439.         } else {
  440.             // the larger existing mesh should absorb the smaller new mesh
  441.             larger  = mesh2;
  442.             smaller = mesh1;
  443.         }

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

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

  449.             while (estimateAlongMotion(node, insideNode.getV()) > +mesh1.getAlongGap()) {
  450.                 // the node is before desired index in the along direction
  451.                 // we need to create intermediates nodes up to the desired index
  452.                 node = larger.addNode(node.getAlongIndex() + 1, node.getAcrossIndex());
  453.             }

  454.             while (estimateAlongMotion(node, insideNode.getV()) < -mesh1.getAlongGap()) {
  455.                 // the node is after desired index in the along direction
  456.                 // we need to create intermediates nodes up to the desired index
  457.                 node = larger.addNode(node.getAlongIndex() - 1, node.getAcrossIndex());
  458.             }

  459.             while (estimateAcrossMotion(node, insideNode.getV()) > +mesh1.getAcrossGap()) {
  460.                 // the node is before desired index in the across direction
  461.                 // we need to create intermediates nodes up to the desired index
  462.                 node = larger.addNode(node.getAlongIndex(), node.getAcrossIndex() + 1);
  463.             }

  464.             while (estimateAcrossMotion(node, insideNode.getV()) < -mesh1.getAcrossGap()) {
  465.                 // the node is after desired index in the across direction
  466.                 // we need to create intermediates nodes up to the desired index
  467.                 node = larger.addNode(node.getAlongIndex(), node.getAcrossIndex() - 1);
  468.             }

  469.             // now we are close to the inside node,
  470.             // make sure the four surrounding nodes are available
  471.             final int otherAlong  = (estimateAlongMotion(node, insideNode.getV()) < 0.0) ?
  472.                                     node.getAlongIndex()  - 1 : node.getAlongIndex() + 1;
  473.             final int otherAcross = (estimateAcrossMotion(node, insideNode.getV()) < 0.0) ?
  474.                                     node.getAcrossIndex()  - 1 : node.getAcrossIndex() + 1;
  475.             addNode(node.getAlongIndex(), node.getAcrossIndex(), larger, mergingSeeds);
  476.             addNode(node.getAlongIndex(), otherAcross,           larger, mergingSeeds);
  477.             addNode(otherAlong,           node.getAcrossIndex(), larger, mergingSeeds);
  478.             addNode(otherAlong,           otherAcross,           larger, mergingSeeds);

  479.         }

  480.         return larger;

  481.     }

  482.     /** Ensure all 8 neighbors of a node are in the mesh.
  483.      * @param base base node
  484.      * @param mesh complete mesh containing nodes
  485.      * @param newNodes queue where new node must be put
  486.      */
  487.     private void addAllNeighborsIfNeeded(final Mesh.Node base, final Mesh mesh,
  488.                                          final Collection<Mesh.Node> newNodes) {
  489.         addNode(base.getAlongIndex() - 1, base.getAcrossIndex() - 1, mesh, newNodes);
  490.         addNode(base.getAlongIndex() - 1, base.getAcrossIndex(),     mesh, newNodes);
  491.         addNode(base.getAlongIndex() - 1, base.getAcrossIndex() + 1, mesh, newNodes);
  492.         addNode(base.getAlongIndex(),     base.getAcrossIndex() - 1, mesh, newNodes);
  493.         addNode(base.getAlongIndex(),     base.getAcrossIndex() + 1, mesh, newNodes);
  494.         addNode(base.getAlongIndex() + 1, base.getAcrossIndex() - 1, mesh, newNodes);
  495.         addNode(base.getAlongIndex() + 1, base.getAcrossIndex(),     mesh, newNodes);
  496.         addNode(base.getAlongIndex() + 1, base.getAcrossIndex() + 1, mesh, newNodes);
  497.     }

  498.     /** Add a node to a mesh if not already present.
  499.      * @param alongIndex index in the along direction
  500.      * @param acrossIndex index in the across direction
  501.      * @param mesh complete mesh containing nodes
  502.      * @param newNodes queue where new node must be put
  503.      */
  504.     private void addNode(final int alongIndex, final int acrossIndex,
  505.                          final Mesh mesh, final Collection<Mesh.Node> newNodes) {

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

  507.         if (!node.isEnabled()) {
  508.             // enable the node
  509.             node.setEnabled();
  510.             newNodes.add(node);
  511.         }

  512.     }

  513.     /** Convert a point on the unit 2-sphere to geodetic coordinates.
  514.      * @param point point on the unit 2-sphere
  515.      * @return geodetic point (arbitrarily set at altitude 0)
  516.      */
  517.     protected GeodeticPoint toGeodetic(final S2Point point) {
  518.         return new GeodeticPoint(0.5 * FastMath.PI - point.getPhi(), point.getTheta(), 0.0);
  519.     }

  520.     /** Build a simple zone (connected zone without holes).
  521.      * <p>
  522.      * In order to build more complex zones (not connected or with
  523.      * holes), the user should directly call Hipparchus
  524.      * {@link SphericalPolygonsSet} constructors and
  525.      * {@link RegionFactory region factory} if set operations
  526.      * are needed (union, intersection, difference ...).
  527.      * </p>
  528.      * <p>
  529.      * Take care that the vertices boundary points must be given <em>counterclockwise</em>.
  530.      * Using the wrong order defines the complementary of the real zone,
  531.      * and will often result in tessellation failure as the zone is too
  532.      * wide.
  533.      * </p>
  534.      * @param tolerance angular separation below which points are considered
  535.      * equal (typically 1.0e-10)
  536.      * @param points vertices of the boundary, in <em>counterclockwise</em>
  537.      * order, each point being a two-elements arrays with latitude at index 0
  538.      * and longitude at index 1
  539.      * @return a zone defined on the unit 2-sphere
  540.      */
  541.     public static SphericalPolygonsSet buildSimpleZone(final double tolerance,
  542.                                                        final double[]... points) {
  543.         final S2Point[] vertices = new S2Point[points.length];
  544.         for (int i = 0; i < points.length; ++i) {
  545.             vertices[i] = new S2Point(points[i][1], 0.5 * FastMath.PI - points[i][0]);
  546.         }
  547.         return new SphericalPolygonsSet(tolerance, vertices);
  548.     }

  549.     /** Build a simple zone (connected zone without holes).
  550.      * <p>
  551.      * In order to build more complex zones (not connected or with
  552.      * holes), the user should directly call Hipparchus
  553.      * {@link SphericalPolygonsSet} constructors and
  554.      * {@link RegionFactory region factory} if set operations
  555.      * are needed (union, intersection, difference ...).
  556.      * </p>
  557.      * <p>
  558.      * Take care that the vertices boundary points must be given <em>counterclockwise</em>.
  559.      * Using the wrong order defines the complementary of the real zone,
  560.      * and will often result in tessellation failure as the zone is too
  561.      * wide.
  562.      * </p>
  563.      * @param tolerance angular separation below which points are considered
  564.      * equal (typically 1.0e-10)
  565.      * @param points vertices of the boundary, in <em>counterclockwise</em>
  566.      * order
  567.      * @return a zone defined on the unit 2-sphere
  568.      */
  569.     public static SphericalPolygonsSet buildSimpleZone(final double tolerance,
  570.                                                        final GeodeticPoint... points) {
  571.         final S2Point[] vertices = new S2Point[points.length];
  572.         for (int i = 0; i < points.length; ++i) {
  573.             vertices[i] = new S2Point(points[i].getLongitude(),
  574.                                       0.5 * FastMath.PI - points[i].getLatitude());
  575.         }
  576.         return new SphericalPolygonsSet(tolerance, vertices);
  577.     }

  578.     /** Estimate an approximate motion in the along direction.
  579.      * @param start node at start of motion
  580.      * @param end desired point at end of motion
  581.      * @return approximate motion in the along direction
  582.      */
  583.     private double estimateAlongMotion(final Mesh.Node start, final Vector3D end) {
  584.         return Vector3D.dotProduct(start.getAlong(), end.subtract(start.getV()));
  585.     }

  586.     /** Estimate an approximate motion in the across direction.
  587.      * @param start node at start of motion
  588.      * @param end desired point at end of motion
  589.      * @return approximate motion in the across direction
  590.      */
  591.     private double estimateAcrossMotion(final Mesh.Node start, final Vector3D end) {
  592.         return Vector3D.dotProduct(start.getAcross(), end.subtract(start.getV()));
  593.     }

  594.     /** Check if an arc meets the inside of a zone.
  595.      * @param s1 first point
  596.      * @param s2 second point
  597.      * @param zone zone to check arc against
  598.      * @return true if the arc meets the inside of the zone
  599.      */
  600.     private boolean meetInside(final S2Point s1, final S2Point s2,
  601.                                final SphericalPolygonsSet zone) {
  602.         final Circle  circle = new Circle(s1, s2, zone.getTolerance());
  603.         final double alpha1  = circle.toSubSpace(s1).getAlpha();
  604.         final double alpha2  = MathUtils.normalizeAngle(circle.toSubSpace(s2).getAlpha(),
  605.                                                         alpha1 + FastMath.PI);
  606.         final SubCircle sub  = new SubCircle(circle,
  607.                                              new ArcsSet(alpha1, alpha2, zone.getTolerance()));
  608.         return recurseMeetInside(zone.getTree(false), sub);

  609.     }

  610.     /** Check if an arc meets the inside of a zone.
  611.      * <p>
  612.      * This method is heavily based on the Characterization class from
  613.      * Hipparchus library, also distributed under the terms
  614.      * of the Apache Software License V2.
  615.      * </p>
  616.      * @param node spherical zone node
  617.      * @param sub arc to characterize
  618.      * @return true if the arc meets the inside of the zone
  619.      */
  620.     private boolean recurseMeetInside(final BSPTree<Sphere2D> node, final SubHyperplane<Sphere2D> sub) {

  621.         if (node.getCut() == null) {
  622.             // we have reached a leaf node
  623.             if (sub.isEmpty()) {
  624.                 return false;
  625.             } else {
  626.                 return (Boolean) node.getAttribute();
  627.             }
  628.         } else {
  629.             final Hyperplane<Sphere2D> hyperplane = node.getCut().getHyperplane();
  630.             final SubHyperplane.SplitSubHyperplane<Sphere2D> split = sub.split(hyperplane);
  631.             switch (split.getSide()) {
  632.                 case PLUS:
  633.                     return recurseMeetInside(node.getPlus(),  sub);
  634.                 case MINUS:
  635.                     return recurseMeetInside(node.getMinus(), sub);
  636.                 case BOTH:
  637.                     if (recurseMeetInside(node.getPlus(), split.getPlus())) {
  638.                         return true;
  639.                     } else {
  640.                         return recurseMeetInside(node.getMinus(), split.getMinus());
  641.                     }
  642.                 default:
  643.                     // this should not happen
  644.                     throw new OrekitInternalError(null);
  645.             }
  646.         }
  647.     }

  648.     /** Get an iterator over mesh nodes indices.
  649.      * @param minIndex minimum node index
  650.      * @param maxIndex maximum node index
  651.      * @param truncateLast true if we can reduce last tile
  652.      * @return iterator over mesh nodes indices
  653.      */
  654.     private Iterable<Range> nodesIndices(final int minIndex, final int maxIndex, final boolean truncateLast) {

  655.         final int first;
  656.         if (truncateLast) {

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

  659.         } else {

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

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

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

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

  668.             first = minIndex - extraBefore;

  669.         }

  670.         return new Iterable<Range>() {

  671.             /** {@inheritDoc} */
  672.             @Override
  673.             public Iterator<Range> iterator() {
  674.                 return new Iterator<Range>() {

  675.                     private int nextLower = first;

  676.                     /** {@inheritDoc} */
  677.                     @Override
  678.                     public boolean hasNext() {
  679.                         return nextLower < maxIndex;
  680.                     }

  681.                     /** {@inheritDoc} */
  682.                     @Override
  683.                     public Range next() {

  684.                         if (nextLower >= maxIndex) {
  685.                             throw new NoSuchElementException();
  686.                         }
  687.                         final int lower = nextLower;

  688.                         nextLower += quantization;
  689.                         if (truncateLast && nextLower > maxIndex && lower < maxIndex) {
  690.                             // truncate last tile
  691.                             nextLower = maxIndex;
  692.                         }

  693.                         return new Range(lower, nextLower);

  694.                     }

  695.                     /** {@inheritDoc} */
  696.                     @Override
  697.                     public void remove() {
  698.                         throw new UnsupportedOperationException();
  699.                     }

  700.                 };
  701.             }
  702.         };

  703.     }

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

  706.         /** Lower index. */
  707.         private final int lower;

  708.         /** Upper index. */
  709.         private final int upper;

  710.         /** Simple constructor.
  711.          * @param lower lower index
  712.          * @param upper upper index
  713.          */
  714.         Range(final int lower, final int upper) {
  715.             this.lower = lower;
  716.             this.upper = upper;
  717.         }

  718.     }

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

  721.         /** Across range. */
  722.         private final Range across;

  723.         /** Along range. */
  724.         private final Range along;

  725.         /** Simple constructor.
  726.          * @param across across range
  727.          * @param along along range
  728.          */
  729.         RangePair(final Range across, final Range along) {
  730.             this.across = across;
  731.             this.along  = along;
  732.         }

  733.     }

  734. }