SimpleTile.java

  1. /* Copyright 2013-2022 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.raster;

  19. import java.util.Arrays;

  21. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  22. import org.hipparchus.util.FastMath;
  23. import org.hipparchus.util.Precision;
  24. import org.orekit.rugged.errors.DumpManager;
  25. import org.orekit.rugged.errors.RuggedException;
  26. import org.orekit.rugged.errors.RuggedMessages;
  27. import org.orekit.rugged.utils.MaxSelector;
  28. import org.orekit.rugged.utils.MinSelector;
  29. import org.orekit.rugged.utils.NormalizedGeodeticPoint;


  32. /** Simple implementation of a {@link Tile}.
  33.  * @see SimpleTileFactory
  34.  * @author Luc Maisonobe
  35.  * @author Guylaine Prat
  36.  */
  37. public class SimpleTile implements Tile {

  39.     /** Tolerance used to interpolate points slightly out of tile (in cells). */
  40.     private static final double TOLERANCE = 1.0 / 8.0;

  42.     /** Minimum latitude. */
  43.     private double minLatitude;

  45.     /** Minimum longitude. */
  46.     private double minLongitude;

  48.     /** Step in latitude (size of one raster element). */
  49.     private double latitudeStep;

  51.     /** Step in longitude (size of one raster element). */
  52.     private double longitudeStep;

  54.     /** Number of latitude rows. */
  55.     private int latitudeRows;

  57.     /** Number of longitude columns. */
  58.     private int longitudeColumns;

  60.     /** Minimum elevation. */
  61.     private double minElevation;

  63.     /** Latitude index of min elevation. */
  64.     private int minElevationLatitudeIndex;

  66.     /** Longitude index of min elevation. */
  67.     private int minElevationLongitudeIndex;

  69.     /** Maximum elevation. */
  70.     private double maxElevation;

  72.     /** Latitude index of max elevation. */
  73.     private int maxElevationLatitudeIndex;

  75.     /** Longitude index of max elevation. */
  76.     private int maxElevationLongitudeIndex;

  78.     /** Elevation array. */
  79.     private double[] elevations;

  81.     /** Simple constructor.
  82.      * <p>
  83.      * Creates an empty tile.
  84.      * </p>
  85.      */
  86.     protected SimpleTile() {
  87.     }

  89.     /** {@inheritDoc} */
  90.     @Override
  91.     public void setGeometry(final double newMinLatitude, final double newMinLongitude,
  92.                             final double newLatitudeStep, final double newLongitudeStep,
  93.                             final int newLatitudeRows, final int newLongitudeColumns) {
  94.         this.minLatitude                = newMinLatitude;
  95.         this.minLongitude               = newMinLongitude;
  96.         this.latitudeStep               = newLatitudeStep;
  97.         this.longitudeStep              = newLongitudeStep;
  98.         this.latitudeRows               = newLatitudeRows;
  99.         this.longitudeColumns           = newLongitudeColumns;
  100.         this.minElevation               = Double.POSITIVE_INFINITY;
  101.         this.minElevationLatitudeIndex  = -1;
  102.         this.minElevationLongitudeIndex = -1;
  103.         this.maxElevation               = Double.NEGATIVE_INFINITY;
  104.         this.maxElevationLatitudeIndex  = -1;
  105.         this.maxElevationLongitudeIndex = -1;

  107.         if (newLatitudeRows < 1 || newLongitudeColumns < 1) {
  108.             throw new RuggedException(RuggedMessages.EMPTY_TILE, newLatitudeRows, newLongitudeColumns);
  109.         }
  110.         this.elevations = new double[newLatitudeRows * newLongitudeColumns];
  111.         Arrays.fill(elevations, Double.NaN);

  113.     }

  115.     /** {@inheritDoc} */
  116.     @Override
  117.     public void tileUpdateCompleted() {
  118.         processUpdatedElevation(elevations);
  119.     }

  121.     /** Process elevation array at completion.
  122.      * <p>
  123.      * This method is called at tile update completion, it is
  124.      * expected to be overridden by subclasses. The default
  125.      * implementation does nothing.
  126.      * </p>
  127.      * @param elevationsArray elevations array
  128.      */
  129.     protected void processUpdatedElevation(final double[] elevationsArray) {
  130.         // do nothing by default
  131.     }

  133.     /** {@inheritDoc} */
  134.     @Override
  135.     public double getMinimumLatitude() {
  136.         return getLatitudeAtIndex(0);
  137.     }

  139.     /** {@inheritDoc} */
  140.     @Override
  141.     public double getLatitudeAtIndex(final int latitudeIndex) {
  142.         return minLatitude + latitudeStep * latitudeIndex;
  143.     }

  145.     /** {@inheritDoc} */
  146.     @Override
  147.     public double getMaximumLatitude() {
  148.         return getLatitudeAtIndex(latitudeRows - 1);
  149.     }

  151.     /** {@inheritDoc} */
  152.     @Override
  153.     public double getMinimumLongitude() {
  154.         return getLongitudeAtIndex(0);
  155.     }

  157.     /** {@inheritDoc} */
  158.     @Override
  159.     public double getLongitudeAtIndex(final int longitudeIndex) {
  160.         return minLongitude + longitudeStep * longitudeIndex;
  161.     }

  163.     /** {@inheritDoc} */
  164.     @Override
  165.     public double getMaximumLongitude() {
  166.         return getLongitudeAtIndex(longitudeColumns - 1);
  167.     }

  169.     /** {@inheritDoc} */
  170.     @Override
  171.     public double getLatitudeStep() {
  172.         return latitudeStep;
  173.     }

  175.     /** {@inheritDoc} */
  176.     @Override
  177.     public double getLongitudeStep() {
  178.         return longitudeStep;
  179.     }

  181.     /** {@inheritDoc} */
  182.     @Override
  183.     public int getLatitudeRows() {
  184.         return latitudeRows;
  185.     }

  187.     /** {@inheritDoc} */
  188.     @Override
  189.     public int getLongitudeColumns() {
  190.         return longitudeColumns;
  191.     }

  193.     /** {@inheritDoc} */
  194.     @Override
  195.     public double getMinElevation() {
  196.         return minElevation;
  197.     }

  199.     /** {@inheritDoc} */
  200.     @Override
  201.     public int getMinElevationLatitudeIndex() {
  202.         return minElevationLatitudeIndex;
  203.     }

  205.     /** {@inheritDoc} */
  206.     @Override
  207.     public int getMinElevationLongitudeIndex() {
  208.         return minElevationLongitudeIndex;
  209.     }

  211.     /** {@inheritDoc} */
  212.     @Override
  213.     public double getMaxElevation() {
  214.         return maxElevation;
  215.     }

  217.     /** {@inheritDoc} */
  218.     @Override
  219.     public int getMaxElevationLatitudeIndex() {
  220.         return maxElevationLatitudeIndex;
  221.     }

  223.     /** {@inheritDoc} */
  224.     @Override
  225.     public int getMaxElevationLongitudeIndex() {
  226.         return maxElevationLongitudeIndex;
  227.     }

  229.     /** {@inheritDoc} */
  230.     @Override
  231.     public void setElevation(final int latitudeIndex, final int longitudeIndex, final double elevation) {

  233.         if (latitudeIndex  < 0 || latitudeIndex  > (latitudeRows - 1) ||
  234.             longitudeIndex < 0 || longitudeIndex > (longitudeColumns - 1)) {
  235.             throw new RuggedException(RuggedMessages.OUT_OF_TILE_INDICES,
  236.                                       latitudeIndex, longitudeIndex,
  237.                                       latitudeRows - 1, longitudeColumns - 1);
  238.         }
  239.         if (MinSelector.getInstance().selectFirst(elevation, minElevation)) {
  240.             minElevation               = elevation;
  241.             minElevationLatitudeIndex  = latitudeIndex;
  242.             minElevationLongitudeIndex = longitudeIndex;
  243.         }
  244.         if (MaxSelector.getInstance().selectFirst(elevation, maxElevation)) {
  245.             maxElevation               = elevation;
  246.             maxElevationLatitudeIndex  = latitudeIndex;
  247.             maxElevationLongitudeIndex = longitudeIndex;
  248.         }
  249.         elevations[latitudeIndex * getLongitudeColumns() + longitudeIndex] = elevation;
  250.     }

  252.     /** {@inheritDoc} */
  253.     @Override
  254.     public double getElevationAtIndices(final int latitudeIndex, final int longitudeIndex) {
  255.         final double elevation = elevations[latitudeIndex * getLongitudeColumns() + longitudeIndex];
  256.         DumpManager.dumpTileCell(this, latitudeIndex, longitudeIndex, elevation);
  257.         return elevation;
  258.     }

  260.     /** {@inheritDoc}
  261.      * <p>
  262.      * This classes uses an arbitrary 1/8 cell tolerance for interpolating
  263.      * slightly out of tile points.
  264.      * </p>
  265.      */
  266.     @Override
  267.     public double interpolateElevation(final double latitude, final double longitude) {

  269.         final double doubleLatitudeIndex  = getDoubleLatitudeIndex(latitude);
  270.         final double doubleLongitudeIndex = getDoubleLontitudeIndex(longitude);
  271.         if (doubleLatitudeIndex  < -TOLERANCE || doubleLatitudeIndex  >= (latitudeRows - 1 + TOLERANCE) ||
  272.             doubleLongitudeIndex < -TOLERANCE || doubleLongitudeIndex >= (longitudeColumns - 1 + TOLERANCE)) {
  273.             throw new RuggedException(RuggedMessages.OUT_OF_TILE_ANGLES,
  274.                                       FastMath.toDegrees(latitude),
  275.                                       FastMath.toDegrees(longitude),
  276.                                       FastMath.toDegrees(getMinimumLatitude()),
  277.                                       FastMath.toDegrees(getMaximumLatitude()),
  278.                                       FastMath.toDegrees(getMinimumLongitude()),
  279.                                       FastMath.toDegrees(getMaximumLongitude()));
  280.         }

  282.         final int latitudeIndex  = FastMath.max(0,
  283.                                                 FastMath.min(latitudeRows - 2,
  284.                                                              (int) FastMath.floor(doubleLatitudeIndex)));
  285.         final int longitudeIndex = FastMath.max(0,
  286.                                                 FastMath.min(longitudeColumns - 2,
  287.                                                              (int) FastMath.floor(doubleLongitudeIndex)));

  289.         // bi-linear interpolation
  290.         final double dLat = doubleLatitudeIndex  - latitudeIndex;
  291.         final double dLon = doubleLongitudeIndex - longitudeIndex;
  292.         final double e00  = getElevationAtIndices(latitudeIndex,     longitudeIndex);
  293.         final double e10  = getElevationAtIndices(latitudeIndex,     longitudeIndex + 1);
  294.         final double e01  = getElevationAtIndices(latitudeIndex + 1, longitudeIndex);
  295.         final double e11  = getElevationAtIndices(latitudeIndex + 1, longitudeIndex + 1);

  297.         return (e00 * (1.0 - dLon) + dLon * e10) * (1.0 - dLat) +
  298.                (e01 * (1.0 - dLon) + dLon * e11) * dLat;

  300.     }

  302.     /** {@inheritDoc} */
  303.     @Override
  304.     public NormalizedGeodeticPoint cellIntersection(final NormalizedGeodeticPoint p, final Vector3D los,
  305.                                                     final int latitudeIndex, final int longitudeIndex) {

  307.         // ensure neighboring cells to not fall out of tile
  308.         final int iLat  = FastMath.max(0, FastMath.min(latitudeRows     - 2, latitudeIndex));
  309.         final int jLong = FastMath.max(0, FastMath.min(longitudeColumns - 2, longitudeIndex));

  311.         // Digital Elevation Mode coordinates at cell vertices
  312.         final double x00 = getLongitudeAtIndex(jLong);
  313.         final double y00 = getLatitudeAtIndex(iLat);
  314.         final double z00 = getElevationAtIndices(iLat,     jLong);
  315.         final double z01 = getElevationAtIndices(iLat + 1, jLong);
  316.         final double z10 = getElevationAtIndices(iLat,     jLong + 1);
  317.         final double z11 = getElevationAtIndices(iLat + 1, jLong + 1);

  319.         // normalize back to tile coordinates
  320.         final NormalizedGeodeticPoint tileP = new NormalizedGeodeticPoint(p.getLatitude(),
  321.                                                                           p.getLongitude(),
  322.                                                                           p.getAltitude(),
  323.                                                                           x00);

  325.         // line-of-sight coordinates at close points
  326.         final double dxA = (tileP.getLongitude() - x00) / longitudeStep;
  327.         final double dyA = (tileP.getLatitude()  - y00) / latitudeStep;
  328.         final double dzA = tileP.getAltitude();
  329.         final double dxB = dxA + los.getX() / longitudeStep;
  330.         final double dyB = dyA + los.getY() / latitudeStep;
  331.         final double dzB = dzA + los.getZ();

  333.         // points along line-of-sight can be defined as a linear progression
  334.         // along the line depending on free variable t: p(t) = p + t * los.
  335.         // As the point latitude and longitude are linear with respect to t,
  336.         // and as Digital Elevation Model is quadratic with respect to latitude
  337.         // and longitude, the altitude of DEM at same horizontal position as
  338.         // point is quadratic in t:
  339.         // z_DEM(t) = u t² + v t + w
  340.         final double u = (dxA - dxB) * (dyA - dyB) * (z00 - z10 - z01 + z11);
  341.         final double v = ((dxA - dxB) * (1 - dyA) + (dyA - dyB) * (1 - dxA)) * z00 +
  342.                          (dxA * (dyA - dyB) - (dxA - dxB) * (1 - dyA)) * z10 +
  343.                          (dyA * (dxA - dxB) - (dyA - dyB) * (1 - dxA)) * z01 +
  344.                          ((dxB - dxA) * dyA + (dyB - dyA) * dxA) * z11;
  345.         final double w = (1 - dxA) * ((1 - dyA) * z00 + dyA * z01) +
  346.                          dxA       * ((1 - dyA) * z10 + dyA * z11);

  348.         // subtract linear z from line-of-sight
  349.         // z_DEM(t) - z_LOS(t) = a t² + b t + c
  350.         final double a = u;
  351.         final double b = v + dzA - dzB;
  352.         final double c = w - dzA;

  354.         // solve the equation
  355.         final double t1;
  356.         final double t2;
  357.         if (FastMath.abs(a) <= Precision.EPSILON * FastMath.abs(c)) {
  358.             // the equation degenerates to a linear (or constant) equation
  359.             final double t = -c / b;
  360.             t1 = Double.isNaN(t) ? 0.0 : t;
  361.             t2 = Double.POSITIVE_INFINITY;
  362.         } else {
  363.             // the equation is quadratic
  364.             final double b2  = b * b;
  365.             final double fac = 4 * a * c;
  366.             if (b2 < fac) {
  367.                 // no intersection at all
  368.                 return null;
  369.             }
  370.             final double s = FastMath.sqrt(b2 - fac);
  371.             t1 = (b < 0) ? (s - b) / (2 * a) : -2 * c / (b + s);
  372.             t2 = c / (a * t1);

  374.         }

  376.         final NormalizedGeodeticPoint p1 = interpolate(t1, tileP, dxA, dyA, los, x00);
  377.         final NormalizedGeodeticPoint p2 = interpolate(t2, tileP, dxA, dyA, los, x00);

  379.         // select the first point along line-of-sight
  380.         if (p1 == null) {
  381.             return p2;
  382.         } else if (p2 == null) {
  383.             return p1;
  384.         } else {
  385.             return t1 <= t2 ? p1 : p2;
  386.         }

  388.     }

  390.     /** Interpolate point along a line.
  391.      * @param t abscissa along the line
  392.      * @param tileP start point, normalized to tile area
  393.      * @param dxP relative coordinate of the start point with respect to current cell
  394.      * @param dyP relative coordinate of the start point with respect to current cell
  395.      * @param los direction of the line-of-sight, in geodetic space
  396.      * @param centralLongitude reference longitude lc such that the point longitude will
  397.      * be normalized between lc-π and lc+π
  398.      * @return interpolated point along the line
  399.      */
  400.     private NormalizedGeodeticPoint interpolate(final double t, final NormalizedGeodeticPoint tileP,
  401.                                                 final double dxP, final double dyP,
  402.                                                 final Vector3D los, final double centralLongitude) {

  404.         if (Double.isInfinite(t)) {
  405.             return null;
  406.         }

  408.         final double dx = dxP + t * los.getX() / longitudeStep;
  409.         final double dy = dyP + t * los.getY() / latitudeStep;
  410.         if (dx >= -TOLERANCE && dx <= 1 + TOLERANCE && dy >= -TOLERANCE && dy <= 1 + TOLERANCE) {
  411.             return new NormalizedGeodeticPoint(tileP.getLatitude()  + t * los.getY(),
  412.                                                tileP.getLongitude() + t * los.getX(),
  413.                                                tileP.getAltitude()  + t * los.getZ(),
  414.                                                centralLongitude);
  415.         } else {
  416.             return null;
  417.         }

  419.     }

  421.     /** {@inheritDoc} */
  422.     @Override
  423.     public int getFloorLatitudeIndex(final double latitude) {
  424.         return (int) FastMath.floor(getDoubleLatitudeIndex(latitude));
  425.     }

  427.     /** {@inheritDoc} */
  428.     @Override
  429.     public int getFloorLongitudeIndex(final double longitude) {
  430.         return (int) FastMath.floor(getDoubleLontitudeIndex(longitude));
  431.     }

  433.     /** Get the latitude index of a point.
  434.      * @param latitude geodetic latitude (rad)
  435.      * @return latitude index (it may lie outside of the tile!)
  436.      */
  437.     private double getDoubleLatitudeIndex(final double latitude) {
  438.         return (latitude  - minLatitude)  / latitudeStep;
  439.     }

  441.     /** Get the longitude index of a point.
  442.      * @param longitude geodetic longitude (rad)
  443.      * @return longitude index (it may lie outside of the tile!)
  444.      */
  445.     private double getDoubleLontitudeIndex(final double longitude) {
  446.         return (longitude - minLongitude) / longitudeStep;
  447.     }

  449.     /** {@inheritDoc} */
  450.     @Override
  451.     public Location getLocation(final double latitude, final double longitude) {
  452.         final int latitudeIndex  = getFloorLatitudeIndex(latitude);
  453.         final int longitudeIndex = getFloorLongitudeIndex(longitude);
  454.         if (longitudeIndex < 0) {
  455.             if (latitudeIndex < 0) {
  456.                 return Location.SOUTH_WEST;
  457.             } else if (latitudeIndex <= (latitudeRows - 2)) {
  458.                 return Location.WEST;
  459.             } else {
  460.                 return Location.NORTH_WEST;
  461.             }
  462.         } else if (longitudeIndex <= (longitudeColumns - 2)) {
  463.             if (latitudeIndex < 0) {
  464.                 return Location.SOUTH;
  465.             } else if (latitudeIndex <= (latitudeRows - 2)) {
  466.                 return Location.HAS_INTERPOLATION_NEIGHBORS;
  467.             } else {
  468.                 return Location.NORTH;
  469.             }
  470.         } else {
  471.             if (latitudeIndex < 0) {
  472.                 return Location.SOUTH_EAST;
  473.             } else if (latitudeIndex <= (latitudeRows - 2)) {
  474.                 return Location.EAST;
  475.             } else {
  476.                 return Location.NORTH_EAST;
  477.             }
  478.         }
  479.     }

  481. }
Created with JaCoCo 0.8.7.202105040129