EllipsoidTessellator.java

/* Copyright 2002-2024 CS GROUP
 * Licensed to CS GROUP (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * CS licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.models.earth.tessellation;

import java.util.ArrayList;
import java.util.Collection;
import java.util.IdentityHashMap;
import java.util.Iterator;
import java.util.LinkedList;
import java.util.List;
import java.util.Map;
import java.util.NoSuchElementException;
import java.util.Queue;

import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.geometry.partitioning.BSPTree;
import org.hipparchus.geometry.partitioning.Hyperplane;
import org.hipparchus.geometry.partitioning.RegionFactory;
import org.hipparchus.geometry.partitioning.SubHyperplane;
import org.hipparchus.geometry.spherical.oned.ArcsSet;
import org.hipparchus.geometry.spherical.twod.Circle;
import org.hipparchus.geometry.spherical.twod.S2Point;
import org.hipparchus.geometry.spherical.twod.Sphere2D;
import org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet;
import org.hipparchus.geometry.spherical.twod.SubCircle;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathUtils;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitInternalError;

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

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

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

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

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

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

    /** Tessellate a zone of interest into tiles.
     * <p>
     * The created tiles will completely cover the zone of interest.
     * </p>
     * <p>
     * The distance between a vertex at a tile corner and the vertex at the same corner
     * in the next vertex are computed by subtracting the overlap width (resp. overlap length)
     * from the full width (resp. full length). If for example the full width is specified to
     * be 55 km and the overlap in width is specified to be +5 km, successive tiles would span
     * as follows:
     * </p>
     * <ul>
     *   <li>tile 1 covering from   0 km to  55 km</li>
     *   <li>tile 2 covering from  50 km to 105 km</li>
     *   <li>tile 3 covering from 100 km to 155 km</li>
     *   <li>...</li>
     * </ul>
     * <p>
     * In order to achieve the same 50 km step but using a 5 km gap instead of an overlap, one would
     * need to specify the full width to be 45 km and the overlap to be -5 km. With these settings,
     * successive tiles would span as follows:
     * </p>
     * <ul>
     *   <li>tile 1 covering from   0 km to  45 km</li>
     *   <li>tile 2 covering from  50 km to  95 km</li>
     *   <li>tile 3 covering from 100 km to 155 km</li>
     *   <li>...</li>
     * </ul>
     * @param zone zone of interest to tessellate
     * @param fullWidth full tiles width as a distance on surface, including overlap (in meters)
     * @param fullLength full tiles length as a distance on surface, including overlap (in meters)
     * @param widthOverlap overlap between adjacent tiles (in meters), if negative the tiles
     * will have a gap between each other instead of an overlap
     * @param lengthOverlap overlap between adjacent tiles (in meters), if negative the tiles
     * will have a gap between each other instead of an overlap
     * @param truncateLastWidth if true, the first tiles strip will be started as close as
     * possible to the zone of interest, and the last tiles strip will have its width reduced
     * to also remain close to the zone of interest; if false all tiles strip will have the
     * same {@code fullWidth} and they will be balanced around zone of interest
     * @param truncateLastLength if true, the first tile in each strip will be started as close as
     * possible to the zone of interest, and the last tile in each strip will have its length reduced
     * to also remain close to the zone of interest; if false all tiles in each strip will have the
     * same {@code fullLength} and they will be balanced around zone of interest
     * @return a list of lists of tiles covering the zone of interest,
     * each sub-list corresponding to a part not connected to the other
     * parts (for example for islands)
     */
    public List<List<Tile>> tessellate(final SphericalPolygonsSet zone,
                                       final double fullWidth, final double fullLength,
                                       final double widthOverlap, final double lengthOverlap,
                                       final boolean truncateLastWidth, final boolean truncateLastLength) {

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

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

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

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

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

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

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

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

                    }
                }

            }

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

            map.put(mesh, tiles);

        }

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

        return tilesLists;

    }

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

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

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

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

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

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

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

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

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

                    }
                }

            }

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

            map.put(mesh, sample);

        }

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

        return sampleLists;

    }

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

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

    }

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

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

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

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

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

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

            }

            // second expansion step: check if the loop of outside nodes
            // completely surrounds the zone, i.e. there are no peaks
            // pointing out of the loop between two nodes
            expanding = false;
            final List<Mesh.Node> boundary = mesh.getTaxicabBoundary(false);
            if (boundary.size() > 1) {
                Mesh.Node previous = boundary.get(boundary.size() - 1);
                for (final Mesh.Node node : boundary) {
                    if (meetInside(previous.getS2P(), node.getS2P(), zone)) {
                        // part of the mesh boundary is still inside the zone!
                        // the mesh must be expanded again
                        addAllNeighborsIfNeeded(previous, mesh, newNodes);
                        addAllNeighborsIfNeeded(node,     mesh, newNodes);
                        expanding = true;
                    }
                    previous = node;
                }
            }

        }

    }

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

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

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

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

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

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

                // apply tile overlap
                final S2Point s2p0 = node0.move(new Vector3D(-0.5 * lengthOverlap, node0.getAlong(),
                                                             -0.5 * widthOverlap,  node0.getAcross()));
                final S2Point s2p1 = node1.move(new Vector3D(+0.5 * lengthOverlap, node1.getAlong(),
                                                             -0.5 * widthOverlap,  node1.getAcross()));
                final S2Point s2p2 = node2.move(new Vector3D(+0.5 * lengthOverlap, node2.getAlong(),
                                                             +0.5 * widthOverlap,  node2.getAcross()));
                final S2Point s2p3 = node3.move(new Vector3D(-0.5 * lengthOverlap, node2.getAlong(),
                                                             +0.5 * widthOverlap,  node2.getAcross()));

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

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

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

                }

            }
        }

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

        return tiles;

    }

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

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

                // count how many points would be selected for the current modulus
                int count = 0;
                for (int across = mesh.getMinAcrossIndex() + acrossModulus;
                        across <= mesh.getMaxAcrossIndex();
                        across += quantization) {
                    for (int along = mesh.getMinAlongIndex() + alongModulus;
                            along <= mesh.getMaxAlongIndex();
                            along += quantization) {
                        final Mesh.Node  node = mesh.getNode(along, across);
                        if (node != null && node.isInside()) {
                            ++count;
                        }
                    }
                }

                if (count > selectedCount) {
                    // current modulus are better than the selected ones
                    selectedAcrossModulus = acrossModulus;
                    selectedAlongModulus  = alongModulus;
                    selectedCount         = count;
                }
            }
        }

        // extract the sample points
        final List<GeodeticPoint> sample = new ArrayList<GeodeticPoint>(selectedCount);
        for (int across = mesh.getMinAcrossIndex() + selectedAcrossModulus;
                across <= mesh.getMaxAcrossIndex();
                across += quantization) {
            for (int along = mesh.getMinAlongIndex() + selectedAlongModulus;
                    along <= mesh.getMaxAlongIndex();
                    along += quantization) {
                final Mesh.Node  node = mesh.getNode(along, across);
                if (node != null && node.isInside()) {
                    sample.add(toGeodetic(node.getS2P()));
                }
            }
        }

        return sample;

    }

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

        // select the way merge will be performed
        final Mesh larger;
        final Mesh smaller;
        if (mesh1.getNumberOfNodes() >= mesh2.getNumberOfNodes()) {
            // the larger new mesh should absorb the smaller existing mesh
            larger  = mesh1;
            smaller = mesh2;
        } else {
            // the larger existing mesh should absorb the smaller new mesh
            larger  = mesh2;
            smaller = mesh1;
        }

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

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

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

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

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

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

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

        }

        return larger;

    }

    /** Ensure all 8 neighbors of a node are in the mesh.
     * @param base base node
     * @param mesh complete mesh containing nodes
     * @param newNodes queue where new node must be put
     */
    private void addAllNeighborsIfNeeded(final Mesh.Node base, final Mesh mesh,
                                         final Collection<Mesh.Node> newNodes) {
        addNode(base.getAlongIndex() - 1, base.getAcrossIndex() - 1, mesh, newNodes);
        addNode(base.getAlongIndex() - 1, base.getAcrossIndex(),     mesh, newNodes);
        addNode(base.getAlongIndex() - 1, base.getAcrossIndex() + 1, mesh, newNodes);
        addNode(base.getAlongIndex(),     base.getAcrossIndex() - 1, mesh, newNodes);
        addNode(base.getAlongIndex(),     base.getAcrossIndex() + 1, mesh, newNodes);
        addNode(base.getAlongIndex() + 1, base.getAcrossIndex() - 1, mesh, newNodes);
        addNode(base.getAlongIndex() + 1, base.getAcrossIndex(),     mesh, newNodes);
        addNode(base.getAlongIndex() + 1, base.getAcrossIndex() + 1, mesh, newNodes);
    }

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

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

        if (!node.isEnabled()) {
            // enable the node
            node.setEnabled();
            newNodes.add(node);
        }

    }

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

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

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

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

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

    /** Check if an arc meets the inside of a zone.
     * @param s1 first point
     * @param s2 second point
     * @param zone zone to check arc against
     * @return true if the arc meets the inside of the zone
     */
    private boolean meetInside(final S2Point s1, final S2Point s2,
                               final SphericalPolygonsSet zone) {
        final Circle  circle = new Circle(s1, s2, zone.getTolerance());
        final double alpha1  = circle.toSubSpace(s1).getAlpha();
        final double alpha2  = MathUtils.normalizeAngle(circle.toSubSpace(s2).getAlpha(),
                                                        alpha1 + FastMath.PI);
        final SubCircle sub  = new SubCircle(circle,
                                             new ArcsSet(alpha1, alpha2, zone.getTolerance()));
        return recurseMeetInside(zone.getTree(false), sub);

    }

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

        if (node.getCut() == null) {
            // we have reached a leaf node
            if (sub.isEmpty()) {
                return false;
            } else {
                return (Boolean) node.getAttribute();
            }
        } else {
            final Hyperplane<Sphere2D> hyperplane = node.getCut().getHyperplane();
            final SubHyperplane.SplitSubHyperplane<Sphere2D> split = sub.split(hyperplane);
            switch (split.getSide()) {
                case PLUS:
                    return recurseMeetInside(node.getPlus(),  sub);
                case MINUS:
                    return recurseMeetInside(node.getMinus(), sub);
                case BOTH:
                    if (recurseMeetInside(node.getPlus(), split.getPlus())) {
                        return true;
                    } else {
                        return recurseMeetInside(node.getMinus(), split.getMinus());
                    }
                default:
                    // this should not happen
                    throw new OrekitInternalError(null);
            }
        }
    }

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

        final int first;
        if (truncateLast) {

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

        } else {

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

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

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

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

            first = minIndex - extraBefore;

        }

        return new Iterable<Range>() {

            /** {@inheritDoc} */
            @Override
            public Iterator<Range> iterator() {
                return new Iterator<Range>() {

                    private int nextLower = first;

                    /** {@inheritDoc} */
                    @Override
                    public boolean hasNext() {
                        return nextLower < maxIndex;
                    }

                    /** {@inheritDoc} */
                    @Override
                    public Range next() {

                        if (nextLower >= maxIndex) {
                            throw new NoSuchElementException();
                        }
                        final int lower = nextLower;

                        nextLower += quantization;
                        if (truncateLast && nextLower > maxIndex && lower < maxIndex) {
                            // truncate last tile
                            nextLower = maxIndex;
                        }

                        return new Range(lower, nextLower);

                    }

                    /** {@inheritDoc} */
                    @Override
                    public void remove() {
                        throw new UnsupportedOperationException();
                    }

                };
            }
        };

    }

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

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

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

        /** Simple constructor.
         * @param lower lower index
         * @param upper upper index
         */
        Range(final int lower, final int upper) {
            this.lower = lower;
            this.upper = upper;
        }

    }

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

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

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

        /** Simple constructor.
         * @param across across range
         * @param along along range
         */
        RangePair(final Range across, final Range along) {
            this.across = across;
            this.along  = along;
        }

    }

}