Fix various bugs with raster data and tiling

src/main/java/fr/ign/voxatile/core/inputs/FloatArrayGeographicDataMatrix2d.java

@@ -38,6 +38,8 @@ public FloatArrayGeographicDataMatrix2d(int sizeX, int sizeY, double offsetX, do
 
     @Override
     public float getFloat(int x, int y) {
+        if (x < 0 || x >= sizeX || y < 0 || y >= sizeY)
+            throw new IndexOutOfBoundsException();
         return data[x + (sizeY - y - 1) * sizeX];
     }
 }

src/main/java/fr/ign/voxatile/core/inputs/GeoTiffDataProvider.java

@@ -9,6 +9,7 @@
 import org.geotools.api.referencing.crs.CoordinateReferenceSystem;
 import org.geotools.coverage.grid.GridCoverage2D;
 import org.geotools.coverage.grid.GridEnvelope2D;
+import org.geotools.coverage.grid.InvalidGridGeometryException;
 import org.geotools.gce.geotiff.GeoTiffReader;
 import org.geotools.geometry.jts.ReferencedEnvelope;
 import org.geotools.util.factory.Hints;
@@ -58,30 +59,56 @@ public Provider.Result<FloatGeographicDataMatrix2d> provide(WorldBBox3d bbox) th
         } catch (IOException e) {
             throw new RetryableException(e);
         }
+
         CoordinateReferenceSystem crs = grid.getCoordinateReferenceSystem();
 
         ReferencedEnvelope envelope;
-        GridEnvelope2D gridEnvelope;
         try {
             envelope = envelopeProvider.computeForCRS(crs, bbox).intersection(grid.getGridGeometry().getEnvelope2D());
+        } catch (FactoryException | TransformException e) {
+            throw new GenerationFailedException(e);
+        }
+
+        // Compute in grid envelope (pixel envelope)
+        GridEnvelope2D gridEnvelope;
+        try {
             gridEnvelope = grid.getGridGeometry().worldToGrid(envelope);
-        } catch (FactoryException | TransformException | org.geotools.api.referencing.operation.TransformException e) {
+        } catch (InvalidGridGeometryException | org.geotools.api.referencing.operation.TransformException e) {
+           throw new GenerationFailedException(e);
+        }
+
+        // This will ensure we have enough pixels for interpolation
+        // - First, we need to be sure we always have four samples around each voxel center.
+        //   I.E: we need to be sure we always include outer samples
+        //   gridEnvelope only includes a sample if original envelope overlaps its pixel.
+        //   this means outer sample is only included if envelope limit is closer to them than to inner sample.
+        // - Then, we have to fix envelope computation which excludes last line/column.
+        // We could do that more accurately. Growing by 2 pixel in every direction will work but may include one or two useless pixel columns or lines.
+        // TODO: grown envelope may overflow out of raster data. Is this correct?
+        gridEnvelope.grow(2, 2);
+
+        // Now we reconvert actual grid envelope to map coordinates
+        // This will be used to consistently adjust pixel positions to map coordinates
+        ReferencedEnvelope mapEnvelope;
+
+        try {
+            mapEnvelope = grid.getGridGeometry().gridToWorld(gridEnvelope);
+        } catch (org.geotools.api.referencing.operation.TransformException e) {
             throw new GenerationFailedException(e);
         }
 
         // RandomIter provides a view of the underlying image to read arbitrary pixel values
         RandomIter data = RandomIterFactory.create(grid.getRenderedImage(), gridEnvelope);
-
         FloatGeographicDataMatrix2d result = new FloatImageGeographicDataMatrix2d(
             data,
             gridEnvelope.x,
             gridEnvelope.y,
             gridEnvelope.width,
             gridEnvelope.height,
-            envelope.getMinX(),
-            envelope.getMinY(),
-            envelope.getWidth() / gridEnvelope.width,
-            envelope.getHeight() / gridEnvelope.height
+            mapEnvelope.getMinX(),
+            mapEnvelope.getMinY(),
+            mapEnvelope.getWidth() / gridEnvelope.width,
+            mapEnvelope.getHeight() / gridEnvelope.height
         );
 
         return new SimpleResult<>(crs, Iterators.iterator(result));

src/main/java/fr/ign/voxatile/core/inputs/WMSFloatBilDataProvider.java

@@ -31,13 +31,18 @@ public class WMSFloatBilDataProvider implements Provider<FloatGeographicDataMatr
     private final EnvelopeProvider envelopeProvider;
     private final String srsName;
 
