You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sis.apache.org by de...@apache.org on 2020/12/11 17:37:30 UTC
[sis] 03/03: First version of a contouring algorithm,
adapted from Johann Sorel work based on "marching squares"
algorithm. This implementation builds isolines for all bands and all levels
in a single iteration over image pixels. Java2D shapes are constructed
on-the-fly during iteration.
This is an automated email from the ASF dual-hosted git repository.
desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git
commit 4f223fe9523e4363fbb97d832fafc5a4df74d069
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Mon Dec 7 19:25:11 2020 +0100
First version of a contouring algorithm, adapted from Johann Sorel work based on "marching squares" algorithm.
This implementation builds isolines for all bands and all levels in a single iteration over image pixels.
Java2D shapes are constructed on-the-fly during iteration.
https://issues.apache.org/jira/browse/SIS-506
---
.../internal/processing/image/IsolineTracer.java | 697 +++++++++++++++++++++
.../sis/internal/processing/image/Isolines.java | 259 ++++++++
.../internal/processing/image/package-info.java | 32 +
.../internal/processing/image/IsolinesTest.java | 308 +++++++++
.../apache/sis/test/suite/FeatureTestSuite.java | 3 +-
5 files changed, 1298 insertions(+), 1 deletion(-)
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/IsolineTracer.java b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/IsolineTracer.java
new file mode 100644
index 0000000..056fb73
--- /dev/null
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/IsolineTracer.java
@@ -0,0 +1,697 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.internal.processing.image;
+
+import java.util.Arrays;
+import java.nio.DoubleBuffer;
+import java.awt.geom.Path2D;
+import java.awt.geom.PathIterator;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.util.ArraysExt;
+
+
+/**
+ * Iterator over contouring grid cells together with an interpolator and an assembler of polyline segments.
+ * A single instance of this class is created by {@code Isolines.generate(…)} for all bands to process in a
+ * given image. {@code IsolineTracer} is used for doing a single iteration over all image pixels.
+ *
+ * @author Johann Sorel (Geomatys)
+ * @author Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ *
+ * @see <a href="https://en.wikipedia.org/wiki/Marching_squares">Marching squares on Wikipedia</a>
+ *
+ * @since 1.1
+ * @module
+ */
+final class IsolineTracer {
+ /**
+ * Mask to apply on {@link Level#isDataAbove} for telling that value in a corner is higher than the level value.
+ * Values are defined in {@code PixelIterator.Window} iteration order: from left to right, then top to bottom.
+ *
+ * <p>Note: there is some hard-coded dependencies to those exact values.
+ * If values are changed, search for example for {@code log2(UPPER_RIGHT)} in comments.</p>
+ */
+ static final int UPPER_LEFT = 1, UPPER_RIGHT = 2, LOWER_LEFT = 4, LOWER_RIGHT = 8;
+
+ /**
+ * The 2×2 window containing pixel values in the 4 corners of current contouring grid cell.
+ * Values are always stored with band index varying fastest, then column index, then row index.
+ * Capacity and limit of data buffer is <var>(number of bands)</var> × 2 (width) × 2 (height).
+ */
+ private final DoubleBuffer window;
+
+ /**
+ * Increment to the position for reading next sample value.
+ * It corresponds to the number of bands in {@link #window}.
+ */
+ private final int pixelStride;
+
+ /**
+ * Pixel coordinate on the left side of the cell where to interpolate.
+ */
+ int x;
+
+ /**
+ * Pixel coordinate on the top side of the cell where to interpolate.
+ */
+ int y;
+
+ /**
+ * Threshold for considering two coordinates as equal.
+ * Shall be a value between 0 and 0.5.
+ */
+ private final double tolerance;
+
+ /**
+ * Final transform to apply on coordinates.
+ */
+ private final MathTransform gridToCRS;
+
+ /**
+ * Creates a new position for the given data window.
+ */
+ IsolineTracer(final DoubleBuffer window, final int pixelStride, double tolerance, final MathTransform gridToCRS) {
+ this.window = window;
+ this.pixelStride = pixelStride;
+ this.tolerance = (tolerance = Math.min(Math.abs(tolerance), 0.5)) >= 0 ? tolerance : 0;
+ this.gridToCRS = gridToCRS;
+ }
+
+ /**
+ * Builder of polylines for a single level. The segments to create are determined by a set
+ * of {@linkplain #isDataAbove four flags} (one for each corner) encoded in an integer.
+ * The meaning of those flags is described in Wikipedia "Marching squares" article,
+ * except that this implementation uses different values.
+ */
+ final class Level {
+ /**
+ * The level value. This is a copy of {@link Isolines#levelValues} at the index of this level.
+ *
+ * @see #interpolate(int, int)
+ */
+ final double value;
+
+ /**
+ * Bitset telling which corners have a value greater than this isoline level {@linkplain #value}.
+ * Each corner is associated to one of the bits illustrated below, where bit (0) is the less significant.
+ * Note that this bit order is different than the order used in Wikipedia "Marching squares" article.
+ * The order used in this class allows more direct bitwise operations as described in next section.
+ *
+ * {@preformat text
+ * (0)╌╌╌(1)
+ * ╎ ╎
+ * (2)╌╌╌(3)
+ * }
+ *
+ * Bits are set to 1 where the data value is above the isoline {@linkplain #value}, and 0 where the data
+ * value is equal or below the isoline value. Data values exactly equal to the isoline value are handled
+ * as if they were greater. It does not matter for interpolations: we could flip this convention randomly,
+ * the interpolated points would still the same. It could change the way line segments are assembled in a
+ * single {@link Polyline}, but the algorithm stay consistent if we always apply the same rule for all points.
+ *
+ * <h4>Reusing bits from previous iteration</h4>
+ * We will iterate on pixels from left to right, then from top to bottom. With that iteration order,
+ * bits 0 and 2 can be obtained from the bit pattern of previous iteration with a simple bit shift.
+ *
+ * @see #UPPER_LEFT
+ * @see #UPPER_RIGHT
+ * @see #LOWER_LEFT
+ * @see #LOWER_RIGHT
+ */
+ int isDataAbove;
+
+ /**
+ * The polyline to be continued on the next column. This is a single instance because iteration happens
+ * from left to right before top to bottom. This instance is non-empty if the cell in previous iteration
+ * was like below (all those examples have a line crossing the right border):
+ *
+ * {@preformat text
+ * ●╌╌╌╌╌╌● ○╌╱╌╌╌╌●╱ ○╌╌╌╌╲╌●
+ * ╎ ╎ ╎╱ ╱ ╎ ╲╎
+ * ─┼──────┼─ ╱ ╱╎ ╎ ╲
+ * ○╌╌╌╌╌╌○ ╱●╌╌╌╌╱╌○ ○╌╌╌╌╌╌○╲
+ * }
+ *
+ * This field {@linkplain Polyline#isEmpty() is empty} if the cell in previous iteration was like below
+ * (no line cross the right border):
+ *
+ * {@preformat text
+ * ○╌╲╌╌╌╌● ○╌╌╌┼╌╌●
+ * ╎ ╲ ╎ ╎ │ ╎
+ * ╎ ╲ ╎ ╎ │ ╎
+ * ○╌╌╌╌╲╌● ○╌╌╌┼╌╌●
+ * }
+ */
+ private final Polyline polylineOnLeft;
+
+ /**
+ * The polylines in each column which need to be continued on the next row.
+ * This array contains empty instances in columns where there is no polyline to continue on next row.
+ * For non-empty element at index <var>x</var>, values on the left border are given by pixels at coordinate
+ * {@code x} and values on the right border are given by pixels at coordinate {@code x+1}. Example:
+ *
+ * {@preformat text
+ * ○╌╌╌╌╌╌●╱
+ * ╎ Top ╱
+ * ╎ [x] ╱╎
+ * ●╌╌╌╌╌╌●╌╌╌╌╱╌○
+ * ╎ Left ╎██████╎ ← Cell where to create a segment
+ * ─┼──────┼██████╎
+ * ○╌╌╌╌╌╌○╌╌╌╌╌╌○
+ * ↑
+ * x coordinate of first pixel (upper-left corner)
+ * }
+ */
+ private final Polyline[] polylinesOnTop;
+
+ /**
+ * The isolines as a Java2D shape, created when first needed. The {@link Polyline} coordinates are copied in
+ * this path when a geometry is closed or when iteration finished on a row and the polyline is not reused by
+ * next row. This is the shape to be returned to user for this level after we finished to process all cells.
+ */
+ Path2D path;
+
+ /**
+ * Creates new isoline levels for the given value.
+ *
+ * @param value the isoline level value.
+ * @param width the contouring grid cell width (one cell smaller than image width).
+ */
+ Level(final double value, final int width) {
+ this.value = value;
+ polylineOnLeft = new Polyline();
+ polylinesOnTop = new Polyline[width];
+ for (int i=0; i<width; i++) {
+ polylinesOnTop[i] = new Polyline();
+ }
+ }
+
+ /**
+ * Initializes the {@link #isDataAbove} value with values for the column on the right side.
+ * After this method call, the {@link #UPPER_RIGHT} and {@link #LOWER_RIGHT} bits still need to be set.
+ */
+ final void nextColumn() {
+ /*
+ * Move bits on the right side to the left side.
+ * The 1 operand in >>> is the hard-coded value
+ * of log2(UPPER_RIGHT) - log2(UPPER_LEFT)
+ * and log2(LOWER_RIGHT) - log2(LOWER_LEFT).
+ */
+ isDataAbove = (isDataAbove & (UPPER_RIGHT | LOWER_RIGHT)) >>> 1;
+ }
+
+ /**
+ * Adds segments computed for values in a single pixel. Interpolations are determined by the 4 lowest bits
+ * of {@link #isDataAbove}. The {@link #polylineOnLeft} and {@code polylinesOnTop[x]} elements are updated
+ * by this method.
+ *
+ * <h4>How NaN values are handled</h4>
+ * This algorithm does not need special attention for {@link Double#NaN} values. Interpolations will produce
+ * {@code NaN} values and append them to the correct polyline (which does not depend on interpolation result)
+ * like real values. Those NaN values will be filtered later in another method, when copying coordinates in
+ * {@link Path2D} objects.
+ */
+ final void interpolate() throws TransformException {
+ switch (isDataAbove) {
+ default: {
+ throw new AssertionError(isDataAbove); // Should never happen.
+ }
+ /* ○╌╌╌╌╌╌○ ●╌╌╌╌╌╌●
+ * ╎ ╎ ╎ ╎
+ * ╎ ╎ ╎ ╎
+ * ○╌╌╌╌╌╌○ ●╌╌╌╌╌╌●
+ */
+ case 0:
+ case UPPER_LEFT | UPPER_RIGHT | LOWER_LEFT | LOWER_RIGHT: {
+ assert polylinesOnTop[x].isEmpty();
+ assert polylineOnLeft .isEmpty();
+ break;
+ }
+ /* ○╌╌╌╌╌╌○ ●╌╌╌╌╌╌●
+ * ─┼──────┼─ ─┼──────┼─
+ * ╎ ╎ ╎ ╎
+ * ●╌╌╌╌╌╌● ○╌╌╌╌╌╌○
+ */
+ case LOWER_LEFT | LOWER_RIGHT:
+ case UPPER_LEFT | UPPER_RIGHT: {
+ assert polylinesOnTop[x].isEmpty();
+ if (polylineOnLeft.isEmpty()) {
+ interpolateOnLeftSide();
+ }
+ interpolateOnRightSide(polylineOnLeft); // Will be the left side of next column.
+ break;
+ }
+ /* ○╌╌╌┼╌╌● ●╌╌╌┼╌╌○
+ * ╎ │ ╎ ╎ │ ╎
+ * ╎ │ ╎ ╎ │ ╎
+ * ○╌╌╌┼╌╌● ●╌╌╌┼╌╌○
+ */
+ case UPPER_RIGHT | LOWER_RIGHT:
+ case UPPER_LEFT | LOWER_LEFT: {
+ assert polylineOnLeft.isEmpty();
+ final Polyline polylineOnTop = polylinesOnTop[x];
+ if (polylineOnTop.isEmpty()) {
+ interpolateOnTopSide(polylineOnTop);
+ }
+ interpolateOnBottomSide(polylineOnTop); // Will be top side of next row.
+ break;
+ }
+ /* ╲○╌╌╌╌╌╌○ ╲●╌╌╌╌╌╌●
+ * ╲ ╎ ╲ ╎
+ * ╎╲ ╎ ╎╲ ╎
+ * ●╌╲╌╌╌╌○ ○╌╲╌╌╌╌●
+ */
+ case LOWER_LEFT:
+ case UPPER_LEFT | UPPER_RIGHT | LOWER_RIGHT: {
+ assert polylinesOnTop[x].isEmpty();
+ if (polylineOnLeft.isEmpty()) {
+ interpolateOnLeftSide();
+ }
+ interpolateOnBottomSide(polylinesOnTop[x].transferFrom(polylineOnLeft));
+ break;
+ }
+ /* ○╌╌╌╌╲╌● ●╌╌╌╌╲╌○
+ * ╎ ╲╎ ╎ ╲╎
+ * ╎ ╲ ╎ ╲
+ * ○╌╌╌╌╌╌○╲ ●╌╌╌╌╌╌●╲
+ */
+ case UPPER_RIGHT:
+ case UPPER_LEFT | LOWER_LEFT | LOWER_RIGHT: {
+ assert polylineOnLeft.isEmpty();
+ if (polylineOnLeft.transferFrom(polylinesOnTop[x]).isEmpty()) {
+ interpolateOnTopSide(polylineOnLeft);
+ }
+ interpolateOnRightSide(polylineOnLeft);
+ break;
+ }
+ /* ○╌╌╌╌╌╌○╱ ●╌╌╌╌╌╌●╱
+ * ╎ ╱ ╎ ╱
+ * ╎ ╱╎ ╎ ╱╎
+ * ○╌╌╌╌╱╌● ●╌╌╌╌╱╌○
+ */
+ case LOWER_RIGHT:
+ case UPPER_LEFT | UPPER_RIGHT | LOWER_LEFT: {
+ assert polylinesOnTop[x].isEmpty();
+ assert polylineOnLeft .isEmpty();
+ interpolateOnRightSide (polylineOnLeft);
+ interpolateOnBottomSide(polylinesOnTop[x].attach(polylineOnLeft));
+ // Bottom of this cell will be top of next row.
+ break;
+ }
+ /* ●╌╱╌╌╌╌○ ○╌╱╌╌╌╌●
+ * ╎╱ ╎ ╎╱ ╎
+ * ╱ ╎ ╱ ╎
+ * ╱○╌╌╌╌╌╌○ ╱●╌╌╌╌╌╌●
+ */
+ case UPPER_LEFT:
+ case UPPER_RIGHT | LOWER_LEFT | LOWER_RIGHT: {
+ if (polylineOnLeft.isEmpty()) {
+ interpolateOnLeftSide();
+ }
+ interpolateOnTopSide(polylineOnLeft);
+ path = close(polylineOnLeft, polylinesOnTop[x], path);
+ break;
+ }
+ /* ○╌╱╌╌╌╌●╱ ╲●╌╌╌╌╲╌○
+ * ╎╱ ╱ ╲ ╲╎
+ * ╱ ╱╎ ╎╲ ╲
+ * ╱●╌╌╌╌╱╌○ ○╌╲╌╌╌╌●╲
+ *
+ * Disambiguation of saddle points: use the average data value for the center of the cell.
+ * If the estimated center value is greater than the isoline value, the above drawings are
+ * okay and we do not need to change `isDataAbove`. This is the left side illustrated below.
+ * But if the center value is below isoline value, then we need to flip `isDataAbove` bits
+ * (conceptually; not really because we need to keep `isDataAbove` value for next iteration).
+ * This is the right side illustrated below.
+ *
+ * ○╱╌╌●╱ ╲●╌╌╲○ ╲○╌╌╲● ●╱╌╌○╱
+ * ╱ ● ╱ ╲ ● ╲ ╲ ○ ╲ ╱ ○ ╱
+ * ╱●╌╌╱○ ○╲╌╌●╲ ●╲╌╌○╲ ╱○╌╌╱●
+ */
+ case UPPER_RIGHT | LOWER_LEFT:
+ case UPPER_LEFT | LOWER_RIGHT: {
+ double average = 0;
+ { // Compute sum of 4 corners.
+ final DoubleBuffer data = window;
+ final int limit = data.limit();
+ int p = data.position();
+ do average += data.get(p);
+ while ((p += pixelStride) < limit);
+ average /= 4;
+ }
+ boolean LLtoUR = isDataAbove == (LOWER_LEFT | UPPER_RIGHT);
+ LLtoUR ^= (average <= value);
+ if (polylineOnLeft.isEmpty()) {
+ interpolateOnLeftSide();
+ }
+ if (LLtoUR) {
+ interpolateOnTopSide(polylineOnLeft);
+ path = close(polylineOnLeft, polylinesOnTop[x], path);
+ interpolateOnRightSide (polylineOnLeft);
+ interpolateOnBottomSide(polylinesOnTop[x].attach(polylineOnLeft));
+ } else {
+ final Polyline swap = new Polyline();
+ final Polyline polylineOnTop = polylinesOnTop[x];
+ if (swap.transferFrom(polylineOnTop).isEmpty()) {
+ interpolateOnTopSide(swap);
+ }
+ interpolateOnRightSide(swap);
+ interpolateOnBottomSide(polylineOnTop.transferFrom(polylineOnLeft));
+ polylineOnLeft.transferFrom(swap);
+ }
+ break;
+ }
+ }
+ }
+
+ /** Appends to {@link #polylineOnLeft} a point interpolated on the left side. */
+ private void interpolateOnLeftSide() {
+ polylineOnLeft.append(x, y + interpolate(0, 2 * pixelStride));
+ }
+
+ /** Appends to the given polyline a point interpolated on the right side. */
+ private void interpolateOnRightSide(final Polyline appendTo) {
+ appendTo.append(x + 1, y + interpolate(pixelStride, 3 * pixelStride));
+ }
+
+ /** Appends to the given polyline a point interpolated on the top side. */
+ private void interpolateOnTopSide(final Polyline appendTo) {
+ appendTo.append(x + interpolate(0, pixelStride), y);
+ }
+
+ /** Appends to the given polyline a point interpolated on the bottom side. */
+ private void interpolateOnBottomSide(final Polyline appendTo) {
+ appendTo.append(x + interpolate(2 * pixelStride, 3 * pixelStride), y + 1);
+ }
+
+ /**
+ * Interpolates the position where the isoline passes between two values.
+ * The {@link #window} buffer position shall be the first sample value
+ * for the band to process.
+ *
+ * @param i1 index of first value in the buffer, ignoring band offset.
+ * @param i2 index of second value in the buffer, ignoring band offset.
+ * @return a value interpolated between the values at the two given indices.
+ */
+ private double interpolate(final int i1, final int i2) {
+ final DoubleBuffer data = window;
+ final int p = data.position();
+ final double v1 = data.get(p + i1);
+ final double v2 = data.get(p + i2);
+ return (value - v1) / (v2 - v1);
+ }
+
+ /**
+ * Invoked after iteration on a single row has been completed. If there is a polyline
+ * finishing on the right image border, that polyline needs to be written now because
+ * it will not be continued by cells on next rows.
+ */
+ final void finishedRow() throws TransformException {
+ if (!polylineOnLeft.transferToOpposite()) {
+ path = writeTo(polylineOnLeft, path);
+ }
+ isDataAbove = 0;
+ }
+
+ /**
+ * Invoked after the iteration has been completed on the full image.
+ * This method flushes all reminding polylines to the {@link #path}.
+ * It assumes that {@link #finishedRow()} has already been invoked.
+ */
+ final void finish() throws TransformException {
+ for (int i=0; i < polylinesOnTop.length; i++) {
+ path = writeTo(polylinesOnTop[i], path);
+ polylinesOnTop[i] = null;
+ }
+ }
+ }
+
+ /**
+ * Coordinates of a polyline under construction. Coordinates can be appended in only one direction.
+ * If the polyline may growth on both directions (which happens if the polyline crosses the bottom
+ * side and the right side of a cell), then the two directions are handled by two distinct instances
+ * connected by their {@link #opposite} field.
+ *
+ * <p>When a polyline has been completed, its content is copied to {@link Level#path}
+ * and the {@code Polyline} object is recycled for a new polyline.</p>
+ */
+ private static final class Polyline {
+ /**
+ * Number of coordinates in a tuple.
+ */
+ private static final int DIMENSION = 2;
+
+ /**
+ * Coordinates as (x,y) tuples. This array is expanded as needed.
+ */
+ private double[] coordinates;
+
+ /**
+ * Number of valid elements in the {@link #coordinates} array.
+ * This is twice the number of points.
+ */
+ private int size;
+
+ /**
+ * If the polyline has points added to its two extremities, the other extremity.
+ * Otherwise {@code null}.
+ */
+ private Polyline opposite;
+
+ /**
+ * Creates an initially empty polyline.
+ */
+ Polyline() {
+ coordinates = ArraysExt.EMPTY_DOUBLE;
+ }
+
+ /**
+ * Discards all coordinates in this polyline.
+ */
+ final void clear() {
+ opposite = null;
+ size = 0;
+ }
+
+ /**
+ * Returns whether this polyline is empty.
+ */
+ final boolean isEmpty() {
+ return size == 0;
+ }
+
+ /**
+ * Declares that the specified polyline will add points in the direction opposite to this polyline.
+ * This happens when the polyline crosses the bottom side and the right side of a cell (assuming an
+ * iteration from left to right and top to bottom).
+ *
+ * @return {@code this} for method calls chaining.
+ */
+ final Polyline attach(final Polyline other) {
+ assert (opposite == null) & (other.opposite == null);
+ other.opposite = this;
+ opposite = other;
+ return this;
+ }
+
+ /**
+ * Transfers all coordinates from given polylines to this polylines, in same order.
+ * This is used when polyline on the left side continues on bottom side,
+ * or conversely when polyline on the top side continues on right side.
+ * This polyline shall be empty before this method is invoked.
+ * The given source will become empty after this method returned.
+ *
+ * @param source the source from which to take data.
+ * @return {@code this} for method calls chaining.
+ */
+ final Polyline transferFrom(final Polyline source) {
+ assert isEmpty();
+ final double[] buffer = coordinates;
+ coordinates = source.coordinates;
+ size = source.size;
+ opposite = source.opposite;
+ if (opposite != null) {
+ opposite.opposite = this;
+ }
+ source.clear();
+ source.coordinates = buffer;
+ return this;
+ }
+
+ /**
+ * Transfers all coordinates from this polyline to the polyline going in opposite direction.
+ * This is used when this polyline reached the right image border, in which case its data
+ * will be lost if we do not copy them somewhere.
+ *
+ * @return {@code true} if coordinates have been transferred,
+ * or {@code false} if there is no opposite direction.
+ */
+ final boolean transferToOpposite() {
+ if (opposite == null) {
+ return false;
+ }
+ final int sum = size + opposite.size;
+ double[] data = opposite.coordinates;
+ if (sum > data.length) {
+ data = new double[sum];
+ }
+ System.arraycopy(opposite.coordinates, 0, data, size, opposite.size);
+ for (int i=0, t=size; (t -= DIMENSION) >= 0;) {
+ data[t ] = coordinates[i++];
+ data[t+1] = coordinates[i++];
+ }
+ opposite.size = sum;
+ opposite.coordinates = data;
+ opposite.opposite = null;
+ clear();
+ return true;
+ }
+
+ /**
+ * Appends given coordinates to this polyline.
+ *
+ * @param x first coordinate of the (x,y) tuple to add.
+ * @param y second coordinate of the (x,y) tuple to add.
+ */
+ final void append(final double x, final double y) {
+ if (size >= coordinates.length) {
+ coordinates = Arrays.copyOf(coordinates, Math.max(size * 2, 32));
+ }
+ coordinates[size++] = x;
+ coordinates[size++] = y;
+ }
+
+ /**
+ * Returns a string representation for debugging purposes.
+ */
+ @Override
+ public String toString() {
+ final StringBuilder b = new StringBuilder(30).append('[');
+ if (size >= DIMENSION) {
+ b.append((float) coordinates[0]).append(", ").append((float) coordinates[1]);
+ final int n = size - DIMENSION;
+ if (n >= DIMENSION) {
+ b.append(", ");
+ if (size >= DIMENSION*3) {
+ b.append(" … (").append(size / DIMENSION).append(" pts) … ");
+ }
+ b.append((float) coordinates[n]).append(", ").append((float) coordinates[n+1]);
+ }
+ }
+ return b.append(']').toString();
+ }
+ }
+
+ /**
+ * Creates a polygon by joining first polyline with the specified polyline, then writes
+ * it to the specified path. The two polylines will become empty after this method call.
+ *
+ * @param join the other polyline to join.
+ * @param path where to write the polygon, or {@code null} if not yet created.
+ * @return the given path, or a newly created path if the argument was null.
+ */
+ final Path2D close(final Polyline polyline, final Polyline join, final Path2D path) throws TransformException {
+ return writeTo(path, polyline.opposite, polyline, join, join.opposite);
+ }
+
+ /**
+ * Writes the content of this polyline to the given path.
+ * This polyline will become empty after this method call.
+ *
+ * @param path where to write the polylines, or {@code null} if not yet created.
+ * @return the given path, or a newly created path if the argument was null.
+ */
+ final Path2D writeTo(final Polyline polyline, final Path2D path) throws TransformException {
+ return writeTo(path, polyline.opposite, polyline);
+ }
+
+ /**
+ * Writes all given polylines to the specified path. Null {@code Polyline} instances are ignored.
+ * {@code Polyline} instances at even index are written with their coordinates in reverse order.
+ *
+ * @param path where to write the polylines, or {@code null} if not yet created.
+ * @param polylines the polylines to write.
+ * @return the given path, or a newly created path if the argument was null.
+ */
+ private Path2D writeTo(Path2D path, final Polyline... polylines) throws TransformException {
+ double xo = Double.NaN;
+ double yo = Double.NaN;
+ double px = Double.NaN;
+ double py = Double.NaN;
+ int state = PathIterator.SEG_MOVETO;
+ for (int pi=0; pi < polylines.length; pi++) {
+ final Polyline p = polylines[pi];
+ if (p == null) {
+ continue;
+ }
+ final boolean reverse = (pi & 1) == 0;
+ final double[] coordinates = p.coordinates;
+ final int size = p.size;
+ gridToCRS.transform(coordinates, 0, coordinates, 0, size / Polyline.DIMENSION);
+ for (int i=0; i<size;) {
+ final double x, y;
+ if (reverse) {
+ y = coordinates[size - ++i];
+ x = coordinates[size - ++i];
+ } else {
+ x = coordinates[i++];
+ y = coordinates[i++];
+ }
+ if (Double.isNaN(x) || Double.isNaN(y) || (Math.abs(x - px) <= tolerance && Math.abs(y - py) <= tolerance)) {
+ continue;
+ }
+ switch (state) {
+ case PathIterator.SEG_MOVETO: {
+ px = xo = x;
+ py = yo = y;
+ state = PathIterator.SEG_LINETO;
+ break;
+ }
+ case PathIterator.SEG_LINETO: {
+ if (path == null) {
+ int s = size - (i - 2*Polyline.DIMENSION);
+ for (int k=pi; ++k < polylines.length;) {
+ final Polyline next = polylines[k];
+ if (next != null) {
+ s = Math.addExact(s, next.size);
+ }
+ }
+ path = new Path2D.Double(Path2D.WIND_NON_ZERO, s / Polyline.DIMENSION);
+ }
+ path.moveTo(xo, yo);
+ path.lineTo(x, y);
+ state = PathIterator.SEG_CLOSE;
+ break;
+ }
+ default: {
+ if (Math.abs(x - xo) <= tolerance && Math.abs(y - yo) <= tolerance) {
+ path.closePath();
+ state = PathIterator.SEG_MOVETO;
+ } else {
+ path.lineTo(x, y);
+ }
+ break;
+ }
+ }
+ }
+ p.clear();
+ }
+ return path;
+ }
+}
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/Isolines.java b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/Isolines.java
new file mode 100644
index 0000000..cb329e3
--- /dev/null
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/Isolines.java
@@ -0,0 +1,259 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.internal.processing.image;
+
+import java.util.Arrays;
+import java.util.TreeMap;
+import java.util.NavigableMap;
+import java.nio.DoubleBuffer;
+import java.awt.Dimension;
+import java.awt.Shape;
+import java.awt.geom.Path2D;
+import java.awt.image.RenderedImage;
+import org.opengis.coverage.grid.SequenceType;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.image.PixelIterator;
+import org.apache.sis.image.TransferType;
+import org.apache.sis.util.ArgumentChecks;
+import org.apache.sis.util.ArraysExt;
+
+import static org.apache.sis.internal.processing.image.IsolineTracer.UPPER_LEFT;
+import static org.apache.sis.internal.processing.image.IsolineTracer.UPPER_RIGHT;
+import static org.apache.sis.internal.processing.image.IsolineTracer.LOWER_RIGHT;
+
+
+/**
+ * Creates isolines at specified levels from grid values provided in a {@link RenderedImage}.
+ * Isolines are created by calls to the {@link #generate(RenderedImage, double[][], double,
+ * MathTransform)} static method.
+ *
+ * @author Johann Sorel (Geomatys)
+ * @author Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ *
+ * @see <a href="https://en.wikipedia.org/wiki/Marching_squares">Marching squares on Wikipedia</a>
+ *
+ * @since 1.1
+ * @module
+ */
+public final class Isolines {
+ /**
+ * Values for which to compute isolines, sorted in ascending order and without NaN values.
+ * This array is used for {@linkplain Arrays#binarySearch(double[], double) binary searches}.
+ * Each value is associated to an instance in the {@link #levels} array.
+ *
+ * @see IsolineTracer.Level#value
+ */
+ private final double[] levelValues;
+
+ /**
+ * Isoline data for each level.
+ * The length of this array is equal to the {@link #levelValues} array length.
+ */
+ private final IsolineTracer.Level[] levels;
+
+ /**
+ * Creates an initially empty set of isolines for the given levels.
+ * The given array should be a clone of user-provided array because
+ * this constructor may modify it in-place.
+ */
+ private Isolines(final IsolineTracer tracer, final double[] values, final int width) {
+ Arrays.sort(values);
+ int n = values.length;
+ while (n > 0 && Double.isNaN(values[n-1])) n--;
+ for (int i=n; --i>0;) {
+ if (values[i] == values[i-1]) {
+ // Remove duplicated elements. May replace -0 by +0.
+ System.arraycopy(values, i, values, i-1, n-- - i);
+ }
+ }
+ levelValues = ArraysExt.resize(values, n);
+ levels = new IsolineTracer.Level[n];
+ for (int i=0; i<n; i++) {
+ levels[i] = tracer.new Level(levelValues[i], width);
+ }
+ }
+
+ /**
+ * Sets the specified bit on {@link IsolineTracer.Level#isDataAbove} for all levels lower than buffer value.
+ * The {@code window} buffer shall be positioned on the value to read, consistently with {@code bit} value.
+ * After this method call, the buffer position will be incremented to the next band if there is more bands
+ * to read, or to the next pixel otherwise.
+ *
+ * <h4>How strict equalities are handled</h4>
+ * Sample values exactly equal to the isoline value are handled as if they were greater. It does not matter
+ * for interpolations: we could flip this convention randomly, the interpolated points would still the same.
+ * However it could change the way line segments are assembled in a single polyline, but the algorithm stay
+ * consistent if we always apply the same rule for all points.
+ *
+ * <h4>How NaN values are handled</h4>
+ * {@link Double#NaN} sample values are considered higher than any level values. The algorithm does not need
+ * special attention for those values; bit patterns will be computed in a consistent way, and interpolations
+ * will produce NaN values and append them to polylines like real values. Those NaN values will be filtered
+ * out in the final stage, when copying coordinates in {@link Path2D} objects.
+ *
+ * @param data a 2×2 view on pixel values in the image.
+ * @param bit {@value IsolineTracer#UPPER_LEFT}, {@value IsolineTracer#UPPER_RIGHT},
+ * {@value IsolineTracer#LOWER_LEFT} or {@value IsolineTracer#LOWER_RIGHT}.
+ */
+ private void setMaskBit(final DoubleBuffer data, final int bit) {
+ int i = Arrays.binarySearch(levelValues, data.get());
+ if (i < 0) i = ~i; // Really tild, not minus.
+ while (--i >= 0) { // Value is higher than all levels below `i`.
+ levels[i].isDataAbove |= bit;
+ }
+ }
+
+ /**
+ * Generates isolines at the specified levels computed from data provided by the given image.
+ * Isolines will be created for every bands in the given image.
+ *
+ * @param data image providing source values.
+ * @param levels values for which to compute isolines. An array should be provided for each band.
+ * If there is more bands than {@code levels.length}, the last array is reused for
+ * all remaining bands.
+ * @param tolerance threshold for considering two coordinates as equal.
+ * @param gridToCRS transform from pixel coordinates to geometry coordinates, or {@code null} if none.
+ * @return the isolines for each band in the given image.
+ * @throws TransformException if an interpolated point can not be transformed using the given transform.
+ */
+ public static Isolines[] generate(final RenderedImage data, final double[][] levels, final double tolerance,
+ MathTransform gridToCRS) throws TransformException
+ {
+ ArgumentChecks.ensureNonNull("data", data);
+ ArgumentChecks.ensureNonNull("levels", levels);
+ { // For keeping variable locale.
+ final MathTransform gridToImage = MathTransforms.translation(data.getMinX(), data.getMinY());
+ if (gridToCRS == null) {
+ gridToCRS = gridToImage;
+ } else {
+ ArgumentChecks.ensureDimensionsMatch("gridToCRS", 2, 2, gridToCRS);
+ gridToCRS = MathTransforms.concatenate(gridToImage, gridToCRS);
+ }
+ }
+ final PixelIterator iterator = new PixelIterator.Builder().setIteratorOrder(SequenceType.LINEAR)
+ .setWindowSize(new Dimension(2,2)).create(data);
+ /*
+ * A window of size 2×2 pixels over pixel values.
+ */
+ final PixelIterator.Window<DoubleBuffer> window = iterator.createWindow(TransferType.DOUBLE);
+ final DoubleBuffer buffer = window.values;
+ final int numBands = iterator.getNumBands();
+ final IsolineTracer tracer = new IsolineTracer(buffer, numBands, tolerance, gridToCRS);
+ /*
+ * Prepare the set of isolines for each band in the image.
+ * The number of cells on the horizontal axis is one less
+ * than the image width.
+ */
+ final int width = data.getWidth() - 1;
+ final Isolines[] isolines = new Isolines[numBands];
+ { // For keeping variable locale.
+ double[] levelValues = ArraysExt.EMPTY_DOUBLE;
+ for (int b=0; b<numBands; b++) {
+ if (b < levels.length) {
+ levelValues = levels[b];
+ ArgumentChecks.ensureNonNullElement("levels", b, levelValues);
+ levelValues = levelValues.clone();
+ }
+ isolines[b] = new Isolines(tracer, levelValues, width);
+ }
+ }
+ /*
+ * Compute isolines for all bands. Iteration over bands must be the innermost loop because
+ * data layout in buffer is band index varying fastest, then column index, then row index.
+ */
+ final int lastPixel = numBands * 3;
+abort: while (iterator.next()) {
+ /*
+ * First pixel of a new row.
+ */
+ window.update(); // Also reset buffer position to zero.
+ for (int flag = UPPER_LEFT; flag <= LOWER_RIGHT; flag <<= 1) {
+ for (int b=0; b<numBands; b++) { // Must be the inner loop (see above comment).
+ isolines[b].setMaskBit(buffer, flag);
+ }
+ }
+ for (int b=0; b<numBands; b++) {
+ buffer.position(b);
+ for (final IsolineTracer.Level level : isolines[b].levels) {
+ level.interpolate();
+ }
+ }
+ /*
+ * All pixels on a row after the first column. We can reuse the bitmask of previous
+ * iteration with a simple bit shift operation.
+ */
+ for (tracer.x = 1; tracer.x < width; tracer.x++) {
+ if (!iterator.next()) break abort;
+ window.update();
+ for (int b=0; b<numBands; b++) {
+ final Isolines iso = isolines[b];
+ for (final IsolineTracer.Level level : iso.levels) {
+ level.nextColumn();
+ }
+ // TODO! remove cast on JDK9.
+ iso.setMaskBit((DoubleBuffer) buffer.position(numBands + b), UPPER_RIGHT);
+ iso.setMaskBit((DoubleBuffer) buffer.position(lastPixel + b), LOWER_RIGHT);
+ buffer.position(b);
+ for (final IsolineTracer.Level level : iso.levels) {
+ level.interpolate();
+ }
+ }
+ }
+ /*
+ * Finished iteration on a row. Clear flags and update position
+ * before to move to next row.
+ */
+ for (int b=0; b<numBands; b++) {
+ for (final IsolineTracer.Level level : isolines[b].levels) {
+ level.finishedRow();
+ }
+ }
+ tracer.x = 0;
+ tracer.y++;
+ }
+ /*
+ * Finished iteration over the whole image.
+ */
+ for (int b=0; b<numBands; b++) {
+ for (final IsolineTracer.Level level : isolines[b].levels) {
+ level.finish();
+ }
+ }
+ return isolines;
+ }
+
+ /**
+ * Returns the polylines for each level specified to the
+ * {@link #generate(RenderedImage, double[][], double, MathTransform)} method.
+ *
+ * @return the polylines for each level.
+ */
+ public final NavigableMap<Double,Shape> polylines() {
+ final TreeMap<Double,Shape> paths = new TreeMap<>();
+ for (final IsolineTracer.Level level : levels) {
+ final Path2D path = level.path;
+ if (path != null) {
+ // TODO: invoke path.trimToSize() with JDK10.
+ paths.put(level.value, path);
+ }
+ }
+ return paths;
+ }
+}
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/package-info.java b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/package-info.java
new file mode 100644
index 0000000..a087331
--- /dev/null
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/package-info.java
@@ -0,0 +1,32 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/**
+ * Image processing operations.
+ *
+ * <p><strong>Do not use!</strong></p>
+ *
+ * This package is for internal use by SIS only. Classes in this package
+ * may change in incompatible ways in any future version without notice.
+ *
+ * @author Johann Sorel (Geomatys)
+ * @author Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since 1.1
+ * @module
+ */
+package org.apache.sis.internal.processing.image;
diff --git a/core/sis-feature/src/test/java/org/apache/sis/internal/processing/image/IsolinesTest.java b/core/sis-feature/src/test/java/org/apache/sis/internal/processing/image/IsolinesTest.java
new file mode 100644
index 0000000..d0aaf2d
--- /dev/null
+++ b/core/sis-feature/src/test/java/org/apache/sis/internal/processing/image/IsolinesTest.java
@@ -0,0 +1,308 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.internal.processing.image;
+
+import java.util.Map;
+import java.awt.Shape;
+import java.awt.geom.PathIterator;
+import java.awt.image.DataBuffer;
+import java.awt.image.BufferedImage;
+import java.awt.image.WritableRaster;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.internal.coverage.j2d.RasterFactory;
+import org.apache.sis.test.TestCase;
+import org.junit.Test;
+
+import static org.junit.Assert.*;
+
+
+/**
+ * Tests {@link Isolines}.
+ *
+ * @author Johann Sorel (Geomatys)
+ * @author Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since 1.1
+ * @module
+ */
+public final strictfp class IsolinesTest extends TestCase {
+ /**
+ * Tolerance threshold for rounding errors.
+ * Could be zero, but we nevertheless keep a margin for safety.
+ */
+ private static final double TOLERANCE = 1E-14;
+
+ /**
+ * The threshold value. This is used by {@code generateFromXXX(…)} methods.
+ */
+ private double threshold;
+
+ /**
+ * The isoline being tested. This is set by {@code generateFromXXX(…)} methods.
+ */
+ private Shape isoline;
+
+ /**
+ * Tests isolines computed in a contouring grid having only one cell.
+ * The cell may have zero, one or two line segments.
+ *
+ * @throws TransformException if a point can not be transformed to its final coordinate space.
+ */
+ @Test
+ public void testSingleCell() throws TransformException {
+ threshold = -5;
+ generateFromCell(0, 0, 0, 0);
+ assertNull(isoline);
+
+ threshold = 2;
+ generateFromCell(0, 0, 0, 0);
+ assertNull(isoline);
+ /*
+ * ▘5╌╌╌╌1
+ * ╎ ╎
+ * 0╌╌╌╌0
+ */
+ generateFromCell(5, 1, 0, 0);
+ assertSegmentEquals(0, 0.6, 0.75, 0);
+ /*
+ * 0╌╌╌╌5▝
+ * ╎ ╎
+ * ▖10╌╌20▗
+ */
+ generateFromCell(0, 5, 10, 20);
+ assertSegmentEquals(0, 0.2, 0.4, 0);
+ /*
+ * 1╌╌╌╌5▝
+ * ╎ ╎
+ * 0╌╌╌╌0
+ */
+ generateFromCell(1, 5, 0, 0);
+ assertSegmentEquals(0.25, 0, 1, 0.6);
+ /*
+ * ▘5╌╌╌╌0
+ * ╎ ╎
+ * ▖20╌╌10▗
+ */
+ generateFromCell(5, 0, 20, 10);
+ assertSegmentEquals(0.6, 0, 1, 0.2);
+ /*
+ * 1╌╌╌╌0
+ * ╎ ╎
+ * ▖5╌╌╌10▗
+ */
+ generateFromCell(1, 0, 5, 10);
+ assertSegmentEquals(0, 0.25, 1, 0.2);
+ /*
+ * ▘5╌╌╌10▝
+ * ╎ ╎
+ * 1╌╌╌╌0
+ */
+ generateFromCell(5, 10, 1, 0);
+ assertSegmentEquals(0, 0.75, 1, 0.8);
+ /*
+ * 0╌╌╌╌0
+ * ╎ ╎
+ * 1╌╌╌╌5▗
+ */
+ generateFromCell(0, 0, 1, 5);
+ assertSegmentEquals(1, 0.4, 0.25, 1);
+ /*
+ * ▘20╌╌10▝
+ * ╎ ╎
+ * ▖5╌╌╌╌0
+ */
+ generateFromCell(20, 10, 5, 0);
+ assertSegmentEquals(1, 0.8, 0.6, 1);
+ /*
+ * 1╌╌╌11▝
+ * ╎ ╎
+ * ▖5╌╌╌╌1
+ */
+ generateFromCell(1, 11, 5, 1);
+ assertSegmentsEqual(0, 0.25, 0.1, 0,
+ 1, 0.9, 0.75, 1);
+ /*
+ * ▘11╌╌╌1
+ * ╎ ╎
+ * 1╌╌╌╌5▗
+ */
+ generateFromCell(11, 1, 1, 5);
+ assertSegmentsEqual(0.9, 0, 1, 0.25,
+ 0, 0.9, 0.25, 1);
+ /*
+ * ▘5╌╌╌╌1
+ * ╎ ╎
+ * ▖10╌╌╌0
+ */
+ generateFromCell(5, 1, 10, 0);
+ assertSegmentEquals(0.75, 0, 0.8, 1);
+ /*
+ * 1╌╌╌╌5▝
+ * ╎ ╎
+ * 0╌╌╌10▗
+ */
+ generateFromCell(1, 5, 0, 10);
+ assertSegmentEquals(0.25, 0, 0.2, 1);
+ /*
+ * 0╌╌╌╌0
+ * ╎ ╎
+ * ▖5╌╌╌╌1
+ */
+ generateFromCell(0, 0, 5, 1);
+ assertSegmentEquals(0, 0.4, 0.75, 1);
+ /*
+ * ▘10╌╌20▝
+ * ╎ ╎
+ * 0╌╌╌╌5▗
+ */
+ generateFromCell(10, 20, 0, 5);
+ assertSegmentEquals(0, 0.8, 0.4, 1);
+ }
+
+ /**
+ * Tests isolines computed in a contouring grid having 2×2 cells.
+ *
+ * @throws TransformException if a point can not be transformed to its final coordinate space.
+ */
+ @Test
+ public void testMultiCells() throws TransformException {
+ threshold = 0.5;
+ generateFromImage(3,
+ 0,0,0,
+ 0,1,0,
+ 0,0,0);
+ /*
+ * Expected coordinates:
+ *
+ * (1): 1.5 1.0 (2)
+ * (2): 1.0 0.5 (3) (1)
+ * (3): 0.5 1.0 (4)
+ * (4): 1.0 1.5
+ */
+ final double[] buffer = new double[2];
+ final PathIterator it = isoline.getPathIterator(null);
+ assertSegmentEquals(it, buffer, 1.5, 1, 1, 0.5);
+
+ assertFalse(it.isDone());
+ assertEquals(PathIterator.SEG_LINETO, it.currentSegment(buffer));
+ assertEquals("x2", 0.5, buffer[0], TOLERANCE);
+ assertEquals("y2", 1, buffer[1], TOLERANCE);
+
+ it.next();
+ assertFalse(it.isDone());
+ assertEquals(PathIterator.SEG_LINETO, it.currentSegment(buffer));
+ assertEquals("x3", 1, buffer[0], TOLERANCE);
+ assertEquals("y3", 1.5, buffer[1], TOLERANCE);
+
+// it.next();
+// assertFalse(it.isDone());
+// assertEquals(PathIterator.SEG_CLOSE, it.currentSegment(buffer));
+
+ it.next();
+ assertTrue(it.isDone());
+ }
+
+ /**
+ * Tests a cell containing a NaN value.
+ *
+ * @throws TransformException if a point can not be transformed to its final coordinate space.
+ */
+ @Test
+ public void testNaN() throws TransformException {
+ threshold = 2;
+ /*
+ * ▘5╌╌NaN▝
+ * ╎ ╎
+ * 0╌╌╌╌0
+ */
+ generateFromCell(5, Float.NaN, 0, 0);
+ assertNull(isoline);
+ }
+
+ /**
+ * Generates isolines from a 2×2 image having the given values.
+ * The result is stored in {@link #isoline}; it may be {@code null}.
+ *
+ * @throws TransformException if a point can not be transformed to its final coordinate space.
+ */
+ private void generateFromCell(float v00, float v10, float v01, float v11) throws TransformException {
+ generateFromImage(2, v00, v10, v01, v11);
+ }
+
+ /**
+ * Generates isolines from a size×size image having the given values.
+ * The result is stored in {@link #isoline}; it may be {@code null}.
+ *
+ * @throws TransformException if a point can not be transformed to its final coordinate space.
+ */
+ private void generateFromImage(final int size, final float... values) throws TransformException {
+ final BufferedImage image = RasterFactory.createGrayScaleImage(DataBuffer.TYPE_FLOAT, size, size, 1, 0, 0, 10);
+ final WritableRaster raster = image.getRaster();
+ for (int i=0; i<values.length; i++) {
+ raster.setSample(i % size, i / size, 0, values[i]);
+ }
+ final Isolines[] isolines = Isolines.generate(image, new double[][] {{threshold}}, 0, null);
+ assertEquals("Number of bands", 1, isolines.length);
+ final Map<Double, Shape> polylines = isolines[0].polylines();
+ assertTrue(polylines.size() <= 1);
+ isoline = polylines.get(threshold);
+ }
+
+ /**
+ * Asserts that {@link #isoline} is a segment having the given coordinates.
+ */
+ private void assertSegmentEquals(final double x0, final double y0, final double x1, final double y1) {
+ assertNotNull("isoline", isoline);
+ final double[] buffer = new double[2];
+ final PathIterator it = isoline.getPathIterator(null);
+ assertSegmentEquals(it, buffer, x0, y0, x1, y1);
+ assertTrue(it.isDone());
+ }
+
+ /**
+ * Asserts that the iterator contains a {@code SEG_MOVETO} followed by a {@code SEG_LINETO} instruction
+ * with the given coordinates. The segment is positioned on next segment after this method call.
+ */
+ private static void assertSegmentEquals(final PathIterator it, final double[] buffer,
+ final double x0, final double y0, final double x1, final double y1)
+ {
+ assertFalse(it.isDone());
+ assertEquals(PathIterator.SEG_MOVETO, it.currentSegment(buffer));
+ assertEquals("x0", x0, buffer[0], TOLERANCE);
+ assertEquals("y0", y0, buffer[1], TOLERANCE);
+ it.next();
+ assertFalse(it.isDone());
+ assertEquals(PathIterator.SEG_LINETO, it.currentSegment(buffer));
+ assertEquals("x1", x1, buffer[0], TOLERANCE);
+ assertEquals("y1", y1, buffer[1], TOLERANCE);
+ it.next();
+ }
+
+ /**
+ * Asserts that {@link #isoline} is two segments having the given coordinates.
+ */
+ private void assertSegmentsEqual(final double x0, final double y0, final double x1, final double y1,
+ final double x2, final double y2, final double x3, final double y3)
+ {
+ assertNotNull("isoline", isoline);
+ final double[] buffer = new double[2];
+ final PathIterator it = isoline.getPathIterator(null);
+ assertSegmentEquals(it, buffer, x0, y0, x1, y1);
+ assertSegmentEquals(it, buffer, x2, y2, x3, y3);
+ assertTrue(it.isDone());
+ }
+}
diff --git a/core/sis-feature/src/test/java/org/apache/sis/test/suite/FeatureTestSuite.java b/core/sis-feature/src/test/java/org/apache/sis/test/suite/FeatureTestSuite.java
index fa96385..67cea45 100644
--- a/core/sis-feature/src/test/java/org/apache/sis/test/suite/FeatureTestSuite.java
+++ b/core/sis-feature/src/test/java/org/apache/sis/test/suite/FeatureTestSuite.java
@@ -108,7 +108,8 @@ import org.junit.runners.Suite;
org.apache.sis.coverage.grid.ConvertedGridCoverageTest.class,
org.apache.sis.coverage.grid.ResampledGridCoverageTest.class,
- // Index
+ // Index and processing
+ org.apache.sis.internal.processing.image.IsolinesTest.class,
org.apache.sis.index.tree.PointTreeNodeTest.class,
org.apache.sis.index.tree.PointTreeTest.class
})