BasicScanAlgorithm.java

  1. /* Copyright 2013-2017 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 org.hipparchus.geometry.euclidean.threed.Vector3D;
  20. import org.hipparchus.util.FastMath;
  21. import java.util.ArrayList;
  22. import java.util.List;

  24. import org.orekit.bodies.GeodeticPoint;
  25. import org.orekit.errors.OrekitException;
  26. import org.orekit.rugged.api.AlgorithmId;
  27. import org.orekit.rugged.errors.DumpManager;
  28. import org.orekit.rugged.errors.RuggedException;
  29. import org.orekit.rugged.raster.SimpleTile;
  30. import org.orekit.rugged.raster.SimpleTileFactory;
  31. import org.orekit.rugged.raster.Tile;
  32. import org.orekit.rugged.raster.TileUpdater;
  33. import org.orekit.rugged.raster.TilesCache;
  34. import org.orekit.rugged.utils.ExtendedEllipsoid;
  35. import org.orekit.rugged.utils.NormalizedGeodeticPoint;

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

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

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

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

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

  66.     /** {@inheritDoc} */
  67.     @Override
  68.     public NormalizedGeodeticPoint intersection(final ExtendedEllipsoid ellipsoid,
  69.                                                 final Vector3D position, final Vector3D los)
  70.         throws RuggedException {
  71.         try {

  73.             DumpManager.dumpAlgorithm(AlgorithmId.BASIC_SLOW_EXHAUSTIVE_SCAN_FOR_TESTS_ONLY);

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

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

  99.                 exitPoint = ellipsoid.transform(ellipsoid.pointAtAltitude(position, los, Double.isInfinite(hMin) ? 0.0 : hMin),
  100.                                                 ellipsoid.getBodyFrame(), null, centralLongitude);
  101.                 final SimpleTile exitTile = cache.getTile(exitPoint.getLatitude(), exitPoint.getLongitude());
  102.                 addIfNotPresent(scannedTiles, exitTile);

  104.                 minLatitude  = FastMath.min(entryPoint.getLatitude(),  exitPoint.getLatitude());
  105.                 maxLatitude  = FastMath.max(entryPoint.getLatitude(),  exitPoint.getLatitude());
  106.                 minLongitude = FastMath.min(entryPoint.getLongitude(), exitPoint.getLongitude());
  107.                 maxLongitude = FastMath.max(entryPoint.getLongitude(), exitPoint.getLongitude());

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

  123.             }

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

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

  147.                             if (gpImproved != null) {
  148.                                 final Vector3D point = ellipsoid.transform(gpImproved);
  149.                                 final double dot = Vector3D.dotProduct(point.subtract(position), los);
  150.                                 if (dot < intersectionDot) {
  151.                                     intersectionGP  = gpImproved;
  152.                                     intersectionDot = dot;
  153.                                 }
  154.                             }

  156.                         }
  157.                     }
  158.                 }
  159.             }

  161.             return intersectionGP;

  163.         } catch (OrekitException oe) {
  164.             // this should never happen
  165.             throw new RuggedException(oe, oe.getSpecifier(), oe.getParts());
  166.         }
  167.     }

  169.     /** {@inheritDoc} */
  170.     @Override
  171.     public NormalizedGeodeticPoint refineIntersection(final ExtendedEllipsoid ellipsoid,
  172.                                                       final Vector3D position, final Vector3D los,
  173.                                                       final NormalizedGeodeticPoint closeGuess)
  174.         throws RuggedException {
  175.         try {
  176.             DumpManager.dumpAlgorithm(AlgorithmId.BASIC_SLOW_EXHAUSTIVE_SCAN_FOR_TESTS_ONLY);
  177.             final Vector3D      delta     = ellipsoid.transform(closeGuess).subtract(position);
  178.             final double        s         = Vector3D.dotProduct(delta, los) / los.getNormSq();
  179.             final GeodeticPoint projected = ellipsoid.transform(new Vector3D(1, position, s, los),
  180.                                                                 ellipsoid.getBodyFrame(), null);
  181.             final NormalizedGeodeticPoint normalizedProjected = new NormalizedGeodeticPoint(projected.getLatitude(),
  182.                                                                                             projected.getLongitude(),
  183.                                                                                             projected.getAltitude(),
  184.                                                                                             closeGuess.getLongitude());
  185.             final Tile          tile      = cache.getTile(normalizedProjected.getLatitude(),
  186.                                                           normalizedProjected.getLongitude());
  187.             return tile.cellIntersection(normalizedProjected,
  188.                                          ellipsoid.convertLos(normalizedProjected, los),
  189.                                          tile.getFloorLatitudeIndex(normalizedProjected.getLatitude()),
  190.                                          tile.getFloorLongitudeIndex(normalizedProjected.getLongitude()));
  191.         } catch (OrekitException oe) {
  192.             throw new RuggedException(oe, oe.getSpecifier(), oe.getParts());
  193.         }
  194.     }

  196.     /** {@inheritDoc} */
  197.     @Override
  198.     public double getElevation(final double latitude, final double longitude)
  199.         throws RuggedException {
  200.         DumpManager.dumpAlgorithm(AlgorithmId.BASIC_SLOW_EXHAUSTIVE_SCAN_FOR_TESTS_ONLY);
  201.         final Tile tile = cache.getTile(latitude, longitude);
  202.         return tile.interpolateElevation(latitude, longitude);
  203.     }

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

  211.         boolean changedMinMax = false;

  213.         for (final SimpleTile tile : tiles) {

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

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

  227.         }

  229.         return changedMinMax;

  231.     }

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

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

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

  249.     }

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

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

  271. }
Created with JaCoCo 0.7.9.201702052155