BasicScanAlgorithm.java

  1. /* Copyright 2013-2025 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.rugged.intersection;

  19. import java.util.ArrayList;
  20. import java.util.List;

  22. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  23. import org.hipparchus.util.FastMath;
  24. import org.orekit.bodies.GeodeticPoint;
  25. import org.orekit.rugged.api.AlgorithmId;
  26. import org.orekit.rugged.errors.DumpManager;
  27. import org.orekit.rugged.raster.SimpleTile;
  28. import org.orekit.rugged.raster.SimpleTileFactory;
  29. import org.orekit.rugged.raster.Tile;
  30. import org.orekit.rugged.raster.TileUpdater;
  31. import org.orekit.rugged.raster.TilesCache;
  32. import org.orekit.rugged.utils.ExtendedEllipsoid;
  33. import org.orekit.rugged.utils.NormalizedGeodeticPoint;

  35. /** Intersection computation using a basic algorithm based on exhaustive scan.
  36.  * <p>
  37.  * The algorithm simply computes entry and exit points at high and low altitudes,
  38.  * and scans all Digital Elevation Models in the sub-tiles defined by these two
  39.  * corner points. It is not designed for operational use.
  40.  * </p>
  41.  * @author Luc Maisonobe
  42.  * @author Guylaine Prat
  43.  */
  44. public class BasicScanAlgorithm implements IntersectionAlgorithm {

  46.     /** Cache for DEM tiles. */
  47.     private final TilesCache<SimpleTile> cache;

  49.     /** Minimum altitude encountered. */
  50.     private double hMin;

  52.     /** Maximum altitude encountered. */
  53.     private double hMax;

  55.     /** Algorithm Id.
  56.      * @since 2.2 */
  57.     private final AlgorithmId algorithmId;

  59.     /** Simple constructor.
  60.      * @param updater updater used to load Digital Elevation Model tiles
  61.      * @param maxCachedTiles maximum number of tiles stored in the cache
  62.      * @param isOverlappingTiles flag to tell if the DEM tiles are overlapping:
  63.      *                          true if overlapping; false otherwise.
  64.      */
  65.     public BasicScanAlgorithm(final TileUpdater updater, final int maxCachedTiles, final boolean isOverlappingTiles) {
  66.         this.cache = new TilesCache<>(new SimpleTileFactory(), updater, maxCachedTiles, isOverlappingTiles);
  67.         this.hMin  = Double.POSITIVE_INFINITY;
  68.         this.hMax  = Double.NEGATIVE_INFINITY;
  69.         this.algorithmId = AlgorithmId.BASIC_SLOW_EXHAUSTIVE_SCAN_FOR_TESTS_ONLY;
  70.     }

  72.     /** {@inheritDoc} */
  73.     @Override
  74.     public NormalizedGeodeticPoint intersection(final ExtendedEllipsoid ellipsoid,
  75.             final Vector3D position, final Vector3D los) {

  77.         DumpManager.dumpAlgorithm(this.algorithmId);

  79.         // find the tiles between the entry and exit point in the Digital Elevation Model
  80.         NormalizedGeodeticPoint entryPoint = null;
  81.         NormalizedGeodeticPoint exitPoint  = null;
  82.         double minLatitude  = Double.NaN;
  83.         double maxLatitude  = Double.NaN;
  84.         double minLongitude = Double.NaN;
  85.         double maxLongitude = Double.NaN;
  86.         final List<SimpleTile> scannedTiles = new ArrayList<>();
  87.         double centralLongitude = Double.NaN;
  88.         for (boolean changedMinMax = true; changedMinMax; changedMinMax = checkMinMax(scannedTiles)) {

  90.             scannedTiles.clear();
  91.             // compute entry and exit points
  92.             entryPoint = ellipsoid.transform(ellipsoid.pointAtAltitude(position, los, Double.isInfinite(hMax) ? 0.0 : hMax),
  93.                     ellipsoid.getBodyFrame(), null,
  94.                     Double.isNaN(centralLongitude) ? 0.0 : centralLongitude);
  95.             final SimpleTile entryTile = cache.getTile(entryPoint.getLatitude(), entryPoint.getLongitude());
  96.             if (Double.isNaN(centralLongitude)) {
  97.                 centralLongitude = entryTile.getMinimumLongitude();
  98.                 entryPoint = new NormalizedGeodeticPoint(entryPoint.getLatitude(), entryPoint.getLongitude(),
  99.                         entryPoint.getAltitude(), centralLongitude);
  100.             }
  101.             addIfNotPresent(scannedTiles, entryTile);

  103.             exitPoint = ellipsoid.transform(ellipsoid.pointAtAltitude(position, los, Double.isInfinite(hMin) ? 0.0 : hMin),
  104.                     ellipsoid.getBodyFrame(), null, centralLongitude);
  105.             final SimpleTile exitTile = cache.getTile(exitPoint.getLatitude(), exitPoint.getLongitude());
  106.             addIfNotPresent(scannedTiles, exitTile);

  108.             minLatitude  = FastMath.min(entryPoint.getLatitude(),  exitPoint.getLatitude());
  109.             maxLatitude  = FastMath.max(entryPoint.getLatitude(),  exitPoint.getLatitude());
  110.             minLongitude = FastMath.min(entryPoint.getLongitude(), exitPoint.getLongitude());
  111.             maxLongitude = FastMath.max(entryPoint.getLongitude(), exitPoint.getLongitude());

  113.             if (scannedTiles.size() > 1) {
  114.                 // the entry and exit tiles are different, maybe other tiles should be added on the way
  115.                 // in the spirit of simple and exhaustive, we add all tiles in a rectangular area
  116.                 final double latStep = 0.5 * FastMath.min(entryTile.getLatitudeStep()  * entryTile.getLatitudeRows(),
  117.                         exitTile.getLatitudeStep()   * exitTile.getLatitudeRows());
  118.                 final double lonStep = 0.5 * FastMath.min(entryTile.getLongitudeStep() * entryTile.getLongitudeColumns(),
  119.                         exitTile.getLongitudeStep()  * exitTile.getLongitudeColumns());
  120.                 for (double latitude = minLatitude; latitude <= maxLatitude; latitude += latStep) {
  121.                     for (double longitude = minLongitude; longitude < maxLongitude; longitude += lonStep) {
  122.                         addIfNotPresent(scannedTiles, cache.getTile(latitude, longitude));
  123.                     }
  124.                 }
  125.             }

  127.         }

  129.         // scan the tiles
  130.         NormalizedGeodeticPoint intersectionGP = null;
  131.         double intersectionDot = Double.POSITIVE_INFINITY;
  132.         for (final SimpleTile tile : scannedTiles) {
  133.             for (int i = latitudeIndex(tile, minLatitude); i <= latitudeIndex(tile, maxLatitude); ++i) {
  134.                 for (int j = longitudeIndex(tile, minLongitude); j <= longitudeIndex(tile, maxLongitude); ++j) {
  135.                     final NormalizedGeodeticPoint gp = tile.cellIntersection(entryPoint, ellipsoid.convertLos(entryPoint, los), i, j);
  136.                     if (gp != null) {

  138.                         // improve the point, by projecting it back on the 3D line, fixing the small body curvature at cell level
  139.                         final Vector3D      delta     = ellipsoid.transform(gp).subtract(position);
  140.                         final double        s         = Vector3D.dotProduct(delta, los) / los.getNormSq();
  141.                         final GeodeticPoint projected = ellipsoid.transform(new Vector3D(1, position, s, los),
  142.                                 ellipsoid.getBodyFrame(), null);
  143.                         final NormalizedGeodeticPoint normalizedProjected = new NormalizedGeodeticPoint(projected.getLatitude(),
  144.                                 projected.getLongitude(),
  145.                                 projected.getAltitude(),
  146.                                 gp.getLongitude());
  147.                         final NormalizedGeodeticPoint gpImproved = tile.cellIntersection(normalizedProjected,
  148.                                 ellipsoid.convertLos(normalizedProjected, los),
  149.                                 i, j);

  151.                         if (gpImproved != null) {
  152.                             final Vector3D point = ellipsoid.transform(gpImproved);
  153.                             final double dot = Vector3D.dotProduct(point.subtract(position), los);
  154.                             if (dot < intersectionDot) {
  155.                                 intersectionGP  = gpImproved;
  156.                                 intersectionDot = dot;
  157.                             }
  158.                         }

  160.                     }
  161.                 }
  162.             }
  163.         }

  165.         return intersectionGP;
  166.     }

  168.     /** {@inheritDoc} */
  169.     @Override
  170.     public NormalizedGeodeticPoint refineIntersection(final ExtendedEllipsoid ellipsoid,
  171.                                                       final Vector3D position, final Vector3D los,
  172.                                                       final NormalizedGeodeticPoint closeGuess) {
  173.         DumpManager.dumpAlgorithm(this.algorithmId);
  174.         final Vector3D      delta     = ellipsoid.transform(closeGuess).subtract(position);
  175.         final double        s         = Vector3D.dotProduct(delta, los) / los.getNormSq();
  176.         final GeodeticPoint projected = ellipsoid.transform(new Vector3D(1, position, s, los),
  177.                 ellipsoid.getBodyFrame(), null);
  178.         final NormalizedGeodeticPoint normalizedProjected = new NormalizedGeodeticPoint(projected.getLatitude(),
  179.                 projected.getLongitude(),
  180.                 projected.getAltitude(),
  181.                 closeGuess.getLongitude());
  182.         final Tile          tile      = cache.getTile(normalizedProjected.getLatitude(),
  183.                 normalizedProjected.getLongitude());
  184.         return tile.cellIntersection(normalizedProjected,
  185.                 ellipsoid.convertLos(normalizedProjected, los),
  186.                 tile.getFloorLatitudeIndex(normalizedProjected.getLatitude()),
  187.                 tile.getFloorLongitudeIndex(normalizedProjected.getLongitude()));
  188.     }

  190.     /** {@inheritDoc} */
  191.     @Override
  192.     public double getElevation(final double latitude, final double longitude) {
  193.         DumpManager.dumpAlgorithm(this.algorithmId);
  194.         final Tile tile = cache.getTile(latitude, longitude);
  195.         return tile.interpolateElevation(latitude, longitude);
  196.     }

  198.     /** {@inheritDoc} */
  199.     @Override
  200.     public AlgorithmId getAlgorithmId() {
  201.         return this.algorithmId;
  202.     }

  204.     /** Check the overall min and max altitudes.
  205.      * @param tiles tiles to check
  206.      * @return true if the tile changed either min or max altitude
  207.      */
  208.     private boolean checkMinMax(final List<SimpleTile> tiles) {

  210.         boolean changedMinMax = false;

  212.         for (final SimpleTile tile : tiles) {

  214.             // check minimum altitude
  215.             if (tile.getMinElevation() < hMin) {
  216.                 hMin          = tile.getMinElevation();
  217.                 changedMinMax = true;
  218.             }

  220.             // check maximum altitude
  221.             if (tile.getMaxElevation() > hMax) {
  222.                 hMax          = tile.getMaxElevation();
  223.                 changedMinMax = true;
  224.             }

  226.         }

  228.         return changedMinMax;

  230.     }

  232.     /** Add a tile to a list if not already present.
  233.      * @param list tiles list
  234.      * @param tile new tile to consider
  235.      */
  236.     private void addIfNotPresent(final List<SimpleTile> list, final SimpleTile tile) {

  238.         // look for existing tiles in the list
  239.         for (final SimpleTile existing : list) {
  240.             if (existing == tile) {
  241.                 return;
  242.             }
  243.         }

  245.         // the tile was not there, add it
  246.         list.add(tile);

  248.     }

  250.     /** Get latitude index.
  251.      * @param tile current tile
  252.      * @param latitude current latitude
  253.      * @return index of latitude, truncated at tiles limits
  254.      */
  255.     private int latitudeIndex(final SimpleTile tile, final double latitude) {
  256.         final int rawIndex = tile.getFloorLatitudeIndex(latitude);
  257.         return FastMath.min(FastMath.max(0, rawIndex), tile.getLatitudeRows());
  258.     }

  260.     /** Get longitude index.
  261.      * @param tile current tile
  262.      * @param longitude current longitude
  263.      * @return index of longitude, truncated at tiles limits
  264.      */
  265.     private int longitudeIndex(final SimpleTile tile, final double longitude) {
  266.         final int rawIndex = tile.getFloorLongitudeIndex(longitude);
  267.         return FastMath.min(FastMath.max(0, rawIndex), tile.getLongitudeColumns());
  268.     }
  269. }
Created with JaCoCo 0.8.12.202403310830