+    // Voxel size in map coordinates
+    private final double voxelSizeX;
+    private final double voxelSizeY;
+
     /**
      * Creates a new {@code WMSDataProvider}.
      *
      * @param baseURL base URL of the service
      * @param layer name of the WMS layer to query
      * @param crs coordinate reference system to use for this source
      * @param envelopeProvider function to use to compute envelopes from bounding boxes
+     * @throws GenerationFailedException if voxel size could not be computed in layer CRS
      */
     public WMSFloatBilDataProvider(String baseURL, String layer, CoordinateReferenceSystem crs, EnvelopeProvider envelopeProvider) {
         this.crs = crs;
@@ -55,6 +60,18 @@ public WMSFloatBilDataProvider(String baseURL, String layer, CoordinateReference
             .parameter("FORMAT", "image/x-bil;bits=32")
             .parameter("STYLES", "")
             .build();
+
+        // Compute voxel size in map coordinates
+        try {
+            // Beware, we use world center voxel to get voxel size.
+            // This would cause stitching problems if we ever want to merge two adjacent generated worlds.
+            ReferencedEnvelope voxelEnvelope = envelopeProvider.computeForCRS(crs, WorldBBox3d.ORIGIN);
+            voxelSizeX = voxelEnvelope.getWidth();
+            voxelSizeY = voxelEnvelope.getHeight();
+
+        } catch (FactoryException | TransformException e) {
+            throw new IllegalStateException("Unable to compute voxel size in map coordinate", e);
+        }
     }
 
     @Override
@@ -72,27 +89,21 @@ public Result<FloatGeographicDataMatrix2d> provide(WorldBBox3d bbox) throws Gene
             throw new GenerationFailedException(e);
         }
 
-        // Pixel size in map units
-        // TODO: Should be computed from capabilities and voxel size in realworld
-        // (we don't need information more accurate than voxel size neither information more
-        // accurate than capabilities)
-        double pixelSize = 1;
-
         // This is the WMS bbox expressed in map coordinates.
         // It is used below to deduce matrix offset and cell size.
         // WMS matrix is aligned in the same way in all tiles (use of floor/ceil).
         // This prevents glitches between tiles.
 
-        // We need margin for interpolation (-1/+1 expressed in pixelSize)
+        // We need margin for interpolation (-1/+1 expressed in voxel size)
         // TODO: Margin size should come from processor (may be with PR#123?)
-        double minX = Rounding.floor(envelope.getMinX(), pixelSize, -1);
-        double minY = Rounding.floor(envelope.getMinY(), pixelSize, -1);
-        double maxX = Rounding.ceil(envelope.getMaxX(), pixelSize, 1);
-        double maxY = Rounding.ceil(envelope.getMaxY(), pixelSize, 1);
+        double minX = Rounding.floor(envelope.getMinX(), voxelSizeX, -1);
+        double minY = Rounding.floor(envelope.getMinY(), voxelSizeX, -1);
+        double maxX = Rounding.ceil(envelope.getMaxX(), voxelSizeY, 1);
+        double maxY = Rounding.ceil(envelope.getMaxY(), voxelSizeY, 1);
 
         // Formulas give integer numbers, we round them to avoid surprises with floating points
-        int width  = (int) Math.round((maxX - minX) / pixelSize);
-        int height = (int) Math.round((maxY - minY) / pixelSize);
+        int width  = (int) Math.round((maxX - minX) / voxelSizeX);
+        int height = (int) Math.round((maxY - minY) / voxelSizeY);
 
         // Perform WMS query
         ParameterizedURL url = baseURL.builder()
@@ -130,7 +141,14 @@ public Result<FloatGeographicDataMatrix2d> provide(WorldBBox3d bbox) throws Gene
         if (total != size)
             throw new RetryableException("Incomplete data read from stream");
 
-        FloatArrayGeographicDataMatrix2d result = new FloatArrayGeographicDataMatrix2d(width, height, minX, minY, pixelSize, pixelSize);
+        FloatArrayGeographicDataMatrix2d result = new FloatArrayGeographicDataMatrix2d(
+            width,
+            height,
+            minX,
+            minY,
+            voxelSizeX,
+            voxelSizeY
+        );
 
         // Decode binary data into float matrix
         ByteBuffer.wrap(data).order(ByteOrder.LITTLE_ENDIAN).asFloatBuffer().get(result.data());

src/main/java/fr/ign/voxatile/core/models/FloatMatrixModel.java

