BasicScanAlgorithm.java

  1. /* Copyright 2013-2019 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.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.  */
  43. public class BasicScanAlgorithm implements IntersectionAlgorithm {

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

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

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

  54.     /** Simple constructor.
  55.      * @param updater updater used to load Digital Elevation Model tiles
  56.      * @param maxCachedTiles maximum number of tiles stored in the cache
  57.      */
  58.     public BasicScanAlgorithm(final TileUpdater updater, final int maxCachedTiles) {
  59.         cache = new TilesCache<SimpleTile>(new SimpleTileFactory(), updater, maxCachedTiles);
  60.         hMin  = Double.POSITIVE_INFINITY;
  61.         hMax  = Double.NEGATIVE_INFINITY;
  62.     }

  64.     /** {@inheritDoc} */
  65.     @Override
  66.     public NormalizedGeodeticPoint intersection(final ExtendedEllipsoid ellipsoid,
  67.             final Vector3D position, final Vector3D los) {

  69.         DumpManager.dumpAlgorithm(AlgorithmId.BASIC_SLOW_EXHAUSTIVE_SCAN_FOR_TESTS_ONLY);

  71.         // find the tiles between the entry and exit point in the Digital Elevation Model
  72.         NormalizedGeodeticPoint entryPoint = null;
  73.         NormalizedGeodeticPoint exitPoint  = null;
  74.         double minLatitude  = Double.NaN;
  75.         double maxLatitude  = Double.NaN;
  76.         double minLongitude = Double.NaN;
  77.         double maxLongitude = Double.NaN;
  78.         final List<SimpleTile> scannedTiles = new ArrayList<SimpleTile>();
  79.         double centralLongitude = Double.NaN;
  80.         for (boolean changedMinMax = true; changedMinMax; changedMinMax = checkMinMax(scannedTiles)) {

  82.             scannedTiles.clear();
  83.             // compute entry and exit points
  84.             entryPoint = ellipsoid.transform(ellipsoid.pointAtAltitude(position, los, Double.isInfinite(hMax) ? 0.0 : hMax),
  85.                     ellipsoid.getBodyFrame(), null,
  86.                     Double.isNaN(centralLongitude) ? 0.0 : centralLongitude);
  87.             final SimpleTile entryTile = cache.getTile(entryPoint.getLatitude(), entryPoint.getLongitude());
  88.             if (Double.isNaN(centralLongitude)) {
  89.                 centralLongitude = entryTile.getMinimumLongitude();
  90.                 entryPoint = new NormalizedGeodeticPoint(entryPoint.getLatitude(), entryPoint.getLongitude(),
  91.                         entryPoint.getAltitude(), centralLongitude);
  92.             }
  93.             addIfNotPresent(scannedTiles, entryTile);

  95.             exitPoint = ellipsoid.transform(ellipsoid.pointAtAltitude(position, los, Double.isInfinite(hMin) ? 0.0 : hMin),
  96.                     ellipsoid.getBodyFrame(), null, centralLongitude);
  97.             final SimpleTile exitTile = cache.getTile(exitPoint.getLatitude(), exitPoint.getLongitude());
  98.             addIfNotPresent(scannedTiles, exitTile);

  100.             minLatitude  = FastMath.min(entryPoint.getLatitude(),  exitPoint.getLatitude());
  101.             maxLatitude  = FastMath.max(entryPoint.getLatitude(),  exitPoint.getLatitude());
  102.             minLongitude = FastMath.min(entryPoint.getLongitude(), exitPoint.getLongitude());
  103.             maxLongitude = FastMath.max(entryPoint.getLongitude(), exitPoint.getLongitude());

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

  119.         }

  121.         // scan the tiles
  122.         NormalizedGeodeticPoint intersectionGP = null;
  123.         double intersectionDot = Double.POSITIVE_INFINITY;
  124.         for (final SimpleTile tile : scannedTiles) {
  125.             for (int i = latitudeIndex(tile, minLatitude); i <= latitudeIndex(tile, maxLatitude); ++i) {
  126.                 for (int j = longitudeIndex(tile, minLongitude); j <= longitudeIndex(tile, maxLongitude); ++j) {
  127.                     final NormalizedGeodeticPoint gp = tile.cellIntersection(entryPoint, ellipsoid.convertLos(entryPoint, los), i, j);
  128.                     if (gp != null) {

  130.                         // improve the point, by projecting it back on the 3D line, fixing the small body curvature at cell level
  131.                         final Vector3D      delta     = ellipsoid.transform(gp).subtract(position);
  132.                         final double        s         = Vector3D.dotProduct(delta, los) / los.getNormSq();
  133.                         final GeodeticPoint projected = ellipsoid.transform(new Vector3D(1, position, s, los),
  134.                                 ellipsoid.getBodyFrame(), null);
  135.                         final NormalizedGeodeticPoint normalizedProjected = new NormalizedGeodeticPoint(projected.getLatitude(),
  136.                                 projected.getLongitude(),
  137.                                 projected.getAltitude(),
  138.                                 gp.getLongitude());
  139.                         final NormalizedGeodeticPoint gpImproved = tile.cellIntersection(normalizedProjected,
  140.                                 ellipsoid.convertLos(normalizedProjected, los),
  141.                                 i, j);

  143.                         if (gpImproved != null) {
  144.                             final Vector3D point = ellipsoid.transform(gpImproved);
  145.                             final double dot = Vector3D.dotProduct(point.subtract(position), los);
  146.                             if (dot < intersectionDot) {
  147.                                 intersectionGP  = gpImproved;
  148.                                 intersectionDot = dot;
  149.                             }
  150.                         }

  152.                     }
  153.                 }
  154.             }
  155.         }

  157.         return intersectionGP;
  158.     }

  160.     /** {@inheritDoc} */
  161.     @Override
  162.     public NormalizedGeodeticPoint refineIntersection(final ExtendedEllipsoid ellipsoid,
  163.                                                       final Vector3D position, final Vector3D los,
  164.                                                       final NormalizedGeodeticPoint closeGuess) {
  165.         DumpManager.dumpAlgorithm(AlgorithmId.BASIC_SLOW_EXHAUSTIVE_SCAN_FOR_TESTS_ONLY);
  166.         final Vector3D      delta     = ellipsoid.transform(closeGuess).subtract(position);
  167.         final double        s         = Vector3D.dotProduct(delta, los) / los.getNormSq();
  168.         final GeodeticPoint projected = ellipsoid.transform(new Vector3D(1, position, s, los),
  169.                 ellipsoid.getBodyFrame(), null);
  170.         final NormalizedGeodeticPoint normalizedProjected = new NormalizedGeodeticPoint(projected.getLatitude(),
  171.                 projected.getLongitude(),
  172.                 projected.getAltitude(),
  173.                 closeGuess.getLongitude());
  174.         final Tile          tile      = cache.getTile(normalizedProjected.getLatitude(),
  175.                 normalizedProjected.getLongitude());
  176.         return tile.cellIntersection(normalizedProjected,
  177.                 ellipsoid.convertLos(normalizedProjected, los),
  178.                 tile.getFloorLatitudeIndex(normalizedProjected.getLatitude()),
  179.                 tile.getFloorLongitudeIndex(normalizedProjected.getLongitude()));
  180.     }

  182.     /** {@inheritDoc} */
  183.     @Override
  184.     public double getElevation(final double latitude, final double longitude) {
  185.         DumpManager.dumpAlgorithm(AlgorithmId.BASIC_SLOW_EXHAUSTIVE_SCAN_FOR_TESTS_ONLY);
  186.         final Tile tile = cache.getTile(latitude, longitude);
  187.         return tile.interpolateElevation(latitude, longitude);
  188.     }

  190.     /** Check the overall min and max altitudes.
  191.      * @param tiles tiles to check
  192.      * @return true if the tile changed either min or max altitude
  193.      */
  194.     private boolean checkMinMax(final List<SimpleTile> tiles) {

  196.         boolean changedMinMax = false;

  198.         for (final SimpleTile tile : tiles) {

  200.             // check minimum altitude
  201.             if (tile.getMinElevation() < hMin) {
  202.                 hMin          = tile.getMinElevation();
  203.                 changedMinMax = true;
  204.             }

  206.             // check maximum altitude
  207.             if (tile.getMaxElevation() > hMax) {
  208.                 hMax          = tile.getMaxElevation();
  209.                 changedMinMax = true;
  210.             }

  212.         }

  214.         return changedMinMax;

  216.     }

  218.     /** Add a tile to a list if not already present.
  219.      * @param list tiles list
  220.      * @param tile new tile to consider
  221.      */
  222.     private void addIfNotPresent(final List<SimpleTile> list, final SimpleTile tile) {

  224.         // look for existing tiles in the list
  225.         for (final SimpleTile existing : list) {
  226.             if (existing == tile) {
  227.                 return;
  228.             }
  229.         }

  231.         // the tile was not there, add it
  232.         list.add(tile);

  234.     }

  236.     /** Get latitude index.
  237.      * @param tile current tile
  238.      * @param latitude current latitude
  239.      * @return index of latitude, truncated at tiles limits
  240.      */
  241.     private int latitudeIndex(final SimpleTile tile, final double latitude) {
  242.         final int rawIndex = tile.getFloorLatitudeIndex(latitude);
  243.         return FastMath.min(FastMath.max(0, rawIndex), tile.getLatitudeRows());
  244.     }

  246.     /** Get longitude index.
  247.      * @param tile current tile
  248.      * @param longitude current longitude
  249.      * @return index of longitude, truncated at tiles limits
  250.      */
  251.     private int longitudeIndex(final SimpleTile tile, final double longitude) {
  252.         final int rawIndex = tile.getFloorLongitudeIndex(longitude);
  253.         return FastMath.min(FastMath.max(0, rawIndex), tile.getLongitudeColumns());
  254.     }

  256. }
Created with JaCoCo 0.8.2.201808211720