diff --git a/src/main/java/fr/ign/voxatile/core/inputs/FloatArrayGeographicDataMatrix2d.java b/src/main/java/fr/ign/voxatile/core/inputs/FloatArrayGeographicDataMatrix2d.java index 47fd51bb5..31e621420 100644 --- a/src/main/java/fr/ign/voxatile/core/inputs/FloatArrayGeographicDataMatrix2d.java +++ b/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]; } } diff --git a/src/main/java/fr/ign/voxatile/core/inputs/GeoTiffDataProvider.java b/src/main/java/fr/ign/voxatile/core/inputs/GeoTiffDataProvider.java index b8ba6d74c..ab8206d50 100644 --- a/src/main/java/fr/ign/voxatile/core/inputs/GeoTiffDataProvider.java +++ b/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 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)); diff --git a/src/main/java/fr/ign/voxatile/core/inputs/WMSFloatBilDataProvider.java b/src/main/java/fr/ign/voxatile/core/inputs/WMSFloatBilDataProvider.java index e20135f1d..58ca4c3fe 100644 --- a/src/main/java/fr/ign/voxatile/core/inputs/WMSFloatBilDataProvider.java +++ b/src/main/java/fr/ign/voxatile/core/inputs/WMSFloatBilDataProvider.java @@ -31,6 +31,10 @@ public class WMSFloatBilDataProvider implements Provider 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 +140,14 @@ public Result 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()); diff --git a/src/main/java/fr/ign/voxatile/core/models/FloatMatrixModel.java b/src/main/java/fr/ign/voxatile/core/models/FloatMatrixModel.java index d683cf8cf..1c839a709 100644 --- a/src/main/java/fr/ign/voxatile/core/models/FloatMatrixModel.java +++ b/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)); diff --git a/src/test/java/fr/ign/voxatile/core/models/FloatMatrixModelTest.java b/src/test/java/fr/ign/voxatile/core/models/FloatMatrixModelTest.java index 948ed9864..db78d6dd0 100644 --- a/src/test/java/fr/ign/voxatile/core/models/FloatMatrixModelTest.java +++ b/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))); } }