@@ -35,9 +35,9 @@ public FloatMatrixModel(FloatGeographicDataMatrix2d data, MapToWorldConverter co
 
         this.data = data;
 
-        double maxX = data.offsetX() + data.sizeX() * data.cellSizeX() - 1.0;
-        double maxY = data.offsetY() + data.sizeY() * data.cellSizeY() - 1.0;
-
+        // Size - 1 is the maximum coordinate value.
+        double maxX = data.offsetX() + (data.sizeX() - 1) * data.cellSizeX();
+        double maxY = data.offsetY() + (data.sizeY() - 1) * data.cellSizeY();
         this.bbox = new WorldBBox2d(
             mapToWorld.convert(new MapCoordinates(data.offsetX(), data.offsetY())),
             mapToWorld.convert(new MapCoordinates(data.offsetX(), maxY)),
@@ -63,8 +63,8 @@ public Float get(WorldCoords2d coords) {
         }
 
         // We interpolate between cell centers (not cell upper left corner)
-        float x = (float) ((coordinates.x() + 0.5 - data.offsetX()) / data.cellSizeX() - 0.5);
-        float y = (float) ((coordinates.y() + 0.5 - data.offsetY()) / data.cellSizeY() - 0.5);
+        double x = (coordinates.x() - data.offsetX()) / data.cellSizeX();
+        double y = (coordinates.y() - data.offsetY()) / data.cellSizeY();
 
         // Using separate ceil & floor allows a good management of integer coordinates
         int xf = (int) Math.floor(x);
@@ -75,8 +75,8 @@ public Float get(WorldCoords2d coords) {
             return null;
 
         // Basic bilinear interpolation
-        float fx = x - xf;
-        float fy = y - yf;
+        float fx = (float) (x - xf);
+        float fy = (float) (y - yf);
 
         return (1 - fy) * ((1 - fx) * data.getFloat(xf, yf) + fx * data.getFloat(xc, yf))
             + fy * ((1 - fx) * data.getFloat(xf, yc) + fx * data.getFloat(xc, yc));

src/test/java/fr/ign/voxatile/core/models/FloatMatrixModelTest.java

@@ -66,21 +66,21 @@ public void testGetInterpolation() throws TransformException {
         // Beware, interpolation is between cells centers, not cells upper left corner, at voxel center, not voxel upper left corner
         // So we have an offset of -5/-5 for a cell size of 10
         // and 0.5/0.5 for voxel size in order to have valid interpolable values between [0-10],[0-10].
-        FloatGeographicDataMatrix2d data = new FloatArrayGeographicDataMatrix2d(values, 2, 2, -4.5, -4.5, 10.0, 10.0);
+        FloatGeographicDataMatrix2d data = new FloatArrayGeographicDataMatrix2d(values, 2, 2, -5.0, -5.0, 10.0, 10.0);
 
         FloatMatrixModel model = new FloatMatrixModel(data, converter);
         // Borders
-        assertNull(model.get(new WorldCoords2d(0, 11)));
-        assertNull(model.get(new WorldCoords2d(11, 0)));
-        assertEquals(-3.0f, model.get(new WorldCoords2d(0, 0)));
-        assertEquals(-1.0f, model.get(new WorldCoords2d(0, 10)));
-        assertEquals(5.0f, model.get(new WorldCoords2d(10, 0)));
-        assertEquals(1.0f, model.get(new WorldCoords2d(10, 10)));
+        assertNull(model.get(new WorldCoords2d(0, 6)));
+        assertNull(model.get(new WorldCoords2d(6, 0)));
+        assertEquals(-3.0f, model.get(new WorldCoords2d(-5, -5)));
+        assertEquals(-1.0f, model.get(new WorldCoords2d(-5, 5)));
+        assertEquals(5.0f, model.get(new WorldCoords2d(5, -5)));
+        assertEquals(1.0f, model.get(new WorldCoords2d(5, 5)));
         // Centers
-        assertEquals(1.0f, model.get(new WorldCoords2d(5, 0)));
-        assertEquals(0.0f, model.get(new WorldCoords2d(5, 10)));
-        assertEquals(-2.0f, model.get(new WorldCoords2d(0, 5)));
-        assertEquals(3.0f, model.get(new WorldCoords2d(10, 5)));
-        assertEquals(0.5f, model.get(new WorldCoords2d(5, 5)));
+        assertEquals(1.0f, model.get(new WorldCoords2d(0, -5)));
+        assertEquals(0.0f, model.get(new WorldCoords2d(0, 5)));
+        assertEquals(-2.0f, model.get(new WorldCoords2d(-5, 0)));
+        assertEquals(3.0f, model.get(new WorldCoords2d(5, 0)));
+        assertEquals(0.5f, model.get(new WorldCoords2d(0, 0)));
     }
 }