sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [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.
Date Fri, 11 Dec 2020 17:37:30 GMT
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 <martin.desruisseaux@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
 })


Mime
View raw message