Merge fixes from branch 'geoapi-4.0' into geoapi-3.1;

* GML and WKT parsing fixes.
* https://issues.apache.org/jira/browse/SIS-555
* https://issues.apache.org/jira/browse/SIS-166 partially reverted.
diff --git a/application/sis-javafx/src/main/java/org/apache/sis/gui/coverage/CoverageCanvas.java b/application/sis-javafx/src/main/java/org/apache/sis/gui/coverage/CoverageCanvas.java
index daf3850..eab02b0 100644
--- a/application/sis-javafx/src/main/java/org/apache/sis/gui/coverage/CoverageCanvas.java
+++ b/application/sis-javafx/src/main/java/org/apache/sis/gui/coverage/CoverageCanvas.java
@@ -70,7 +70,7 @@
 import org.apache.sis.portrayal.RenderException;
 import org.apache.sis.internal.map.coverage.RenderingWorkaround;
 import org.apache.sis.internal.coverage.j2d.TileErrorHandler;
-import org.apache.sis.internal.processing.image.Isolines;
+import org.apache.sis.internal.processing.isoline.Isolines;
 import org.apache.sis.internal.gui.BackgroundThreads;
 import org.apache.sis.internal.gui.ExceptionReporter;
 import org.apache.sis.internal.gui.GUIUtilities;
diff --git a/application/sis-javafx/src/main/java/org/apache/sis/gui/coverage/IsolineRenderer.java b/application/sis-javafx/src/main/java/org/apache/sis/gui/coverage/IsolineRenderer.java
index 6391717..d816ab6 100644
--- a/application/sis-javafx/src/main/java/org/apache/sis/gui/coverage/IsolineRenderer.java
+++ b/application/sis-javafx/src/main/java/org/apache/sis/gui/coverage/IsolineRenderer.java
@@ -38,7 +38,7 @@
 import org.opengis.referencing.operation.MathTransform;
 import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.internal.gui.control.ColorRamp;
-import org.apache.sis.internal.processing.image.Isolines;
+import org.apache.sis.internal.processing.isoline.Isolines;
 import org.apache.sis.internal.coverage.j2d.ImageUtilities;
 import org.apache.sis.internal.gui.control.ValueColorMapper.Step;
 import org.apache.sis.internal.feature.j2d.EmptyShape;
diff --git a/application/sis-javafx/src/main/java/org/apache/sis/gui/coverage/StyledRenderingData.java b/application/sis-javafx/src/main/java/org/apache/sis/gui/coverage/StyledRenderingData.java
index 0709fe8..ed19c7d 100644
--- a/application/sis-javafx/src/main/java/org/apache/sis/gui/coverage/StyledRenderingData.java
+++ b/application/sis-javafx/src/main/java/org/apache/sis/gui/coverage/StyledRenderingData.java
@@ -24,7 +24,7 @@
 import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.storage.DataStoreException;
 import org.apache.sis.image.ErrorHandler;
-import org.apache.sis.internal.processing.image.Isolines;
+import org.apache.sis.internal.processing.isoline.Isolines;
 import org.apache.sis.internal.map.coverage.RenderingData;
 
 
diff --git a/core/sis-feature/src/main/java/org/apache/sis/image/ImageProcessor.java b/core/sis-feature/src/main/java/org/apache/sis/image/ImageProcessor.java
index ef4923f..23575b6 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/image/ImageProcessor.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/image/ImageProcessor.java
@@ -51,7 +51,7 @@
 import org.apache.sis.util.collection.WeakHashSet;
 import org.apache.sis.internal.system.Modules;
 import org.apache.sis.internal.coverage.j2d.TiledImage;
-import org.apache.sis.internal.processing.image.Isolines;
+import org.apache.sis.internal.processing.isoline.Isolines;
 import org.apache.sis.internal.feature.Resources;
 import org.apache.sis.measure.NumberRange;
 import org.apache.sis.measure.Units;
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/MultiPolylines.java b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/MultiPolylines.java
index aea1b7a..1139819 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/MultiPolylines.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/MultiPolylines.java
@@ -23,6 +23,7 @@
 import java.awt.geom.PathIterator;
 import java.awt.geom.AffineTransform;
 import org.apache.sis.internal.referencing.j2d.IntervalRectangle;
+import org.apache.sis.util.Classes;
 
 
 /**
@@ -36,7 +37,7 @@
  * </ul>
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.1
+ * @version 1.3
  * @since   1.1
  * @module
  */
@@ -178,4 +179,14 @@
         }
         return null;
     }
+
+    /**
+     * Returns a string representation for debugging purposes.
+     *
+     * @return a debug string representation.
+     */
+    @Override
+    public String toString() {
+        return Classes.getShortClassName(this) + '[' + (polylines.length / 2) + " polylines]";
+    }
 }
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/PathBuilder.java b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/PathBuilder.java
index 3977234..f4ae2b2 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/PathBuilder.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/PathBuilder.java
@@ -221,7 +221,7 @@
      * The {@link #createPolyline(boolean)} method should be invoked before this method
      * for making sure that there are no pending polylines.
      *
-     * @return the polyline, polygon or collector of polylines.
+     * @return the polyline, polygon or collection of polylines.
      *         May be {@code null} if no polyline or polygon has been created.
      */
     public final Shape build() {
@@ -233,6 +233,17 @@
     }
 
     /**
+     * Returns a snapshot of currently added polylines or polygons without modifying the state of this builder.
+     * It is safe to continue building the shape and invoke this method again later for progressive rendering.
+     *
+     * @return the polyline, polygon or collection of polylines added so far.
+     *         May be {@code null} if no polyline or polygon has been created.
+     */
+    public final Shape snapshot() {
+        return build();
+    }
+
+    /**
      * Returns a string representation of the polyline under construction for debugging purposes.
      * Current implementation formats only the first and last points, and tells how many points are between.
      */
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/Polyline.java b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/Polyline.java
index 2272236..5b9c8ac 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/Polyline.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/Polyline.java
@@ -24,6 +24,7 @@
 import java.awt.geom.PathIterator;
 import java.awt.geom.AffineTransform;
 import org.apache.sis.internal.util.Numerics;
+import org.apache.sis.util.Classes;
 
 
 /**
@@ -57,7 +58,7 @@
  * </ol>
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.1
+ * @version 1.3
  * @since   1.1
  * @module
  */
@@ -304,4 +305,14 @@
             return (position == 0) ? SEG_MOVETO : SEG_LINETO;
         }
     }
+
+    /**
+     * Returns a string representation for debugging purposes.
+     *
+     * @return a debug string representation.
+     */
+    @Override
+    public String toString() {
+        return Classes.getShortClassName(this) + '[' + (coordinates.length / 2) + " points]";
+    }
 }
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/package-info.java b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/package-info.java
index d7ee3c9..f889486 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/package-info.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/j2d/package-info.java
@@ -25,7 +25,7 @@
  * may change in incompatible ways in any future version without notice.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.2
+ * @version 1.3
  *
  * @see org.apache.sis.internal.referencing.j2d
  *
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
deleted file mode 100644
index 63dcb38..0000000
--- a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/IsolineTracer.java
+++ /dev/null
@@ -1,1261 +0,0 @@
-/*
- * 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.ArrayList;
-import java.util.Map;
-import java.util.HashMap;
-import java.util.IdentityHashMap;
-import java.util.Collections;
-import java.awt.Point;
-import java.awt.Rectangle;
-import java.awt.Shape;
-import org.opengis.referencing.operation.MathTransform;
-import org.opengis.referencing.operation.TransformException;
-import org.apache.sis.internal.feature.j2d.PathBuilder;
-import org.apache.sis.internal.util.Numerics;
-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.
-     * The length of this array is <var>(number of bands)</var> × 2 (width) × 2 (height).
-     */
-    private final double[] 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;
-
-    /**
-     * Translation to apply on coordinates. For isolines computed sequentially, this is the image origin
-     * (often 0,0 but not necessarily). For isolines computed in parallel, the translations are different
-     * for each computation tile.
-     */
-    private final double translateX, translateY;
-
-    /**
-     * Final transform to apply on coordinates (integer source coordinates at pixel centers).
-     * Can be {@code null} if none.
-     */
-    private final MathTransform gridToCRS;
-
-    /**
-     * Creates a new position for the given data window.
-     *
-     * @param  window       the 2×2 window containing pixel values in the 4 corners of current contouring grid cell.
-     * @param  pixelStride  increment to the position in {@code window} for reading next sample value.
-     * @param  domain       pixel coordinates where iteration will happen.
-     * @param  gridToCRS    final transform to apply on coordinates (integer source coordinates at pixel centers).
-     */
-    IsolineTracer(final double[] window, final int pixelStride, final Rectangle domain, final MathTransform gridToCRS) {
-        this.window      = window;
-        this.pixelStride = pixelStride;
-        this.translateX  = domain.x;
-        this.translateY  = domain.y;
-        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 {
-        /**
-         * Band number where to read values in the {@link #window} array.
-         */
-        private final int band;
-
-        /**
-         * The level value.
-         *
-         * @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 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 are no polylines 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;
-
-        /**
-         * Paths that have not yet been closed. The {@link Polyline} coordinates are copied in this map when iteration
-         * finished on a row and the polyline is not reused by next row, or when the {@link #closeLeftWithTop(Polyline)}
-         * method has been invoked but the geometry to close is still not complete. This map accumulates those partial
-         * shapes for assembling them later when missing parts become available.
-         *
-         * <h4>Map keys</h4>
-         * Keys are grid coordinates rounded toward 0. The coordinate having fraction digits has its bits inverted
-         * by the {@code ~} operator. For each point, there is at most one coordinate having such fraction digits.
-         *
-         * <h4>Map values</h4>
-         * {@code Unclosed} instances are list of {@code double[]} arrays to be concatenated in a single polygon later.
-         * For a given {@code Unclosed} list, all {@code double[]} arrays at even indices shall have their points read
-         * in reverse order and all {@code double[]} arrays at odd indices shall have their points read in forward order.
-         * The list may contain null elements when there is no data in the corresponding iteration order.
-         *
-         * @see #closeLeftWithTop(Polyline)
-         */
-        private final Map<Point,Unclosed> partialPaths;
-
-        /**
-         * Builder of isolines as a Java2D shape, created when first needed.
-         * The {@link Polyline} coordinates are copied in this path when a geometry is closed.
-         *
-         * @see #writeTo(Joiner, Polyline[], boolean)
-         */
-        private Joiner path;
-
-        /**
-         * The isolines as a Java2D shape, created by {@link #finish()}.
-         * This is the shape to be returned to user for this level after we finished to process all cells.
-         */
-        Shape shape;
-
-        /**
-         * Creates new isoline levels for the given value.
-         *
-         * @param  band   band number where to read values in the {@link #window} array.
-         * @param  value  the isoline level value.
-         * @param  width  the contouring grid cell width (one cell smaller than image width).
-         */
-        Level(final int band, final double value, final int width) {
-            this.band      = band;
-            this.value     = value;
-            partialPaths   = new HashMap<>();
-            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.
-         *
-         * @see Isolines#setMaskBit(double, int)
-         */
-        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
-         * the {@link PathBuilder}.
-         */
-        @SuppressWarnings("AssertWithSideEffects")
-        final void interpolate() throws TransformException {
-            /*
-             * Note: `interpolateMissingLeftSide()` and `interpolateMissingTopSide(…)` should do interpolations
-             * only for cells in the first column and first row respectively. We could avoid those method calls
-             * for all other cells if we add two flags in the `isDataAbove` bitmask: FIRST_ROW and FIRST_COLUMN.
-             * The switch cases then become something like below:
-             *
-             *     case <bitmask> | FIRST_COLUMN | FIRST_ROW:
-             *     case <bitmask> | FIRST_COLUMN: {
-             *         interpolateMissingLeftSide();
-             *         // Fall through
-             *     }
-             *     case <bitmask> | FIRST_ROW:
-             *     case <bitmask>: {
-             *         // Interpolations on other borders.
-             *         break;
-             *     }
-             *
-             * We tried that approach, but benchmarking on Java 15 suggested a small performance decrease
-             * instead of an improvement. It may be worth to try again in the future, after advancement
-             * in compiler technology.
-             */
-            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();
-                    interpolateMissingLeftSide();
-                    interpolateOnRightSide();                   // 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];
-                    interpolateMissingTopSide(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();
-                    interpolateMissingLeftSide();
-                    interpolateOnBottomSide(polylinesOnTop[x].transferFrom(polylineOnLeft));
-                    break;
-                }
-                /*     ○╌╌╌╌╲╌●        ●╌╌╌╌╲╌○
-                 *     ╎     ╲╎        ╎     ╲╎
-                 *     ╎      ╲        ╎      ╲
-                 *     ○╌╌╌╌╌╌○╲       ●╌╌╌╌╌╌●╲
-                 */
-                case UPPER_RIGHT:
-                case UPPER_LEFT | LOWER_LEFT | LOWER_RIGHT: {
-                    assert polylineOnLeft.isEmpty();
-                    interpolateMissingTopSide(polylineOnLeft.transferFrom(polylinesOnTop[x]));
-                    interpolateOnRightSide();
-                    break;
-                }
-                /*     ○╌╌╌╌╌╌○╱       ●╌╌╌╌╌╌●╱
-                 *     ╎      ╱        ╎      ╱
-                 *     ╎     ╱╎        ╎     ╱╎
-                 *     ○╌╌╌╌╱╌●        ●╌╌╌╌╱╌○
-                 */
-                case LOWER_RIGHT:
-                case UPPER_LEFT | UPPER_RIGHT | LOWER_LEFT: {
-                    assert polylinesOnTop[x].isEmpty();
-                    assert polylineOnLeft   .isEmpty();
-                    interpolateOnRightSide();
-                    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: {
-                    closeLeftWithTop(polylinesOnTop[x]);
-                    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 double[] data = window;
-                        int p = band;
-                        do average += data[p];
-                        while ((p += pixelStride) < data.length);
-                        assert (p -= band) == pixelStride * 4 : p;
-                        average /= 4;
-                    }
-                    boolean LLtoUR = isDataAbove == (LOWER_LEFT | UPPER_RIGHT);
-                    LLtoUR ^= (average <= value);
-                    final Polyline polylineOnTop = polylinesOnTop[x];
-                    if (LLtoUR) {
-                        closeLeftWithTop(polylineOnTop);
-                        interpolateOnRightSide();
-                        interpolateOnBottomSide(polylineOnTop.attach(polylineOnLeft));
-                    } else {
-                        interpolateMissingLeftSide();
-                        final Polyline swap = new Polyline().transferFrom(polylineOnTop);
-                        interpolateOnBottomSide(polylineOnTop.transferFrom(polylineOnLeft));
-                        interpolateMissingTopSide(polylineOnLeft.transferFrom(swap));
-                        interpolateOnRightSide();
-                    }
-                    break;
-                }
-            }
-        }
-
-        /**
-         * Appends to {@link #polylineOnLeft} a point interpolated on the left side if that point is missing.
-         * This interpolation should happens only in the first column.
-         */
-        private void interpolateMissingLeftSide() {
-            if (polylineOnLeft.size == 0) {
-                polylineOnLeft.append(translateX + (x),
-                                      translateY + (y + interpolate(0, 2*pixelStride)));
-            }
-        }
-
-        /**
-         * Appends to {@code polylineOnTop} a point interpolated on the top side if that point is missing.
-         * This interpolation should happens only in the first row.
-         */
-        private void interpolateMissingTopSide(final Polyline polylineOnTop) {
-            if (polylineOnTop.size == 0) {
-                interpolateOnTopSide(polylineOnTop);
-            }
-        }
-
-        /**
-         * Appends to the given polyline a point interpolated on the top side.
-         */
-        private void interpolateOnTopSide(final Polyline appendTo) {
-            appendTo.append(translateX + (x + interpolate(0, pixelStride)),
-                            translateY + (y));
-        }
-
-        /**
-         * Appends to {@link #polylineOnLeft} a point interpolated on the right side.
-         * The polyline on right side will become {@code polylineOnLeft} in next column.
-         */
-        private void interpolateOnRightSide() {
-            polylineOnLeft.append(translateX + (x + 1),
-                                  translateY + (y + interpolate(pixelStride, 3*pixelStride)));
-        }
-
-        /**
-         * Appends to the given polyline a point interpolated on the bottom side.
-         * The polyline on top side will become a {@code polylineOnBottoù} in next row.
-         */
-        private void interpolateOnBottomSide(final Polyline polylineOnTop) {
-            polylineOnTop.append(translateX + (x + interpolate(2*pixelStride, 3*pixelStride)),
-                                 translateY + (y + 1));
-        }
-
-        /**
-         * Interpolates the position where the isoline passes between two values.
-         *
-         * @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 double[] data = window;
-            final int    p  = band;
-            final double v1 = data[p + i1];
-            final double v2 = data[p + i2];
-            return (value - v1) / (v2 - v1);
-        }
-
-        /**
-         * Joins {@link #polylineOnLeft} with {@code polylineOnTop}, saves their coordinates and clear
-         * those {@code Polyline} instances for use in next cell. The coordinates are written directly
-         * to {@link #path} if we got a closed polygon, or otherwise are saved in {@link #partialPaths}
-         * for later processing. This method is invoked for cells like below:
-         *
-         * {@preformat text
-         *     ●╌╱╌╌╌╌○        ○╌╱╌╌╌╌●        ○╌╱╌╌╌╌●╱
-         *     ╎╱     ╎        ╎╱     ╎        ╎╱     ╱
-         *     ╱      ╎        ╱      ╎        ╱     ╱╎
-         *    ╱○╌╌╌╌╌╌○       ╱●╌╌╌╌╌╌●       ╱●╌╌╌╌╱╌○
-         * }
-         *
-         * This method does itself the interpolations on left side and top side. The two polylines
-         * {@link #polylineOnLeft} and {@code polylineOnTop} will become empty after this method call.
-         *
-         * @param  polylineOnTop  value of {@code polylinesOnTop[x]}.
-         * @throws TransformException if the {@link IsolineTracer#gridToCRS} transform can not be applied.
-         */
-        private void closeLeftWithTop(final Polyline polylineOnTop) throws TransformException {
-            interpolateMissingLeftSide();
-            interpolateMissingTopSide(polylineOnTop);
-            final Polyline[] polylines;
-            if (polylineOnLeft.opposite == polylineOnTop) {
-                assert polylineOnTop.opposite == polylineOnLeft;
-                /*
-                 * We have a loop: the polygon can be closed now, without copying coordinates to temporary buffers.
-                 * Points in the two `Polyline` instances will be iterated in (reverse, forward) order respectively.
-                 * Consequently the points we just interpolated will be first point and last point before closing.
-                 */
-                polylines = new Polyline[] {polylineOnTop, polylineOnLeft};     // (reverse, forward) point order.
-            } else {
-                /*
-                 * Joining left and top polylines do not yet create a closed shape. Consequently we may not write
-                 * in the `path` now. But maybe we can close the polygon later after more polylines are attached.
-                 */
-                final Unclosed fragment = new Unclosed(polylineOnLeft, polylineOnTop);
-                if (fragment.isEmpty()) {
-                    /*
-                     * Fragment starts and ends with NaN values. We will not be able to complete a polygon.
-                     * Better to write the polylines now for avoiding temporary copies of their coordinates.
-                     */
-                    polylines = new Polyline[] {
-                        polylineOnLeft.opposite, polylineOnLeft, polylineOnTop, polylineOnTop.opposite
-                    };
-                } else if (fragment.addOrMerge(partialPaths)) {
-                    /*
-                     * The fragment has been merged with previously existing fragments and became a polygon.
-                     * We can write the polygon immediately. There are no more references to those coordinates
-                     * in the `partialPaths` map.
-                     */
-                    polylines = fragment.toPolylines();
-                } else {
-                    return;
-                }
-            }
-            path = writeTo(path, polylines, true);
-        }
-
-        /**
-         * Writes the content of given polyline without closing it as a polygon.
-         * The given polyline will become empty after this method call.
-         */
-        private void writeUnclosed(final Polyline polyline) throws TransformException {
-            final Unclosed fragment = new Unclosed(polyline, null);
-            final Polyline[] polylines;
-            final boolean close;
-            if (fragment.isEmpty()) {
-                close = false;
-                polylines = new Polyline[] {polyline.opposite, polyline};       // (reverse, forward) point order.
-            } else {
-                close = fragment.addOrMerge(partialPaths);
-                if (!close) {
-                    // Keep in `partialPaths`. Maybe it can be closed later.
-                    return;
-                }
-                polylines = fragment.toPolylines();
-            }
-            path = writeTo(path, polylines, close);
-        }
-
-        /**
-         * Invoked after iteration on a single row has been completed. If there is a polyline
-         * finishing on the right image border, the coordinates needs to be saved somewhere
-         * because that {@code Polyline} will not be continued by cells on next rows.
-         */
-        final void finishedRow() throws TransformException {
-            if (!polylineOnLeft.transferToOpposite()) {
-                writeUnclosed(polylineOnLeft);
-            }
-            isDataAbove = 0;
-        }
-
-        /**
-         * Invoked after the iteration has been completed on the full area of interest.
-         * This method writes all remaining polylines to {@link #partialPaths}.
-         * It assumes that {@link #finishedRow()} has already been invoked.
-         * This {@code Isoline} can not be used anymore after this call.
-         */
-        final void finish() throws TransformException {
-            assert polylineOnLeft.isEmpty();
-            polylineOnLeft.coordinates = null;
-            /*
-             * This method sets various values to null for letting the garbage collector do its work.
-             * This is okay for a `Level` instance which is not going to be used anymore, except for
-             * reading the `shape` field.
-             */
-            for (int i=0; i < polylinesOnTop.length; i++) {
-                writeUnclosed(polylinesOnTop[i]);
-                polylinesOnTop[i] = null;
-            }
-            assert isConsistent();
-        }
-
-        /**
-         * Verifies that {@link #partialPaths} consistency. Used for assertions only.
-         */
-        private boolean isConsistent() {
-            for (final Map.Entry<Point,Unclosed> entry : partialPaths.entrySet()) {
-                if (!entry.getValue().isExtremity(entry.getKey())) return false;
-            }
-            return true;
-        }
-
-        /**
-         * Transfers all {@code other} polylines into this instance. The {@code other} instance should be a neighbor,
-         * i.e. an instance sharing a border with this instance. The {@code other} instance will become empty after
-         * this method call.
-         *
-         * @param  other  a neighbor level (on top, left, right or bottom) to merge with this level.
-         * @throws TransformException if an error occurred during polylines creation.
-         */
-        final void merge(final Level other) throws TransformException {
-            assert other != this && other.value == value;
-            if (path == null) {
-                path = other.path;
-            } else {
-                path.append(other.path);
-            }
-            other.path = null;
-            assert  this.isConsistent();
-            assert other.isConsistent();
-            final IdentityHashMap<Unclosed,Boolean> done = new IdentityHashMap<>(other.partialPaths.size() / 2);
-            for (final Map.Entry<Point,Unclosed> entry : other.partialPaths.entrySet()) {
-                final Unclosed fragment = entry.getValue();
-                if (done.put(fragment, Boolean.TRUE) == null) {
-                    assert fragment.isExtremity(entry.getKey());
-                    if (fragment.addOrMerge(partialPaths)) {
-                        path = writeTo(path, fragment.toPolylines(), true);
-                        fragment.clear();
-                    }
-                }
-                entry.setValue(null);       // Let the garbage collector do its work.
-            }
-        }
-
-        /**
-         * Flushes any pending {@link #partialPaths} to {@link #path}. This method is invoked after
-         * {@link #finish()} has been invoked for all sub-regions (many sub-regions may exist if
-         * isoline generation has been splitted for parallel computation).
-         *
-         * @throws TransformException if an error occurred during polylines creation.
-         */
-        final void flush() throws TransformException {
-            for (final Map.Entry<Point,Unclosed> entry : partialPaths.entrySet()) {
-                final Unclosed fragment = entry.getValue();
-                assert fragment.isExtremity(entry.getKey());
-                if (!fragment.isEmpty()) {
-                    path = writeTo(path, fragment.toPolylines(), false);
-                    fragment.clear();       // Necessary because the same list appears twice in the map.
-                }
-                entry.setValue(null);       // Let the garbage collector do its work.
-            }
-            if (path != null) {
-                shape = path.build();
-                path  = 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.
-         */
-        static final int DIMENSION = 2;
-
-        /**
-         * Coordinates as (x,y) tuples. This array is expanded as needed.
-         */
-        double[] coordinates;
-
-        /**
-         * Number of valid elements in the {@link #coordinates} array.
-         * This is twice the number of points.
-         */
-        int size;
-
-        /**
-         * If the polyline has points added to its two extremities, the other extremity. Otherwise {@code null}.
-         * The first point of {@code opposite} polyline is connected to the first point of this polyline.
-         * Consequently when those two polylines are joined in a single polyline, the coordinates of either
-         * {@code this} or {@code opposite} must be iterated in reverse order.
-         */
-        Polyline opposite;
-
-        /**
-         * Creates an initially empty polyline.
-         */
-        Polyline() {
-            coordinates = ArraysExt.EMPTY_DOUBLE;
-        }
-
-        /**
-         * Creates a new polyline wrapping the given coordinates. Used for wrapping {@link Unclosed}
-         * instances in objects expected by {@link IsolineTracer#writeTo(Joiner, Polyline[], boolean)}.
-         * Those {@code Polyline} instances are temporary.
-         */
-        Polyline(final double[] data) {
-            coordinates = data;
-            size = data.length;
-        }
-
-        /**
-         * Discards all coordinates in this polyline. This method does not clear
-         * the {@link #opposite} polyline; it is caller's responsibility to do so.
-         */
-        final void clear() {
-            opposite = null;
-            size = 0;
-        }
-
-        /**
-         * Returns whether this polyline is empty. This method is used only for {@code assert isEmpty()}
-         * statement because of the check for {@code opposite == null}: an empty polyline should not have
-         * a non-null {@link #opposite} value.
-         */
-        final boolean isEmpty() {
-            return size == 0 & (opposite == null);
-        }
-
-        /**
-         * 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).
-         *
-         * <p>This method is typically invoked in the following pattern (but this is not mandatory).
-         * An important aspect is that {@code this} and {@code other} should be on perpendicular axes:</p>
-         *
-         * {@preformat java
-         *     interpolateOnBottomSide(polylinesOnTop[x].attach(polylineOnLeft));
-         * }
-         *
-         * @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[] swap = coordinates;
-            coordinates = source.coordinates;
-            size        = source.size;
-            opposite    = source.opposite;
-            if (opposite != null) {
-                opposite.opposite = this;
-            }
-            source.clear();
-            source.coordinates = swap;
-            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(Math.multiplyExact(size, 2), 32));
-            }
-            coordinates[size++] = x;
-            coordinates[size++] = y;
-        }
-
-        /**
-         * Returns a string representation of this {@code Polyline} for debugging purposes.
-         */
-        @Override
-        public String toString() {
-            return PathBuilder.toString(coordinates, size);
-        }
-    }
-
-    /**
-     * List of {@code Polyline} coordinates that have not yet been closed. Each {@code double[]} in this list is
-     * a copy of a {@code Polyline} used by {@link Level}. Those copies are performed for saving data before they
-     * are overwritten by next iterated cell.
-     *
-     * <h2>List indices and ordering of points</h2>
-     * For a given {@code Unclosed} list, all {@code double[]} arrays at even indices shall have their points read
-     * in reverse order and all {@code double[]} arrays at odd indices shall have their points read in forward order.
-     * The list size must be even and the list may contain null elements when there is no data in the corresponding
-     * iteration order. This convention makes easy to reverse the order of all points, simply by reversing the order
-     * of {@code double[]} arrays: because even indices become odd and odd indices become even, points order are
-     * implicitly reverted without the need to rewrite all {@code double[]} array contents.
-     *
-     * @see Level#partialPaths
-     */
-    @SuppressWarnings({"CloneableImplementsClone", "serial"})           // Not intended to be cloned or serialized.
-    private static final class Unclosed extends ArrayList<double[]> {
-        /**
-         * The first points and last point in this list of polylines. By convention the coordinate having fraction
-         * digits has all its bits inverted by the {@code ~} operator. May be {@code null} if a coordinate is NaN.
-         * Do not modify {@link Point} field values, because those instances are keys in {@link Level#partialPaths}.
-         */
-        private Point firstPoint, lastPoint;
-
-        /**
-         * Creates a list of polylines initialized to the given items.
-         * The given polylines and their opposite directions are cleared by this method.
-         *
-         * @param  polylineOnLeft  first polyline with points in forward order. Shall not be null.
-         * @param  polylineOnTop    next polyline with points in reverse order, or {@code null} if none.
-         */
-        Unclosed(final Polyline polylineOnLeft, final Polyline polylineOnTop) {
-            /*
-             * Search for first point and last point by inspecting `Polyline`s in the order shown below.
-             * The first 4 rows and the last 4 rows search for first point and last point respectively.
-             * The empty rows in the middle are an intentional gap for creating a regular pattern that
-             * we can exploit for 3 decisions that need to be done during the loop:
-             *
-             *     ✓ (index & 2) = 0    if using `polylineOnLeft` (otherwise `polylineOnTop`).
-             *     ✓ (index % 3) = 0    if using `opposite` value of polyline (may be null).
-             *     ✓ (index & 1) = 0    if fetching last point (otherwise fetch first point).
-             *
-             *  Index   Polyline   (iteration order)  !(i & 2)  !(i % 3)  !(i & 1)   Comment
-             *  ────────────────────────────────────────────────────────────────────────────
-             *   [0]    polylineOnLeft.opposite  (←)      ✓         ✓         ✓        (1)
-             *   [1]    polylineOnLeft           (→)      ✓                            (2)
-             *   [2]    polylineOnTop            (←)                          ✓        (1)
-             *   [3]    polylineOnTop.opposite   (→)                ✓                  (2)
-             *   [4]                                      ✓                   ✓
-             *   |5]                                      ✓
-             *   [6]    polylineOnTop.opposite   (→)                ✓         ✓        (3)
-             *   [7]    polylineOnTop            (←)                                   (4)
-             *   [8]    polylineOnLeft           (→)      ✓                   ✓        (3)
-             *   [9]    polylineOnLeft.opposite  (←)      ✓         ✓                  (4)
-             *
-             * Comments:
-             *   (1) Last  `Polyline` point is first `Unclosed` point because of reverse iteration order.
-             *   (2) First `Polyline` point is first `Unclosed` point because of forward iteration order.
-             *   (3) Last  `Polyline` point is last  `Unclosed` point because of forward iteration order.
-             *   (4) First `Polyline` point is last  `Unclosed` point because of reverse iteration order.
-             */
-            int index = 0;
-            do {
-                Polyline polyline = ((index & 2) == 0) ? polylineOnLeft : polylineOnTop;  // See above table (column 4).
-                if (index % 3 == 0 && polyline != null) polyline = polyline.opposite;     // See above table (column 5).
-                if (polyline != null) {
-                    int n = polyline.size;
-                    if (n != 0) {
-                        final double[] coordinates = polyline.coordinates;
-                        final double x, y;
-                        if (((index & 1) == 0)) {                          // See above table in comment (column 6).
-                            y = coordinates[--n];
-                            x = coordinates[--n];
-                        } else {
-                            x = coordinates[0];
-                            y = coordinates[1];
-                        }
-                        final boolean isLastPoint = (index >= 6);          // See row [6] in above table.
-                        if (Double.isFinite(x) && Double.isFinite(y)) {
-                            final Point p = new Point((int) x, (int) y);
-                            if (!Numerics.isInteger(x)) p.x = ~p.x;
-                            if (!Numerics.isInteger(y)) p.y = ~p.y;
-                            if (isLastPoint) {
-                                lastPoint = p;
-                                break;                                     // Done searching both points.
-                            }
-                            firstPoint = p;
-                        } else if (isLastPoint) {
-                            /*
-                             * If the last point was NaN, check if it was also the case of first point.
-                             * If yes, we will not be able to store this `Unclosed` in `partialPaths`
-                             * because we have no point that we can use as key (it would be pointless
-                             * to search for another point further in the `coordinates` array because
-                             * that point could never be matched with another `Unclosed`). Leave this
-                             * list empty for avoiding the copies done by `take(…)` calls. Instead,
-                             * callers should write polylines in `Level.path` immediately.
-                             */
-                            if (firstPoint == null) return;
-                            break;
-                        }
-                        /*
-                         * Done searching the first point (may still be null if that point is NaN).
-                         * Row [6] in above table is the first row for the search of last point.
-                         */
-                        index = 6;
-                        continue;
-                    }
-                }
-                if (++index == 4) {
-                    // Found no non-empty polylines during search for first point. No need to continue searching.
-                    return;
-                }
-            } while (index <= 9);
-            /*
-             * Copies coordinates only if at least one of `firstPoint` or `lastPoint` is a valid point.
-             */
-            take(polylineOnLeft.opposite);          // Point will be iterated in reverse order.
-            take(polylineOnLeft);                   // Point will be iterated in forward order.
-            if (polylineOnTop != null) {
-                Polyline suffix = polylineOnTop.opposite;
-                take(polylineOnTop);                // Inverse order. Set `polylineOnTop.opposite` to null.
-                take(suffix);                       // Forward order.
-            }
-        }
-
-        /**
-         * Takes a copy of coordinate values of given polyline, then clears that polyline.
-         */
-        private void take(final Polyline polyline) {
-            if (polyline != null && polyline.size != 0) {
-                add(Arrays.copyOf(polyline.coordinates, polyline.size));
-                polyline.clear();
-            } else {
-                add(null);                  // No data for iteration order at this position.
-            }
-        }
-
-        /**
-         * Returns {@code true} if the given point is equal to the start point or end point.
-         * This is used in assertions for checking key validity in {@link Level#partialPaths}.
-         */
-        final boolean isExtremity(final Point key) {
-            return key.equals(firstPoint) || key.equals(lastPoint);
-        }
-
-        /**
-         * Associates this polyline to its two extremities in the given map. If other polylines already exist
-         * for one or both extremities, then this polyline will be merged with previously existing polylines.
-         * This method returns {@code true} if the polyline has been closed, in which case caller should store
-         * the coordinates in {@link Level#path} immediately.
-         *
-         * @param  partialPaths  where to add or merge polylines.
-         * @return {@code true} if this polyline became a closed polygon as a result of merge operation.
-         */
-        final boolean addOrMerge(final Map<Point,Unclosed> partialPaths) {
-            final Unclosed before = partialPaths.remove(firstPoint);
-            final Unclosed after  = partialPaths.remove(lastPoint);
-            if (before != null) partialPaths.remove(addAll(before, true));
-            if (after  != null) partialPaths.remove(addAll(after, false));
-            if (firstPoint != null && firstPoint.equals(lastPoint)) {       // First/last points may have changed.
-                partialPaths.remove(firstPoint);
-                partialPaths.remove(lastPoint);
-                return true;
-            } else {
-                // Intentionally replace previous values.
-                if (firstPoint != null) partialPaths.put(firstPoint, this);
-                if (lastPoint  != null) partialPaths.put(lastPoint,  this);
-                return false;
-            }
-        }
-
-        /**
-         * Prepends or appends the given polylines to this list of polylines.
-         * Points order will be changed as needed in order to match extremities.
-         * The {@code other} instance should be forgotten after this method call.
-         *
-         * @param  other    the other polyline to append or prepend to this polyline.
-         * @param  prepend  {@code true} for prepend operation, {@code false} for append.
-         * @return extremity of {@code other} which has not been assigned to {@code this}.
-         */
-        private Point addAll(final Unclosed other, final boolean prepend) {
-            assert ((size() | other.size()) & 1) == 0;      // Must have even number of elements in both lists.
-            /*
-             * In figures below, ● are the extremities to attach together.
-             * `r` is a bitmask telling which polylines to reverse:
-             * 1=this, 2=other, together with combinations 0=none and 3=other.
-             */
-            int r; if ( lastPoint != null &&  lastPoint.equals(other.firstPoint)) r = 0;    // ○──────● ●──────○
-            else   if (firstPoint != null && firstPoint.equals(other.firstPoint)) r = 1;    // ●──────○ ●──────○
-            else   if ( lastPoint != null &&  lastPoint.equals(other. lastPoint)) r = 2;    // ○──────● ○──────●
-            else   if (firstPoint != null && firstPoint.equals(other. lastPoint)) r = 3;    // ●──────○ ○──────●
-            else {
-                // Should never happen because `other` has been obtained using a point of `this`.
-                throw new AssertionError();
-            }
-            if (prepend) r ^= 3;                      // Swap order in above  ○──○ ○──○  figures.
-            if ((r & 1) != 0)  this.reverse();
-            if ((r & 2) != 0) other.reverse();
-            if (prepend) {
-                addAll(0, other);
-                firstPoint = other.firstPoint;
-                return other.lastPoint;
-            } else {
-                addAll(other);
-                lastPoint = other.lastPoint;
-                return other.firstPoint;
-            }
-        }
-
-        /**
-         * Reverse the order of all points. The last polyline will become the first polyline and vice-versa.
-         * For each polyline, points will be iterated in opposite order. The trick on point order is done by
-         * moving polylines at even indices to odd indices, and conversely (see class javadoc for convention
-         * about even/odd indices).
-         */
-        private void reverse() {
-            Collections.reverse(this);
-            final Point swap = firstPoint;
-            firstPoint = lastPoint;
-            lastPoint = swap;
-        }
-
-        /**
-         * Returns the content of this list as an array of {@link Polyline} instances.
-         * {@code Polyline} instances at even index should be written with their points in reverse order.
-         *
-         * @see #writeTo(Joiner, Polyline[], boolean)
-         */
-        final Polyline[] toPolylines() {
-            final Polyline[] polylines = new Polyline[size()];
-            for (int i=0; i<polylines.length; i++) {
-                final double[] coordinates = get(i);
-                if (coordinates != null) {
-                    polylines[i] = new Polyline(coordinates);
-                }
-            }
-            return polylines;
-        }
-    }
-
-    /**
-     * Assembles arbitrary amount of {@link Polyline}s in a single Java2D {@link Shape} for a specific isoline level.
-     * This class extends {@link PathBuilder} with two additional features: remove spikes caused by ambiguities, then
-     * apply a {@link MathTransform} on all coordinate values.
-     *
-     * <h2>Spikes</h2>
-     * If the shape delimited by given polylines has a part with zero width or height ({@literal i.e.} a spike),
-     * truncates the polylines for removing that spike. This situation happens when some pixel values are exactly
-     * equal to isoline value, as in the picture below:
-     *
-     * {@preformat text
-     *     ●╌╌╌╲╌╌○╌╌╌╌╌╌○╌╌╌╌╌╌○╌╌╌╌╌╌○
-     *     ╎    ╲ ╎      ╎      ╎      ╎
-     *     ╎     ╲╎      ╎   →  ╎      ╎
-     *     ●╌╌╌╌╌╌●──────●──────●⤸╌╌╌╌╌○
-     *     ╎     ╱╎      ╎   ←  ╎      ╎
-     *     ╎    ╱ ╎      ╎      ╎      ╎
-     *     ●╌╌╌╱╌╌○╌╌╌╌╌╌○╌╌╌╌╌╌○╌╌╌╌╌╌○
-     * }
-     *
-     * The spike may appear or not depending on the convention adopted for strictly equal values.
-     * In above picture, the spike appears because the convention used in this implementation is:
-     *
-     * <ul>
-     *   <li>○: {@literal pixel value < isoline value}.</li>
-     *   <li>●: {@literal pixel value ≥ isoline value}.</li>
-     * </ul>
-     *
-     * If the following convention was used instead, the spike would not appear in above figure
-     * (but would appear in different situations):
-     *
-     * <ul>
-     *   <li>○: {@literal pixel value ≤ isoline value}.</li>
-     *   <li>●: {@literal pixel value > isoline value}.</li>
-     * </ul>
-     *
-     * This class detects and removes those spikes for avoiding convention-dependent results.
-     * We assume that spikes can appear only at the junction between two {@link Polyline} instances.
-     * Rational: having a spike require that we move forward then backward on the same coordinates,
-     * which is possible only with a non-null {@link Polyline#opposite} field.
-     */
-    private static final class Joiner extends PathBuilder {
-        /**
-         * Final transform to apply on coordinates, or {@code null} if none.
-         */
-        private final MathTransform gridToCRS;
-
-        /**
-         * Creates an initially empty set of isoline shapes.
-         */
-        Joiner(final MathTransform gridToCRS) {
-            this.gridToCRS = gridToCRS;
-        }
-
-        /**
-         * Detects and removes spikes for avoiding convention-dependent results.
-         * See {@link Joiner} class-javadoc for a description of the problem.
-         *
-         * <p>We perform the analysis in this method instead of in {@link #filterFull(double[], int)} on the
-         * the assumption that spikes can appear only between two calls to {@code append(…)} (because having a
-         * spike require that we move forward then backward on the same coordinates, which happen only with two
-         * distinct {@link Polyline} instances). It reduce the amount of coordinates to examine since we can check
-         * only the extremities instead of looking for spikes anywhere in the array.</p>
-         *
-         * @param  coordinates  the coordinates to filter. Values can be modified in-place.
-         * @param  lower        index of first coordinate to filter. Always even.
-         * @param  upper        index after the last coordinate to filter. Always even.
-         * @return number of valid coordinates after filtering.
-         */
-        @Override
-        protected int filterChunk(final double[] coordinates, final int lower, final int upper) {
-            int spike0 = lower;         // Will be index where (x,y) become different than (xo,yo).
-            int spike1 = lower;         // Idem, but searching forward instead of backward.
-            if (spike1 < upper) {
-                final double xo = coordinates[spike1++];
-                final double yo = coordinates[spike1++];
-                int equalityMask = 3;                   // Bit 1 set if (x == xo), bit 2 set if (y == yo).
-                while (spike1 < upper) {
-                    final int before = equalityMask;
-                    if (coordinates[spike1++] != xo) equalityMask &= ~1;
-                    if (coordinates[spike1++] != yo) equalityMask &= ~2;
-                    if (equalityMask == 0) {
-                        equalityMask = before;              // For keeping same search criterion.
-                        spike1 -= Polyline.DIMENSION;       // Restore previous position before mismatch.
-                        break;
-                    }
-                }
-                while (spike0 > 0) {
-                    if (coordinates[--spike0] != yo) equalityMask &= ~2;
-                    if (coordinates[--spike0] != xo) equalityMask &= ~1;
-                    if (equalityMask == 0) {
-                        spike0 += Polyline.DIMENSION;
-                        break;
-                    }
-                }
-            }
-            /*
-             * Here we have range of indices where the polygon has a width or height of zero.
-             * Search for a common point, then truncate at that point. Indices are like below:
-             *
-             *     0       spike0    lower          spike1         upper
-             *     ●────●────●────●────●────●────●────●────●────●────●
-             *                    └╌╌╌╌remove╌╌╌╌┘
-             * where:
-             *  - `lower` and `spike0` are inclusive.
-             *  - `upper` and `spike1` are exclusive.
-             *  - the region to remove are sowewhere between `spike0` and `spike1`.
-             */
-            final int limit = spike1;
-            int base;
-            while ((base = spike0 + 2*Polyline.DIMENSION) < limit) {    // Spikes exist only with at least 3 points.
-                final double xo = coordinates[spike0++];
-                final double yo = coordinates[spike0++];
-                spike1 = limit;
-                do {
-                    if (coordinates[spike1 - 2] == xo && coordinates[spike1 - 1] == yo) {
-                        /*
-                         * Remove points between the common point (xo,yo). The common point is kept on the
-                         * left side (`spike0` is already after that point) and removed on the right side.
-                         */
-                        System.arraycopy(coordinates, spike1, coordinates, spike0, upper - spike1);
-                        return upper - (spike1 - spike0);
-                    }
-                } while ((spike1 -= Polyline.DIMENSION) > base);
-            }
-            return upper;       // Nothing to remove.
-        }
-
-        /**
-         * Applies user-specified coordinate transform on all points of the whole polyline.
-         * This method is invoked after {@link #filterChunk(double[], int, int)}.
-         */
-        @Override
-        protected int filterFull(final double[] coordinates, final int upper) throws TransformException {
-            if (gridToCRS != null) {
-                gridToCRS.transform(coordinates, 0, coordinates, 0, upper / Polyline.DIMENSION);
-            }
-            return upper;
-        }
-    }
-
-    /**
-     * Writes all given polylines to the specified path builder. Null {@code Polyline} instances are ignored.
-     * {@code Polyline} instances at even index are written with their points in reverse order.
-     * All given polylines are cleared by this method.
-     *
-     * @param  path       where to write the polylines, or {@code null} if not yet created.
-     * @param  polylines  the polylines to write.
-     * @param  close      whether to close the polygon.
-     * @return the given path builder, or a newly created builder if the argument was null.
-     * @throws TransformException if the {@link #gridToCRS} transform can not be applied.
-     */
-    private Joiner writeTo(Joiner path, final Polyline[] polylines, final boolean close) throws TransformException {
-        for (int pi=0; pi < polylines.length; pi++) {
-            final Polyline p = polylines[pi];
-            if (p == null) {
-                continue;
-            }
-            final int size = p.size;
-            if (size == 0) {
-                assert p.isEmpty();
-                continue;
-            }
-            if (path == null) {
-                path = new Joiner(gridToCRS);
-            }
-            path.append(p.coordinates, size, (pi & 1) == 0);
-            p.clear();
-        }
-        if (path != null) {
-            path.createPolyline(close);
-        }
-        return path;
-    }
-}
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/TiledProcess.java b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/TiledProcess.java
index 8adf66e..79e3313 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/TiledProcess.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/TiledProcess.java
@@ -56,14 +56,14 @@
  * Consequently tile size will be determined by other considerations such as the number of processors.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.1
+ * @version 1.3
  *
  * @param  <R>  the type of value computed as a result of this process.
  *
  * @since 1.1
  * @module
  */
-abstract class TiledProcess<R> {
+public abstract class TiledProcess<R> {
     /**
      * Minimal "tile" size for sub-tasks computation. That size should not be too small because the
      * fork/join processes have some extra cost compared to processing the whole image as one single tile.
@@ -221,7 +221,7 @@
      * <p>This class implements {@link Callable} for {@code TiledProcess} convenience.
      * This implementation details should be ignored; it may change in any future version.</p>
      */
-    abstract class Task implements Callable<R> {
+    protected abstract class Task implements Callable<R> {
         /**
          * Index of this {@code Task} instance in the {@link TiledProcess#tasks} array.
          */
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
index a087331..6dd6be1 100644
--- 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
@@ -25,7 +25,7 @@
  *
  * @author  Johann Sorel (Geomatys)
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.1
+ * @version 1.3
  * @since   1.1
  * @module
  */
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Fragments.java b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Fragments.java
new file mode 100644
index 0000000..2392046
--- /dev/null
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Fragments.java
@@ -0,0 +1,274 @@
+/*
+ * 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.isoline;
+
+import java.awt.Point;
+import java.util.Map;
+import java.util.Arrays;
+import java.util.ArrayList;
+import java.util.Collections;
+import org.apache.sis.internal.util.Numerics;
+
+
+/**
+ * List of {@code PolylineBuffer} coordinates that have not yet been closed.
+ * Each {@code double[]} in this list is a copy of a {@link PolylineBuffer} used by {@link Tracer.Level}.
+ * Those copies are performed for saving data before they are overwritten by next iterated cell.
+ *
+ * <h2>List indices and ordering of points</h2>
+ * For a given {@code Fragments} list, all {@code double[]} arrays at even indices shall have their points read
+ * in reverse order and all {@code double[]} arrays at odd indices shall have their points read in forward order.
+ * The list size must be even and the list may contain null elements when there is no data in the corresponding
+ * iteration order. This convention makes easy to reverse the order of all points, simply by reversing the order
+ * of {@code double[]} arrays: because even indices become odd and odd indices become even, points order are
+ * implicitly reverted without the need to rewrite all {@code double[]} array contents.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.3
+ *
+ * @see Tracer.Level#partialPaths
+ *
+ * @since 1.1
+ * @module
+ */
+@SuppressWarnings({"CloneableImplementsClone", "serial"})           // Not intended to be cloned or serialized.
+final class Fragments extends ArrayList<double[]> {
+    /**
+     * The first points and last point in this list of polylines. By convention the coordinate having fraction
+     * digits has all its bits inverted by the {@code ~} operator. May be {@code null} if a coordinate is NaN.
+     * Do not modify {@link Point} field values, because those instances are keys in {@link Tracer.Level#partialPaths}.
+     */
+    private Point firstPoint, lastPoint;
+
+    /**
+     * Creates a list of polylines initialized to the given items.
+     * The given polylines and their opposite directions are cleared by this method.
+     *
+     * @param  polylineOnLeft  first polyline with points in forward order. Shall not be null.
+     * @param  polylineOnTop    next polyline with points in reverse order, or {@code null} if none.
+     */
+    Fragments(final PolylineBuffer polylineOnLeft, final PolylineBuffer polylineOnTop) {
+        /*
+         * Search for first and last point by inspecting `PolylineBuffer` instances in the order shown below.
+         * The first 4 rows and the last 4 rows search for first point and last point respectively.
+         * The empty rows in the middle are an intentional gap for creating a regular pattern that
+         * we can exploit for 3 decisions that need to be done during the loop:
+         *
+         *     ✓ (index & 2) = 0    if using `polylineOnLeft` (otherwise `polylineOnTop`).
+         *     ✓ (index % 3) = 0    if using `opposite` value of polyline (may be null).
+         *     ✓ (index & 1) = 0    if fetching last point (otherwise fetch first point).
+         *
+         *  Index   PolylineBuffer        (order) !(i & 2)  !(i % 3)  !(i & 1)   Comment
+         *  ────────────────────────────────────────────────────────────────────────────
+         *   [0]    polylineOnLeft.opposite  (←)      ✓         ✓         ✓        (1)
+         *   [1]    polylineOnLeft           (→)      ✓                            (2)
+         *   [2]    polylineOnTop            (←)                          ✓        (1)
+         *   [3]    polylineOnTop.opposite   (→)                ✓                  (2)
+         *   [4]                                      ✓                   ✓
+         *   |5]                                      ✓
+         *   [6]    polylineOnTop.opposite   (→)                ✓         ✓        (3)
+         *   [7]    polylineOnTop            (←)                                   (4)
+         *   [8]    polylineOnLeft           (→)      ✓                   ✓        (3)
+         *   [9]    polylineOnLeft.opposite  (←)      ✓         ✓                  (4)
+         *
+         * Comments:
+         *   (1) Last  `PolylineBuffer` point is first `Fragments` point because of reverse iteration order.
+         *   (2) First `PolylineBuffer` point is first `Fragments` point because of forward iteration order.
+         *   (3) Last  `PolylineBuffer` point is last  `Fragments` point because of forward iteration order.
+         *   (4) First `PolylineBuffer` point is last  `Fragments` point because of reverse iteration order.
+         */
+        int index = 0;
+        do {
+            PolylineBuffer polyline = ((index & 2) == 0) ? polylineOnLeft : polylineOnTop;  // See above table (column 4).
+            if (index % 3 == 0 && polyline != null) polyline = polyline.opposite;           // See above table (column 5).
+            if (polyline != null) {
+                int n = polyline.size;
+                if (n != 0) {
+                    final double[] coordinates = polyline.coordinates;
+                    final double x, y;
+                    if (((index & 1) == 0)) {                          // See above table in comment (column 6).
+                        y = coordinates[--n];
+                        x = coordinates[--n];
+                    } else {
+                        x = coordinates[0];
+                        y = coordinates[1];
+                    }
+                    final boolean isLastPoint = (index >= 6);          // See row [6] in above table.
+                    if (Double.isFinite(x) && Double.isFinite(y)) {
+                        final Point p = new Point((int) x, (int) y);
+                        if (!Numerics.isInteger(x)) p.x = ~p.x;
+                        if (!Numerics.isInteger(y)) p.y = ~p.y;
+                        if (isLastPoint) {
+                            lastPoint = p;
+                            break;                                     // Done searching both points.
+                        }
+                        firstPoint = p;
+                    } else if (isLastPoint) {
+                        /*
+                         * If the last point was NaN, check if it was also the case of first point.
+                         * If yes, we will not be able to store this `Fragments` in `partialPaths`
+                         * because we have no point that we can use as key (it would be pointless
+                         * to search for another point further in the `coordinates` array because
+                         * that point could never be matched with another `Fragments`). Leave this
+                         * list empty for avoiding the copies done by `take(…)` calls. Instead,
+                         * callers should write polylines in `Tracer.Level.path` immediately.
+                         */
+                        if (firstPoint == null) return;
+                        break;
+                    }
+                    /*
+                     * Done searching the first point (may still be null if that point is NaN).
+                     * Row [6] in above table is the first row for the search of last point.
+                     */
+                    index = 6;
+                    continue;
+                }
+            }
+            if (++index == 4) {
+                // Found no non-empty polylines during search for first point. No need to continue searching.
+                return;
+            }
+        } while (index <= 9);
+        /*
+         * Copies coordinates only if at least one of `firstPoint` or `lastPoint` is a valid point.
+         */
+        take(polylineOnLeft.opposite);          // Point will be iterated in reverse order.
+        take(polylineOnLeft);                   // Point will be iterated in forward order.
+        if (polylineOnTop != null) {
+            PolylineBuffer suffix = polylineOnTop.opposite;
+            take(polylineOnTop);                // Inverse order. Set `polylineOnTop.opposite` to null.
+            take(suffix);                       // Forward order.
+        }
+    }
+
+    /**
+     * Takes a copy of coordinate values of given polyline, then clears that polyline.
+     */
+    private void take(final PolylineBuffer polyline) {
+        if (polyline != null && polyline.size != 0) {
+            add(Arrays.copyOf(polyline.coordinates, polyline.size));
+            polyline.clear();
+        } else {
+            add(null);                  // No data for iteration order at this position.
+        }
+    }
+
+    /**
+     * Returns {@code true} if the given point is equal to the start point or end point.
+     * This is used in assertions for checking key validity in {@link Tracer.Level#partialPaths}.
+     */
+    final boolean isExtremity(final Point key) {
+        return key.equals(firstPoint) || key.equals(lastPoint);
+    }
+
+    /**
+     * Associates this polyline to its two extremities in the given map. If other polylines already exist
+     * for one or both extremities, then this polyline will be merged with previously existing polylines.
+     * This method returns {@code true} if caller should store the coordinates in {@link Tracer.Level#path}
+     * immediately. It may be either because the polyline has been closed as a polygon, or because the two
+     * extremities contains NaN values in which case the polylines are not anymore in {@code partialPaths}.
+     *
+     * @param  partialPaths  where to add or merge polylines.
+     * @return {@code true} if this polyline became a closed polygon as a result of merge operation.
+     */
+    final boolean addOrMerge(final Map<Point,Fragments> partialPaths) {
+        final Fragments before = partialPaths.remove(firstPoint);
+        final Fragments after  = partialPaths.remove(lastPoint);
+        if (before != null) partialPaths.remove(addAll(before, true));
+        if (after  != null) partialPaths.remove(addAll(after, false));
+        if (firstPoint != null && firstPoint.equals(lastPoint)) {       // First/last points may have changed.
+            partialPaths.remove(firstPoint);
+            partialPaths.remove(lastPoint);
+            return true;
+        } else {
+            // Intentionally replace previous values.
+            if (firstPoint != null) partialPaths.put(firstPoint, this);
+            if (lastPoint  != null) partialPaths.put(lastPoint,  this);
+            return firstPoint == null && lastPoint == null;
+        }
+    }
+
+    /**
+     * Prepends or appends the given polylines to this list of polylines.
+     * Points order will be changed as needed in order to match extremities.
+     * The {@code other} instance should be forgotten after this method call.
+     *
+     * @param  other    the other polyline to append or prepend to this polyline.
+     * @param  prepend  {@code true} for prepend operation, {@code false} for append.
+     * @return extremity of {@code other} which has not been assigned to {@code this}.
+     */
+    private Point addAll(final Fragments other, final boolean prepend) {
+        assert ((size() | other.size()) & 1) == 0;      // Must have even number of elements in both lists.
+        /*
+         * In figures below, ● are the extremities to attach together.
+         * `r` is a bitmask telling which polylines to reverse:
+         * 1=this, 2=other, together with combinations 0=none and 3=other.
+         */
+        int r; if ( lastPoint != null &&  lastPoint.equals(other.firstPoint)) r = 0;    // ○──────● ●──────○
+        else   if (firstPoint != null && firstPoint.equals(other.firstPoint)) r = 1;    // ●──────○ ●──────○
+        else   if ( lastPoint != null &&  lastPoint.equals(other. lastPoint)) r = 2;    // ○──────● ○──────●
+        else   if (firstPoint != null && firstPoint.equals(other. lastPoint)) r = 3;    // ●──────○ ○──────●
+        else {
+            // Should never happen because `other` has been obtained using a point of `this`.
+            throw new AssertionError();
+        }
+        if (prepend) r ^= 3;                      // Swap order in above  ○──○ ○──○  figures.
+        if ((r & 1) != 0)  this.reverse();
+        if ((r & 2) != 0) other.reverse();
+        if (prepend) {
+            addAll(0, other);
+            firstPoint = other.firstPoint;
+            return other.lastPoint;
+        } else {
+            addAll(other);
+            lastPoint = other.lastPoint;
+            return other.firstPoint;
+        }
+    }
+
+    /**
+     * Reverse the order of all points. The last polyline will become the first polyline and vice-versa.
+     * For each polyline, points will be iterated in opposite order. The trick on point order is done by
+     * moving polylines at even indices to odd indices, and conversely (see class javadoc for convention
+     * about even/odd indices).
+     */
+    private void reverse() {
+        Collections.reverse(this);
+        final Point swap = firstPoint;
+        firstPoint = lastPoint;
+        lastPoint = swap;
+    }
+
+    /**
+     * Returns the content of this list as an array of {@link PolylineBuffer} instances.
+     * {@code PolylineBuffer} instances at even index should be written with their points in reverse order.
+     *
+     * @return  elements of this array as polylines. May contain null elements.
+     *
+     * @see #writeTo(Joiner, PolylineBuffer[], boolean)
+     */
+    final PolylineBuffer[] toPolylines() {
+        final PolylineBuffer[] polylines = new PolylineBuffer[size()];
+        for (int i=0; i<polylines.length; i++) {
+            final double[] coordinates = get(i);
+            if (coordinates != null) {
+                polylines[i] = new PolylineBuffer(coordinates);
+            }
+        }
+        return polylines;
+    }
+}
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/isoline/Isolines.java
similarity index 67%
rename from core/sis-feature/src/main/java/org/apache/sis/internal/processing/image/Isolines.java
rename to core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Isolines.java
index abc87cd..76c642c 100644
--- 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/isoline/Isolines.java
@@ -14,16 +14,16 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package org.apache.sis.internal.processing.image;
+package org.apache.sis.internal.processing.isoline;
 
-import java.util.AbstractList;
 import java.util.Arrays;
 import java.util.List;
+import java.util.Map;
+import java.util.EnumMap;
 import java.util.TreeMap;
 import java.util.NavigableMap;
+import java.util.function.BiConsumer;
 import java.util.concurrent.Future;
-import java.util.concurrent.ExecutionException;
-import java.util.concurrent.CompletionException;
 import java.awt.Shape;
 import java.awt.geom.Path2D;
 import java.awt.image.RenderedImage;
@@ -33,11 +33,11 @@
 import org.apache.sis.image.PixelIterator;
 import org.apache.sis.util.ArgumentChecks;
 import org.apache.sis.util.ArraysExt;
-import org.apache.sis.util.resources.Errors;
+import org.apache.sis.util.Debug;
 
-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;
+import static org.apache.sis.internal.processing.isoline.Tracer.UPPER_LEFT;
+import static org.apache.sis.internal.processing.isoline.Tracer.UPPER_RIGHT;
+import static org.apache.sis.internal.processing.isoline.Tracer.LOWER_RIGHT;
 
 
 /**
@@ -46,7 +46,7 @@
  *
  * @author  Johann Sorel (Geomatys)
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.1
+ * @version 1.3
  *
  * @see <a href="https://en.wikipedia.org/wiki/Marching_squares">Marching squares on Wikipedia</a>
  *
@@ -55,16 +55,25 @@
  */
 public final class Isolines {
     /**
-     * Isoline data for each level, sorted in ascending order of {@link IsolineTracer.Level#value}.
+     * Isoline data for each level, sorted in ascending order of {@link Tracer.Level#value}.
      */
-    private final IsolineTracer.Level[] levels;
+    private final Tracer.Level[] levels;
+
+    /**
+     * A consumer to notify about the current state of isoline generation, or {@code null} if none.
+     * This is used for debugging purposes only. This field is temporarily set to a non-null value
+     * when using {@code StepsViewer} (in test package) for following an isoline generation step
+     * by step.
+     */
+    @Debug
+    private static final BiConsumer<String,Isolines> LISTENER = null;
 
     /**
      * Creates an initially empty set of isolines for the given levels. The given {@code values}
      * array should be one of the arrays validated by {@link #cloneAndSort(double[][])}.
      */
-    private Isolines(final IsolineTracer tracer, final int band, final double[] values, final int width) {
-        levels = new IsolineTracer.Level[values.length];
+    private Isolines(final Tracer tracer, final int band, final double[] values, final int width) {
+        levels = new Tracer.Level[values.length];
         for (int i=0; i<values.length; i++) {
             levels[i] = tracer.new Level(band, values[i], width);
         }
@@ -94,7 +103,7 @@
     }
 
     /**
-     * Sets the specified bit on {@link IsolineTracer.Level#isDataAbove} for all levels lower than given value.
+     * Sets the specified bit on {@link Tracer.Level#isDataAbove} for all levels lower than given value.
      *
      * <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
@@ -109,13 +118,13 @@
      * out in the final stage, when copying coordinates in {@link Path2D} objects.
      *
      * @param  value a sample values from the image.
-     * @param  bit   {@value IsolineTracer#UPPER_LEFT}, {@value IsolineTracer#UPPER_RIGHT},
-     *               {@value IsolineTracer#LOWER_LEFT} or {@value IsolineTracer#LOWER_RIGHT}.
+     * @param  bit   {@value Tracer#UPPER_LEFT}, {@value Tracer#UPPER_RIGHT},
+     *               {@value Tracer#LOWER_LEFT} or {@value Tracer#LOWER_RIGHT}.
      *
-     * @see IsolineTracer.Level#nextColumn()
+     * @see Tracer.Level#nextColumn()
      */
     private void setMaskBit(final double value, final int bit) {
-        for (final IsolineTracer.Level level : levels) {
+        for (final Tracer.Level level : levels) {
             if (level.value > value) break;                 // See above javadoc for NaN handling.
             level.isDataAbove |= bit;
         }
@@ -161,14 +170,35 @@
     {
         ArgumentChecks.ensureNonNull("data",   data);
         ArgumentChecks.ensureNonNull("levels", levels);
-        return new Process(data, cloneAndSort(levels), gridToCRS).execute();
+        return new Parallelized(data, cloneAndSort(levels), gridToCRS).execute();
+    }
+
+    /**
+     * Merges results of two partial computations (tiles).
+     * This is invoked only in multi-threaded isoline computation.
+     * The {@code isolines} instances are modified in-place.
+     * The {@code neighbor} instances can be discarded after the merge.
+     *
+     * @param  isolines  the first set of isolines to merge.
+     * @param  neighbor  another set of isolines close to the first set.
+     *
+     * @see Parallelized
+     */
+    static void merge(final Isolines[] isolines, final Isolines[] neighbor) throws TransformException {
+        for (int i=0; i<isolines.length; i++) {
+            final Isolines target = isolines[i];
+            final Isolines source = neighbor[i];
+            for (int j=0; j<target.levels.length; j++) {
+                target.levels[j].merge(source.levels[j]);
+            }
+        }
     }
 
     /**
      * Returns a provider of {@link PixelIterator} suitable to isoline computations.
      * It is critical that iterators use {@link SequenceType#LINEAR} iterator order.
      */
-    private static PixelIterator.Builder iterators() {
+    static PixelIterator.Builder iterators() {
         return new PixelIterator.Builder().setIteratorOrder(SequenceType.LINEAR);
     }
 
@@ -183,9 +213,9 @@
      * @return the {@code isolines} array, returned for convenience.
      * @throws TransformException if an error occurred during polylines creation.
      */
-    private static Isolines[] flush(final Isolines[] isolines) throws TransformException {
+    static Isolines[] flush(final Isolines[] isolines) throws TransformException {
         for (final Isolines isoline : isolines) {
-            for (final IsolineTracer.Level level : isoline.levels) {
+            for (final Tracer.Level level : isoline.levels) {
                 level.flush();
             }
         }
@@ -193,98 +223,21 @@
     }
 
     /**
-     * Wraps {@code Isolines.generate(…)} calculation in a process for parallel execution.
-     * The source image is divided in sub-region and the isolines in each sub-region will
-     * be computed in a different thread.
-     */
-    private static final class Process extends TiledProcess<Isolines[]> {
-        /**
-         * 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.
-         *
-         * @see #cloneAndSort(double[][])
-         */
-        private final double[][] levels;
-
-        /**
-         * Transform from image upper left corner (in pixel coordinates) to geometry coordinates.
-         */
-        private final MathTransform gridToCRS;
-
-        /**
-         * Creates a process for parallel isoline computation.
-         *
-         * @param  data       image providing source values.
-         * @param  levels     values for which to compute isolines.
-         * @param  gridToCRS  transform from pixel coordinates to geometry coordinates, or {@code null} if none.
-         */
-        Process(final RenderedImage data, final double[][] levels, final MathTransform gridToCRS) {
-            super(data, 1, 1, iterators());
-            this.levels = levels;
-            this.gridToCRS = gridToCRS;
-        }
-
-        /**
-         * Invoked by {@link TiledProcess} for creating a sub-task
-         * doing isoline computation on a sub-region of the image.
-         */
-        @Override
-        protected Task createSubTask() {
-            return new Tile();
-        }
-
-        /**
-         * A sub-task doing isoline computation on a sub-region of the image.
-         * The region is determined by the {@link #iterator}.
-         */
-        private final class Tile extends Task {
-            /** Isolines computed in the sub-region of this sub-task. */
-            private Isolines[] isolines;
-
-            /** Creates a new sub-task. */
-            Tile() {
-            }
-
-            /** Invoked in a background thread for performing isoline computation. */
-            @Override protected void execute() throws TransformException {
-                isolines = generate(iterator, levels, gridToCRS);
-            }
-
-            /** Invoked in a background thread for merging results of two sub-tasks. */
-            @Override protected void merge(final Task neighbor) throws TransformException {
-                for (int i=0; i<isolines.length; i++) {
-                    final Isolines target = isolines[i];
-                    final Isolines source = ((Tile) neighbor).isolines[i];
-                    for (int j=0; j<target.levels.length; j++) {
-                        target.levels[j].merge(source.levels[j]);
-                    }
-                }
-            }
-
-            /** Invoked on the last sub-task (after all merges) for getting final result. */
-            @Override protected Isolines[] result() throws TransformException {
-                return flush(isolines);
-            }
-        }
-    }
-
-    /**
      * Generates isolines for the specified levels in the region traversed by the given iterator.
      * This is the implementation of {@link #generate(RenderedImage, double[][], MathTransform)}
      * but potentially on a sub-region for parallel "fork-join" execution.
      */
-    private static Isolines[] generate(final PixelIterator iterator, final double[][] levels,
-                                       final MathTransform gridToCRS) throws TransformException
+    static Isolines[] generate(final PixelIterator iterator, final double[][] levels, final MathTransform gridToCRS)
+            throws TransformException
     {
         /*
          * Prepares a window of size 2×2 pixels over pixel values. Window elements are traversed
          * by incrementing indices in following order: band, column, row. The window content will
-         * be written in this method and read by IsolineTracer.
+         * be written in this method and read by Tracer.
          */
         final int numBands = iterator.getNumBands();
         final double[] window = new double[numBands * 4];
-        final IsolineTracer tracer = new IsolineTracer(window, numBands, iterator.getDomain(), gridToCRS);
+        final Tracer tracer = new Tracer(window, numBands, iterator.getDomain(), gridToCRS);
         /*
          * Prepare the set of isolines for each band in the image.
          * The number of cells on the horizontal axis is one less
@@ -325,7 +278,7 @@
              *
              *  - Get values on the 4 corners.
              *  - Save value of lower-left corner for use by next row.
-             *  - Initialize `IsolineTracer.Level.isDataAbove` bits for all levels.
+             *  - Initialize `Tracer.Level.isDataAbove` bits for all levels.
              *  - Interpolate the first cell.
              */
             System.arraycopy(valuesOnPreviousRow, 0, window, 0, twoPixels);
@@ -339,7 +292,7 @@
                 }
             }
             for (final Isolines iso : isolines) {
-                for (final IsolineTracer.Level level : iso.levels) {
+                for (final Tracer.Level level : iso.levels) {
                     level.interpolate();
                 }
             }
@@ -366,24 +319,29 @@
                 }
                 for (int b=0; b<numBands; b++) {
                     final Isolines iso = isolines[b];
-                    for (final IsolineTracer.Level level : iso.levels) {
+                    for (final Tracer.Level level : iso.levels) {
                         level.nextColumn();
                     }
                     iso.setMaskBit(window[numBands  + b], UPPER_RIGHT);
                     iso.setMaskBit(window[lastPixel + b], LOWER_RIGHT);
-                    for (final IsolineTracer.Level level : iso.levels) {
+                    for (final Tracer.Level level : iso.levels) {
                         level.interpolate();
                     }
                 }
             }
             /*
-             * Finished iteration on a row. Clear flags and update position
-             * before to move to next row.
+             * Finished iteration on a row. Clear flags and update position before to move to next row.
+             * If there is listeners to notify (for debugging purposes), notify them.
              */
             for (int b=0; b<numBands; b++) {
-                for (final IsolineTracer.Level level : isolines[b].levels) {
+                for (final Tracer.Level level : isolines[b].levels) {
                     level.finishedRow();
                 }
+                if (LISTENER != null) {
+                    final int y = tracer.y;
+                    final int h = iterator.getDomain().height;
+                    LISTENER.accept(String.format("After row %d of %d (%3.1f%%)", y, h, 100f*y/h), isolines[b]);
+                }
             }
             tracer.x = 0;
             tracer.y++;
@@ -392,9 +350,12 @@
          * Finished iteration over the whole image.
          */
         for (int b=0; b<numBands; b++) {
-            for (final IsolineTracer.Level level : isolines[b].levels) {
+            for (final Tracer.Level level : isolines[b].levels) {
                 level.finish();
             }
+            if (LISTENER != null) {
+                LISTENER.accept("Finished band " + b, isolines[b]);
+            }
         }
         return isolines;
     }
@@ -406,7 +367,7 @@
      */
     public final NavigableMap<Double,Shape> polylines() {
         final TreeMap<Double,Shape> paths = new TreeMap<>();
-        for (final IsolineTracer.Level level : levels) {
+        for (final Tracer.Level level : levels) {
             final Shape path = level.shape;
             if (path != null) {
                 paths.put(level.value, path);
@@ -421,7 +382,7 @@
      * @param  isolines  result of {@code generate(…)} or {@code parallelGenerate(…)} method call.
      * @return isoline shapes for each values in each band.
      */
-    private static NavigableMap<Double,Shape>[] toArray(final Isolines[] isolines) {
+    static NavigableMap<Double,Shape>[] toArray(final Isolines[] isolines) {
         @SuppressWarnings({"rawtypes", "unchecked"})
         final NavigableMap<Double,Shape>[] result = new NavigableMap[isolines.length];
         for (int i=0; i<result.length; i++) {
@@ -452,55 +413,18 @@
     }
 
     /**
-     * Deferred isoline result, created when computation is continuing in background.
-     * The {@link Future} result is requested the first time that {@link #get(int)} is invoked.
+     * Returns the pixel coordinates of all level, for debugging purposes only.
+     * The {@link #gridToCRS} transform is <em>not</em> applied by this method.
+     * For avoiding confusing behavior, that transform should be null.
+     *
+     * @return the pixel coordinates.
      */
-    private static final class Result extends AbstractList<NavigableMap<Double,Shape>> {
-        /** The task computing isolines result. Reset to {@code null} when no longer needed. */
-        private Future<Isolines[]> task;
-
-        /** The result of {@link Future#get()} fetched when first needed. */
-        private NavigableMap<Double,Shape>[] isolines;
-
-        /** Creates a new list for the given future isolines. */
-        Result(final Future<Isolines[]> task) {
-            this.task = task;
+    @Debug
+    final Map<PolylineStage,Path2D> toRawPath() {
+        final Map<PolylineStage,Path2D> appendTo = new EnumMap<>(PolylineStage.class);
+        for (final Tracer.Level level : levels) {
+            level.toRawPath(appendTo);
         }
-
-        /** Fetches the isolines from the {@link Future} if not already done. */
-        @SuppressWarnings("ReturnOfCollectionOrArrayField")
-        private NavigableMap<Double,Shape>[] isolines() {
-            if (isolines == null) {
-                if (task == null) {
-                    throw new CompletionException(Errors.format(Errors.Keys.BackgroundComputationFailed), null);
-                }
-                try {
-                    isolines = Isolines.toArray(task.get());
-                    task = null;
-                } catch (InterruptedException e) {
-                    // Do not clear `task`: the result may become available later.
-                    throw new CompletionException(Errors.format(Errors.Keys.InterruptedWhileWaitingResult), e);
-                } catch (ExecutionException e) {
-                    task = null;
-                    throw new CompletionException(Errors.format(Errors.Keys.BackgroundComputationFailed), e.getCause());
-                }
-            }
-            return isolines;
-        }
-
-        /** Returns the list length, which is the number of bands. */
-        @Override public int size() {
-            return isolines().length;
-        }
-
-        /** Returns the isolines in the given band. */
-        @Override public NavigableMap<Double,Shape> get(final int band) {
-            return isolines()[band];
-        }
-
-        /** Returns the list content as an array. */
-        @Override public Object[] toArray() {
-            return isolines().clone();
-        }
+        return appendTo;
     }
 }
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Joiner.java b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Joiner.java
new file mode 100644
index 0000000..b9ccfcf
--- /dev/null
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Joiner.java
@@ -0,0 +1,168 @@
+/*
+ * 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.isoline;
+
+import org.apache.sis.internal.feature.j2d.PathBuilder;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.TransformException;
+
+
+/**
+ * Assembles arbitrary amount of {@link PolylineBuffer}s in a single Java2D {@link Shape} for an isoline level.
+ * This class extends {@link PathBuilder} with two additional features: remove spikes caused by ambiguities,
+ * then apply a {@link MathTransform} on all coordinate values.
+ *
+ * <h2>Spikes</h2>
+ * If the shape delimited by given polylines has a part with zero width or height ({@literal i.e.} a spike),
+ * truncates the polylines for removing that spike. This situation happens when some pixel values are exactly
+ * equal to isoline value, as in the picture below:
+ *
+ * {@preformat text
+ *     ●╌╌╌╲╌╌○╌╌╌╌╌╌○╌╌╌╌╌╌○╌╌╌╌╌╌○
+ *     ╎    ╲ ╎      ╎      ╎      ╎
+ *     ╎     ╲╎      ╎   →  ╎      ╎
+ *     ●╌╌╌╌╌╌●──────●──────●⤸╌╌╌╌╌○
+ *     ╎     ╱╎      ╎   ←  ╎      ╎
+ *     ╎    ╱ ╎      ╎      ╎      ╎
+ *     ●╌╌╌╱╌╌○╌╌╌╌╌╌○╌╌╌╌╌╌○╌╌╌╌╌╌○
+ * }
+ *
+ * The spike may appear or not depending on the convention adopted for strictly equal values.
+ * In above picture, the spike appears because the convention used in this implementation is:
+ *
+ * <ul>
+ *   <li>○: {@literal pixel value < isoline value}.</li>
+ *   <li>●: {@literal pixel value ≥ isoline value}.</li>
+ * </ul>
+ *
+ * If the following convention was used instead, the spike would not appear in above figure
+ * (but would appear in different situations):
+ *
+ * <ul>
+ *   <li>○: {@literal pixel value ≤ isoline value}.</li>
+ *   <li>●: {@literal pixel value > isoline value}.</li>
+ * </ul>
+ *
+ * This class detects and removes those spikes for avoiding convention-dependent results.
+ * We assume that spikes can appear only at the junction between two {@link PolylineBuffer} instances.
+ * Rational: having a spike require that we move forward then backward on the same coordinates,
+ * which is possible only with a non-null {@link PolylineBuffer#opposite} field.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.3
+ * @since   1.1
+ * @module
+ */
+final class Joiner extends PathBuilder {
+    /**
+     * Final transform to apply on coordinates, or {@code null} if none.
+     */
+    private final MathTransform gridToCRS;
+
+    /**
+     * Creates an initially empty set of isoline shapes.
+     */
+    Joiner(final MathTransform gridToCRS) {
+        this.gridToCRS = gridToCRS;
+    }
+
+    /**
+     * Detects and removes spikes for avoiding convention-dependent results.
+     * See {@link Joiner} class-javadoc for a description of the problem.
+     *
+     * <p>We perform the analysis in this method instead of in {@link #filterFull(double[], int)} on the
+     * the assumption that spikes can appear only between two calls to {@code append(…)} (because having
+     * a spike requires that we move forward then backward on the same coordinates, which happen only with
+     * two distinct {@link PolylineBuffer} instances). It reduce the amount of coordinates to examine since
+     * we can check only the extremities instead of looking for spikes anywhere in the array.</p>
+     *
+     * @param  coordinates  the coordinates to filter. Values can be modified in-place.
+     * @param  lower        index of first coordinate to filter. Always even.
+     * @param  upper        index after the last coordinate to filter. Always even.
+     * @return number of valid coordinates after filtering.
+     */
+    @Override
+    protected int filterChunk(final double[] coordinates, final int lower, final int upper) {
+        int spike0 = lower;         // Will be index where (x,y) become different than (xo,yo).
+        int spike1 = lower;         // Idem, but searching forward instead of backward.
+        if (spike1 < upper) {
+            final double xo = coordinates[spike1++];
+            final double yo = coordinates[spike1++];
+            int equalityMask = 3;                   // Bit 1 set if (x == xo), bit 2 set if (y == yo).
+            while (spike1 < upper) {
+                final int before = equalityMask;
+                if (coordinates[spike1++] != xo) equalityMask &= ~1;
+                if (coordinates[spike1++] != yo) equalityMask &= ~2;
+                if (equalityMask == 0) {
+                    equalityMask = before;                  // For keeping same search criterion.
+                    spike1 -= PolylineBuffer.DIMENSION;      // Restore previous position before mismatch.
+                    break;
+                }
+            }
+            while (spike0 > 0) {
+                if (coordinates[--spike0] != yo) equalityMask &= ~2;
+                if (coordinates[--spike0] != xo) equalityMask &= ~1;
+                if (equalityMask == 0) {
+                    spike0 += PolylineBuffer.DIMENSION;
+                    break;
+                }
+            }
+        }
+        /*
+         * Here we have range of indices where the polygon has a width or height of zero.
+         * Search for a common point, then truncate at that point. Indices are like below:
+         *
+         *     0       spike0    lower          spike1         upper
+         *     ●────●────●────●────●────●────●────●────●────●────●
+         *                    └╌╌╌╌remove╌╌╌╌┘
+         * where:
+         *  - `lower` and `spike0` are inclusive.
+         *  - `upper` and `spike1` are exclusive.
+         *  - the region to remove are sowewhere between `spike0` and `spike1`.
+         */
+        final int limit = spike1;
+        int base;
+        while ((base = spike0 + 2*PolylineBuffer.DIMENSION) < limit) {    // Spikes exist only with at least 3 points.
+            final double xo = coordinates[spike0++];
+            final double yo = coordinates[spike0++];
+            spike1 = limit;
+            do {
+                if (coordinates[spike1 - 2] == xo && coordinates[spike1 - 1] == yo) {
+                    /*
+                     * Remove points between the common point (xo,yo). The common point is kept on the
+                     * left side (`spike0` is already after that point) and removed on the right side.
+                     */
+                    System.arraycopy(coordinates, spike1, coordinates, spike0, upper - spike1);
+                    return upper - (spike1 - spike0);
+                }
+            } while ((spike1 -= PolylineBuffer.DIMENSION) > base);
+        }
+        return upper;       // Nothing to remove.
+    }
+
+    /**
+     * Applies user-specified coordinate transform on all points of the whole polyline.
+     * This method is invoked after {@link #filterChunk(double[], int, int)}.
+     */
+    @Override
+    protected int filterFull(final double[] coordinates, final int upper) throws TransformException {
+        if (gridToCRS != null) {
+            gridToCRS.transform(coordinates, 0, coordinates, 0, upper / PolylineBuffer.DIMENSION);
+        }
+        return upper;
+    }
+}
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Parallelized.java b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Parallelized.java
new file mode 100644
index 0000000..eab2514
--- /dev/null
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Parallelized.java
@@ -0,0 +1,112 @@
+/*
+ * 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.isoline;
+
+import java.awt.image.RenderedImage;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.internal.processing.image.TiledProcess;
+
+
+/**
+ * Wraps {@code Isolines.generate(…)} calculation in a process for parallel execution.
+ * The source image is divided in sub-region and the isolines in each sub-region will
+ * be computed in a different thread.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.3
+ * @since   1.1
+ * @module
+ */
+final class Parallelized extends TiledProcess<Isolines[]> {
+    /**
+     * 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.
+     *
+     * @see #cloneAndSort(double[][])
+     */
+    private final double[][] levels;
+
+    /**
+     * Transform from image upper left corner (in pixel coordinates) to geometry coordinates.
+     */
+    private final MathTransform gridToCRS;
+
+    /**
+     * Creates a process for parallel isoline computation.
+     *
+     * @param  data       image providing source values.
+     * @param  levels     values for which to compute isolines.
+     * @param  gridToCRS  transform from pixel coordinates to geometry coordinates, or {@code null} if none.
+     */
+    Parallelized(final RenderedImage data, final double[][] levels, final MathTransform gridToCRS) {
+        super(data, 1, 1, Isolines.iterators());
+        this.levels = levels;
+        this.gridToCRS = gridToCRS;
+    }
+
+    /**
+     * Invoked by {@link TiledProcess} for creating a sub-task
+     * doing isoline computation on a sub-region of the image.
+     */
+    @Override
+    protected Task createSubTask() {
+        return new Tile();
+    }
+
+    /**
+     * A sub-task doing isoline computation on a sub-region of the image.
+     * The region is determined by the {@link #iterator}.
+     */
+    private final class Tile extends Task {
+        /**
+         * Isolines computed in the sub-region of this sub-task.
+         */
+        private Isolines[] isolines;
+
+        /**
+         * Creates a new sub-task.
+         */
+        Tile() {
+        }
+
+        /**
+         * Invoked in a background thread for performing isoline computation.
+         */
+        @Override
+        protected void execute() throws TransformException {
+            isolines = Isolines.generate(iterator, levels, gridToCRS);
+        }
+
+        /**
+         * Invoked in a background thread for merging results of two sub-tasks.
+         */
+        @Override
+        protected void merge(final Task neighbor) throws TransformException {
+            Isolines.merge(isolines, ((Tile) neighbor).isolines);
+        }
+
+        /**
+         * Invoked on the last sub-task (after all merges) for getting final result.
+         */
+        @Override
+        protected Isolines[] result() throws TransformException {
+            return Isolines.flush(isolines);
+        }
+    }
+}
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/PolylineBuffer.java b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/PolylineBuffer.java
new file mode 100644
index 0000000..d7c1687
--- /dev/null
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/PolylineBuffer.java
@@ -0,0 +1,210 @@
+/*
+ * 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.isoline;
+
+import java.util.Map;
+import java.util.Arrays;
+import java.awt.geom.Path2D;
+import org.apache.sis.internal.feature.j2d.PathBuilder;
+import org.apache.sis.util.ArraysExt;
+import org.apache.sis.util.Debug;
+
+
+/**
+ * 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 Tracer.Level#path}
+ * and the {@code PolylineBuffer} object is recycled for a new polyline.</p>
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.3
+ * @since   1.1
+ * @module
+ */
+final class PolylineBuffer {
+    /**
+     * Number of coordinates in a tuple.
+     */
+    static final int DIMENSION = 2;
+
+    /**
+     * Coordinates as (x,y) tuples. This array is expanded as needed.
+     */
+    double[] coordinates;
+
+    /**
+     * Number of valid elements in the {@link #coordinates} array.
+     * This is twice the number of points.
+     */
+    int size;
+
+    /**
+     * If the polyline has points added to its two extremities, the other extremity. Otherwise {@code null}.
+     * The first point of {@code opposite} polyline is connected to the first point of this polyline.
+     * Consequently when those two polylines are joined in a single polyline, the coordinates of either
+     * {@code this} or {@code opposite} must be iterated in reverse order.
+     */
+    PolylineBuffer opposite;
+
+    /**
+     * Creates an initially empty polyline.
+     */
+    PolylineBuffer() {
+        coordinates = ArraysExt.EMPTY_DOUBLE;
+    }
+
+    /**
+     * Creates a new polyline wrapping the given coordinates. Used for wrapping {@link Fragments}
+     * instances in objects expected by {@link Tracer#writeTo(Joiner, Polyline[], boolean)}.
+     * Those {@code Polyline} instances are temporary.
+     */
+    PolylineBuffer(final double[] data) {
+        coordinates = data;
+        size = data.length;
+    }
+
+    /**
+     * Discards all coordinates in this polyline. This method does not clear
+     * the {@link #opposite} polyline; it is caller's responsibility to do so.
+     */
+    final void clear() {
+        opposite = null;
+        size = 0;
+    }
+
+    /**
+     * Returns whether this polyline is empty. This method is used only for {@code assert isEmpty()}
+     * statement because of the check for {@code opposite == null}: an empty polyline should not have
+     * a non-null {@link #opposite} value.
+     */
+    final boolean isEmpty() {
+        return size == 0 & (opposite == null);
+    }
+
+    /**
+     * 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).
+     *
+     * <p>This method is typically invoked in the following pattern (but this is not mandatory).
+     * An important aspect is that {@code this} and {@code other} should be on perpendicular axes:</p>
+     *
+     * {@preformat java
+     *     interpolateOnBottomSide(polylinesOnTop[x].attach(polylineOnLeft));
+     * }
+     *
+     * @return {@code this} for method calls chaining.
+     */
+    final PolylineBuffer attach(final PolylineBuffer 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 PolylineBuffer transferFrom(final PolylineBuffer source) {
+        assert isEmpty();
+        final double[] swap = coordinates;
+        coordinates = source.coordinates;
+        size        = source.size;
+        opposite    = source.opposite;
+        if (opposite != null) {
+            opposite.opposite = this;
+        }
+        source.clear();
+        source.coordinates = swap;
+        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(Math.multiplyExact(size, 2), 32));
+        }
+        coordinates[size++] = x;
+        coordinates[size++] = y;
+    }
+
+    /**
+     * Returns a string representation of this {@code Polyline} for debugging purposes.
+     */
+    @Override
+    public String toString() {
+        return PathBuilder.toString(coordinates, size);
+    }
+
+    /**
+     * Appends the pixel coordinates of this polyline to the given path, for debugging purposes only.
+     * The {@link #gridToCRS} transform is <em>not</em> applied by this method.
+     * For avoiding confusing behavior, that transform should be null.
+     *
+     * @param  appendTo  where to append the coordinates.
+     *
+     * @see Tracer.Level#toRawPath(Map)
+     */
+    @Debug
+    final void toRawPath(final Map<PolylineStage,Path2D> appendTo) {
+        PolylineStage.BUFFER.add(appendTo, coordinates, size);
+    }
+}
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/PolylineStage.java b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/PolylineStage.java
new file mode 100644
index 0000000..477aa3a
--- /dev/null
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/PolylineStage.java
@@ -0,0 +1,116 @@
+/*
+ * 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.isoline;
+
+import java.util.Map;
+import java.awt.Shape;
+import java.awt.geom.Path2D;
+import org.apache.sis.util.Debug;
+
+
+/**
+ * Tells at which stage are the polylines represented by a Java2D {@link Shape}.
+ * A set of polylines way still be under construction in {@link PolylineBuffer}
+ * during iteration over pixel values, or the polylines may have been classified
+ * as incomplete after iteration over a row, or the polylines may be final result.
+ *
+ * <p>This is used only for debugging purposes because end users should see only the final result.
+ * This information allows {@code StepsViewer} (in test package) to use different colors for different stages.</p>
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.3
+ * @since   1.3
+ * @module
+ */
+@Debug
+enum PolylineStage {
+    /**
+     * The polylines are under construction in various {@link PolylineBuffer} instances.
+     * This is the first stage, which happens during iteration over pixel values.
+     */
+    BUFFER,
+
+    /**
+     * The polylines are no longer in the buffers filled by the iteration over pixel values,
+     * but are still incomplete. It happens when, after finishing iteration over a row, some
+     * polylines will not be continued by iteration on the next row and those polylines have
+     * not yet been closed as polygons. Those polyline fragments are moved to a "pending" list,
+     * as they may be closed later after more polylines fragments become available.
+     */
+    FRAGMENT,
+
+    /**
+     * The polylines are final result to be show to user.
+     */
+    FINAL;
+
+    /**
+     * Returns the destination where to write polylines for this stage.
+     *
+     * @param  appendTo  map of path for different stages.
+     * @return the path to use for writing polylines at this stage.
+     */
+    private Path2D destination(final Map<PolylineStage,Path2D> appendTo) {
+        return appendTo.computeIfAbsent(this, (k) -> new Path2D.Float());
+    }
+
+    /**
+     * Adds coordinates to the specified map.
+     *
+     * @param  appendTo     where to append the coordinates.
+     * @param  coordinates  (x,y) tuples to append, starting with the coordinate at index 0.
+     * @param  size         number of coordinates to add (twice the number of tuples).
+     */
+    final void add(final Map<PolylineStage,Path2D> appendTo, final double[] coordinates, final int size) {
+        int i = 0;
+        if (i < size) {
+            final Path2D p = destination(appendTo);
+            p.moveTo(coordinates[i++], coordinates[i++]);
+            while (i < size) {
+                p.lineTo(coordinates[i++], coordinates[i++]);
+            }
+        }
+    }
+
+    /**
+     * Adds polylines in the values of the given map. Keys are ignored.
+     *
+     * @param  appendTo      where to append the coordinates.
+     * @param  partialPaths  map of polylines to add.
+     */
+    final void add(final Map<PolylineStage,Path2D> appendTo, final Map<?,Fragments> partialPaths) {
+        for (final Fragments fragment : partialPaths.values()) {
+            for (final double[] coordinates : fragment) {
+                if (coordinates != null) {
+                    add(appendTo, coordinates, coordinates.length);
+                }
+            }
+        }
+    }
+
+    /**
+     * Adds polylines to the specified map.
+     *
+     * @param  appendTo   where to append the polylines.
+     * @param  polylines  the polylines to append to the map, or {@code null} if none.
+     */
+    final void add(final Map<PolylineStage,Path2D> appendTo, final Shape polylines) {
+        if (polylines != null) {
+            destination(appendTo).append(polylines, false);
+        }
+    }
+}
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Result.java b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Result.java
new file mode 100644
index 0000000..9f3e797
--- /dev/null
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Result.java
@@ -0,0 +1,101 @@
+/*
+ * 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.isoline;
+
+import java.awt.Shape;
+import java.util.AbstractList;
+import java.util.NavigableMap;
+import java.util.concurrent.Future;
+import java.util.concurrent.CompletionException;
+import java.util.concurrent.ExecutionException;
+import org.apache.sis.util.resources.Errors;
+
+
+/**
+ * Deferred isoline result, created when computation is continuing in background.
+ * The {@link Future} result is requested the first time that {@link #get(int)} is invoked.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.3
+ * @since   1.1
+ * @module
+ */
+final class Result extends AbstractList<NavigableMap<Double,Shape>> {
+    /**
+     * The task computing isolines result. Reset to {@code null} when no longer needed.
+     */
+    private Future<Isolines[]> task;
+
+    /**
+     * The result of {@link Future#get()} fetched when first needed.
+     */
+    private NavigableMap<Double,Shape>[] isolines;
+
+    /**
+     * Creates a new list for the given future isolines.
+     */
+    Result(final Future<Isolines[]> task) {
+        this.task = task;
+    }
+
+    /**
+     * Fetches the isolines from the {@link Future} if not already done.
+     */
+    @SuppressWarnings("ReturnOfCollectionOrArrayField")
+    private NavigableMap<Double,Shape>[] isolines() {
+        if (isolines == null) {
+            if (task == null) {
+                throw new CompletionException(Errors.format(Errors.Keys.BackgroundComputationFailed), null);
+            }
+            try {
+                isolines = Isolines.toArray(task.get());
+                task = null;
+            } catch (InterruptedException e) {
+                // Do not clear `task`: the result may become available later.
+                throw new CompletionException(Errors.format(Errors.Keys.InterruptedWhileWaitingResult), e);
+            } catch (ExecutionException e) {
+                task = null;
+                throw new CompletionException(Errors.format(Errors.Keys.BackgroundComputationFailed), e.getCause());
+            }
+        }
+        return isolines;
+    }
+
+    /**
+     * Returns the list length, which is the number of bands.
+     */
+    @Override
+    public int size() {
+        return isolines().length;
+    }
+
+    /**
+     * Returns the isolines in the given band.
+     */
+    @Override
+    public NavigableMap<Double,Shape> get(final int band) {
+        return isolines()[band];
+    }
+
+    /**
+     * Returns the list content as an array.
+     */
+    @Override
+    public Object[] toArray() {
+        return isolines().clone();
+    }
+}
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Tracer.java b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Tracer.java
new file mode 100644
index 0000000..e88944a
--- /dev/null
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/Tracer.java
@@ -0,0 +1,739 @@
+/*
+ * 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.isoline;
+
+import java.util.Map;
+import java.util.HashMap;
+import java.util.IdentityHashMap;
+import java.awt.Point;
+import java.awt.Rectangle;
+import java.awt.Shape;
+import java.awt.geom.Path2D;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.internal.feature.j2d.PathBuilder;
+import org.apache.sis.util.Debug;
+
+
+/**
+ * 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 Tracer} is used for doing a single iteration over all image pixels.
+ *
+ * @author  Johann Sorel (Geomatys)
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.3
+ *
+ * @see <a href="https://en.wikipedia.org/wiki/Marching_squares">Marching squares on Wikipedia</a>
+ *
+ * @since 1.1
+ * @module
+ */
+final class Tracer {
+    /**
+     * 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.
+     * The length of this array is <var>(number of bands)</var> × 2 (width) × 2 (height).
+     */
+    private final double[] 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.
+     * The range is 0 inclusive to {@code domain.width} exclusive.
+     */
+    int x;
+
+    /**
+     * Pixel coordinate on the top side of the cell where to interpolate.
+     * The range is 0 inclusive to {@code domain.height} exclusive.
+     */
+    int y;
+
+    /**
+     * Translation to apply on coordinates. For isolines computed sequentially, this is the image origin
+     * (often 0,0 but not necessarily). For isolines computed in parallel, the translations are different
+     * for each computation tile.
+     */
+    private final double translateX, translateY;
+
+    /**
+     * Final transform to apply on coordinates (integer source coordinates at pixel centers).
+     * Can be {@code null} if none.
+     */
+    private final MathTransform gridToCRS;
+
+    /**
+     * Creates a new position for the given data window.
+     *
+     * @param  window       the 2×2 window containing pixel values in the 4 corners of current contouring grid cell.
+     * @param  pixelStride  increment to the position in {@code window} for reading next sample value.
+     * @param  domain       pixel coordinates where iteration will happen.
+     * @param  gridToCRS    final transform to apply on coordinates (integer source coordinates at pixel centers).
+     */
+    Tracer(final double[] window, final int pixelStride, final Rectangle domain, final MathTransform gridToCRS) {
+        this.window      = window;
+        this.pixelStride = pixelStride;
+        this.translateX  = domain.x;
+        this.translateY  = domain.y;
+        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 {
+        /**
+         * Band number where to read values in the {@link #window} array.
+         */
+        private final int band;
+
+        /**
+         * The level value.
+         *
+         * @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 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 be the same. It could change the way line segments are assembled in a single {@link PolylineBuffer},
+         * 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 {@link PolylineBuffer#isEmpty() is empty} if the cell in previous iteration was like below
+         * (no line cross the right border):
+         *
+         * {@preformat text
+         *     ○╌╲╌╌╌╌●              ○╌╌╌┼╌╌●
+         *     ╎  ╲   ╎              ╎   │  ╎
+         *     ╎   ╲  ╎              ╎   │  ╎
+         *     ○╌╌╌╌╲╌●              ○╌╌╌┼╌╌●
+         * }
+         */
+        private final PolylineBuffer polylineOnLeft;
+
+        /**
+         * The polylines in each column which need to be continued on the next row.
+         * This array contains empty instances in columns where there are no polylines 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 PolylineBuffer[] polylinesOnTop;
+
+        /**
+         * Paths that have not yet been closed. The {@link PolylineBuffer} coordinates are copied in this map when
+         * iteration finished on a row but the polyline under construction will not be continued by the next row,
+         * or when the {@link #closeLeftWithTop(PolylineBuffer)} method has been invoked but the geometry to close
+         * is still not complete. This map accumulates those partial shapes for assembling them later when missing
+         * parts become available.
+         *
+         * <h4>Map keys</h4>
+         * Keys are grid coordinates rounded toward 0. The coordinate having fraction digits has its bits inverted
+         * by the {@code ~} operator. For each point, there is at most one coordinate having such fraction digits.
+         *
+         * <h4>Map values</h4>
+         * {@code Fragments} instances are list of {@code double[]} arrays to be concatenated in a single polygon later.
+         * For a given {@code Fragments} list, all {@code double[]} arrays at even indices shall have their points read
+         * in reverse order and all {@code double[]} arrays at odd indices shall have their points read in forward order.
+         * The list may contain null elements when there is no data in the corresponding iteration order.
+         *
+         * @see #closeLeftWithTop(PolylineBuffer)
+         */
+        private final Map<Point,Fragments> partialPaths;
+
+        /**
+         * Builder of isolines as a Java2D shape, created when first needed.
+         * The {@link PolylineBuffer} coordinates are copied in this path when a geometry is closed
+         * and transformed using {@link #gridToCRS}. This is almost final result; the only difference
+         * compared to {@link #shape} is that the coordinates are not yet wrapped in a {@link Shape}.
+         *
+         * @see #writeTo(Joiner, PolylineBuffer[], boolean)
+         * @see PolylineStage#FINAL
+         */
+        private Joiner path;
+
+        /**
+         * The isolines as a Java2D shape, created by {@link #finish()}.
+         * This is the shape to be returned to user for this level after we finished to process all cells.
+         *
+         * @see PolylineStage#FINAL
+         */
+        Shape shape;
+
+        /**
+         * Creates new isoline levels for the given value.
+         *
+         * @param  band   band number where to read values in the {@link #window} array.
+         * @param  value  the isoline level value.
+         * @param  width  the contouring grid cell width (one cell smaller than image width).
+         */
+        Level(final int band, final double value, final int width) {
+            this.band      = band;
+            this.value     = value;
+            partialPaths   = new HashMap<>();
+            polylineOnLeft = new PolylineBuffer();
+            polylinesOnTop = new PolylineBuffer[width];
+            for (int i=0; i<width; i++) {
+                polylinesOnTop[i] = new PolylineBuffer();
+            }
+        }
+
+        /**
+         * 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.
+         *
+         * @see Isolines#setMaskBit(double, int)
+         */
+        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
+         * the {@link PathBuilder}.
+         */
+        @SuppressWarnings("AssertWithSideEffects")
+        final void interpolate() throws TransformException {
+            /*
+             * Note: `interpolateMissingLeftSide()` and `interpolateMissingTopSide(…)` should do interpolations
+             * only for cells in the first column and first row respectively. We could avoid those method calls
+             * for all other cells if we add two flags in the `isDataAbove` bitmask: FIRST_ROW and FIRST_COLUMN.
+             * The switch cases then become something like below:
+             *
+             *     case <bitmask> | FIRST_COLUMN | FIRST_ROW:
+             *     case <bitmask> | FIRST_COLUMN: {
+             *         interpolateMissingLeftSide();
+             *         // Fall through
+             *     }
+             *     case <bitmask> | FIRST_ROW:
+             *     case <bitmask>: {
+             *         // Interpolations on other borders.
+             *         break;
+             *     }
+             *
+             * We tried that approach, but benchmarking on Java 15 suggested a small performance decrease
+             * instead of an improvement. It may be worth to try again in the future, after advancement
+             * in compiler technology.
+             */
+            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();
+                    interpolateMissingLeftSide();
+                    interpolateOnRightSide();                   // Will be the left side of next column.
+                    break;
+                }
+                /*     ○╌╌╌┼╌╌●        ●╌╌╌┼╌╌○
+                 *     ╎   │  ╎        ╎   │  ╎
+                 *     ╎   │  ╎        ╎   │  ╎
+                 *     ○╌╌╌┼╌╌●        ●╌╌╌┼╌╌○
+                 */
+                case UPPER_RIGHT | LOWER_RIGHT:
+                case UPPER_LEFT  | LOWER_LEFT: {
+                    assert polylineOnLeft.isEmpty();
+                    final PolylineBuffer polylineOnTop = polylinesOnTop[x];
+                    interpolateMissingTopSide(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();
+                    interpolateMissingLeftSide();
+                    interpolateOnBottomSide(polylinesOnTop[x].transferFrom(polylineOnLeft));
+                    break;
+                }
+                /*     ○╌╌╌╌╲╌●        ●╌╌╌╌╲╌○
+                 *     ╎     ╲╎        ╎     ╲╎
+                 *     ╎      ╲        ╎      ╲
+                 *     ○╌╌╌╌╌╌○╲       ●╌╌╌╌╌╌●╲
+                 */
+                case UPPER_RIGHT:
+                case UPPER_LEFT | LOWER_LEFT | LOWER_RIGHT: {
+                    assert polylineOnLeft.isEmpty();
+                    interpolateMissingTopSide(polylineOnLeft.transferFrom(polylinesOnTop[x]));
+                    interpolateOnRightSide();
+                    break;
+                }
+                /*     ○╌╌╌╌╌╌○╱       ●╌╌╌╌╌╌●╱
+                 *     ╎      ╱        ╎      ╱
+                 *     ╎     ╱╎        ╎     ╱╎
+                 *     ○╌╌╌╌╱╌●        ●╌╌╌╌╱╌○
+                 */
+                case LOWER_RIGHT:
+                case UPPER_LEFT | UPPER_RIGHT | LOWER_LEFT: {
+                    assert polylinesOnTop[x].isEmpty();
+                    assert polylineOnLeft   .isEmpty();
+                    interpolateOnRightSide();
+                    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: {
+                    closeLeftWithTop(polylinesOnTop[x]);
+                    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 double[] data = window;
+                        int p = band;
+                        do average += data[p];
+                        while ((p += pixelStride) < data.length);
+                        assert (p -= band) == pixelStride * 4 : p;
+                        average /= 4;
+                    }
+                    boolean LLtoUR = isDataAbove == (LOWER_LEFT | UPPER_RIGHT);
+                    LLtoUR ^= (average <= value);
+                    final PolylineBuffer polylineOnTop = polylinesOnTop[x];
+                    if (LLtoUR) {
+                        closeLeftWithTop(polylineOnTop);
+                        interpolateOnRightSide();
+                        interpolateOnBottomSide(polylineOnTop.attach(polylineOnLeft));
+                    } else {
+                        interpolateMissingLeftSide();
+                        final PolylineBuffer swap = new PolylineBuffer().transferFrom(polylineOnTop);
+                        interpolateOnBottomSide(polylineOnTop.transferFrom(polylineOnLeft));
+                        interpolateMissingTopSide(polylineOnLeft.transferFrom(swap));
+                        interpolateOnRightSide();
+                    }
+                    break;
+                }
+            }
+        }
+
+        /**
+         * Appends to {@link #polylineOnLeft} a point interpolated on the left side if that point is missing.
+         * This interpolation should happens only in the first column.
+         */
+        private void interpolateMissingLeftSide() {
+            if (polylineOnLeft.size == 0) {
+                polylineOnLeft.append(translateX + (x),
+                                      translateY + (y + interpolate(0, 2*pixelStride)));
+            }
+        }
+
+        /**
+         * Appends to {@code polylineOnTop} a point interpolated on the top side if that point is missing.
+         * This interpolation should happens only in the first row.
+         */
+        private void interpolateMissingTopSide(final PolylineBuffer polylineOnTop) {
+            if (polylineOnTop.size == 0) {
+                interpolateOnTopSide(polylineOnTop);
+            }
+        }
+
+        /**
+         * Appends to the given polyline a point interpolated on the top side.
+         */
+        private void interpolateOnTopSide(final PolylineBuffer appendTo) {
+            appendTo.append(translateX + (x + interpolate(0, pixelStride)),
+                            translateY + (y));
+        }
+
+        /**
+         * Appends to {@link #polylineOnLeft} a point interpolated on the right side.
+         * The polyline on right side will become {@code polylineOnLeft} in next column.
+         */
+        private void interpolateOnRightSide() {
+            polylineOnLeft.append(translateX + (x + 1),
+                                  translateY + (y + interpolate(pixelStride, 3*pixelStride)));
+        }
+
+        /**
+         * Appends to the given polyline a point interpolated on the bottom side.
+         * The polyline on top side will become a {@code polylineOnBottoù} in next row.
+         */
+        private void interpolateOnBottomSide(final PolylineBuffer polylineOnTop) {
+            polylineOnTop.append(translateX + (x + interpolate(2*pixelStride, 3*pixelStride)),
+                                 translateY + (y + 1));
+        }
+
+        /**
+         * Interpolates the position where the isoline passes between two values.
+         *
+         * @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 double[] data = window;
+            final int    p  = band;
+            final double v1 = data[p + i1];
+            final double v2 = data[p + i2];
+            return (value - v1) / (v2 - v1);
+        }
+
+        /**
+         * Joins {@link #polylineOnLeft} with {@code polylineOnTop}, saves their coordinates
+         * and clear those {@link PolylineBuffer} instances for use in next cell.
+         * The coordinates are written directly to {@link #path} if we got a closed polygon,
+         * or otherwise are saved in {@link #partialPaths} for later processing.
+         * This method is invoked for cells like below:
+         *
+         * {@preformat text
+         *     ●╌╱╌╌╌╌○        ○╌╱╌╌╌╌●        ○╌╱╌╌╌╌●╱
+         *     ╎╱     ╎        ╎╱     ╎        ╎╱     ╱
+         *     ╱      ╎        ╱      ╎        ╱     ╱╎
+         *    ╱○╌╌╌╌╌╌○       ╱●╌╌╌╌╌╌●       ╱●╌╌╌╌╱╌○
+         * }
+         *
+         * This method does itself the interpolations on left side and top side. The two polylines
+         * {@link #polylineOnLeft} and {@code polylineOnTop} will become empty after this method call.
+         *
+         * @param  polylineOnTop  value of {@code polylinesOnTop[x]}.
+         * @throws TransformException if the {@link Tracer#gridToCRS} transform can not be applied.
+         */
+        private void closeLeftWithTop(final PolylineBuffer polylineOnTop) throws TransformException {
+            interpolateMissingLeftSide();
+            interpolateMissingTopSide(polylineOnTop);
+            final PolylineBuffer[] polylines;
+            if (polylineOnLeft.opposite == polylineOnTop) {
+                assert polylineOnTop.opposite == polylineOnLeft;
+                /*
+                 * We have a loop: the polygon can be closed now, without copying coordinates to temporary buffers.
+                 * Points in `PolylineBuffer` instances will be iterated in (reverse, forward) order respectively.
+                 * Consequently the points we just interpolated will be first point and last point before closing.
+                 */
+                polylines = new PolylineBuffer[] {polylineOnTop, polylineOnLeft};    // (reverse, forward) point order.
+            } else {
+                /*
+                 * Joining left and top polylines do not yet create a closed shape. Consequently we may not write
+                 * in the `path` now. But maybe we can close the polygon later after more polylines are attached.
+                 */
+                final Fragments fragment = new Fragments(polylineOnLeft, polylineOnTop);
+                if (fragment.isEmpty()) {
+                    /*
+                     * Fragment starts and ends with NaN values. We will not be able to complete a polygon.
+                     * Better to write the polylines now for avoiding temporary copies of their coordinates.
+                     */
+                    polylines = new PolylineBuffer[] {
+                        polylineOnLeft.opposite, polylineOnLeft, polylineOnTop, polylineOnTop.opposite
+                    };
+                } else if (fragment.addOrMerge(partialPaths)) {
+                    /*
+                     * The fragment has been merged with previously existing fragments and became a polygon.
+                     * We can write the polygon immediately. There are no more references to those coordinates
+                     * in the `partialPaths` map.
+                     */
+                    polylines = fragment.toPolylines();
+                } else {
+                    return;
+                }
+            }
+            path = writeTo(path, polylines, true);
+        }
+
+        /**
+         * Writes the content of given polyline without closing it as a polygon.
+         * The given polyline will become empty after this method call.
+         */
+        private void writeFragment(final PolylineBuffer polyline) throws TransformException {
+            final Fragments fragment = new Fragments(polyline, null);
+            final PolylineBuffer[] polylines;
+            final boolean close;
+            if (fragment.isEmpty()) {
+                close = false;
+                polylines = new PolylineBuffer[] {polyline.opposite, polyline};     // (reverse, forward) point order.
+            } else {
+                close = fragment.addOrMerge(partialPaths);
+                if (!close) {
+                    // Keep in `partialPaths`. Maybe it can be closed later.
+                    return;
+                }
+                polylines = fragment.toPolylines();
+            }
+            path = writeTo(path, polylines, close);
+        }
+
+        /**
+         * Invoked after iteration on a single row has been completed. If there is a polyline
+         * finishing on the right image border, the coordinates needs to be saved somewhere
+         * because that {@link PolylineBuffer} will not be continued by cells on next rows.
+         */
+        final void finishedRow() throws TransformException {
+            if (!polylineOnLeft.transferToOpposite()) {
+                writeFragment(polylineOnLeft);
+            }
+            isDataAbove = 0;
+        }
+
+        /**
+         * Invoked after the iteration has been completed on the full area of interest.
+         * This method writes all remaining polylines to {@link #partialPaths}.
+         * It assumes that {@link #finishedRow()} has already been invoked.
+         * This {@link Level} instance can not be used anymore after this call.
+         */
+        final void finish() throws TransformException {
+            assert polylineOnLeft.isEmpty();
+            polylineOnLeft.coordinates = null;
+            /*
+             * This method sets various values to null for letting the garbage collector do its work.
+             * This is okay for a `Level` instance which is not going to be used anymore, except for
+             * reading the `shape` field.
+             */
+            for (int i=0; i < polylinesOnTop.length; i++) {
+                writeFragment(polylinesOnTop[i]);
+                polylinesOnTop[i] = null;
+            }
+            assert isConsistent();
+        }
+
+        /**
+         * Verifies that {@link #partialPaths} consistency. Used for assertions only.
+         */
+        private boolean isConsistent() {
+            for (final Map.Entry<Point,Fragments> entry : partialPaths.entrySet()) {
+                if (!entry.getValue().isExtremity(entry.getKey())) return false;
+            }
+            return true;
+        }
+
+        /**
+         * Transfers all {@code other} polylines into this instance. The {@code other} instance should be a neighbor,
+         * i.e. an instance sharing a border with this instance. The {@code other} instance will become empty after
+         * this method call.
+         *
+         * @param  other  a neighbor level (on top, left, right or bottom) to merge with this level.
+         * @throws TransformException if an error occurred during polylines creation.
+         */
+        final void merge(final Level other) throws TransformException {
+            assert other != this && other.value == value;
+            if (path == null) {
+                path = other.path;
+            } else {
+                path.append(other.path);
+            }
+            other.path = null;
+            assert  this.isConsistent();
+            assert other.isConsistent();
+            final IdentityHashMap<Fragments,Boolean> done = new IdentityHashMap<>(other.partialPaths.size() / 2);
+            for (final Map.Entry<Point,Fragments> entry : other.partialPaths.entrySet()) {
+                final Fragments fragment = entry.getValue();
+                if (done.put(fragment, Boolean.TRUE) == null) {
+                    assert fragment.isExtremity(entry.getKey());
+                    if (fragment.addOrMerge(partialPaths)) {
+                        path = writeTo(path, fragment.toPolylines(), true);
+                        fragment.clear();
+                    }
+                }
+                entry.setValue(null);       // Let the garbage collector do its work.
+            }
+        }
+
+        /**
+         * Flushes any pending {@link #partialPaths} to {@link #path}. This method is invoked after
+         * {@link #finish()} has been invoked for all sub-regions (many sub-regions may exist if
+         * isoline generation has been splitted for parallel computation).
+         *
+         * @throws TransformException if an error occurred during polylines creation.
+         */
+        final void flush() throws TransformException {
+            for (final Map.Entry<Point,Fragments> entry : partialPaths.entrySet()) {
+                final Fragments fragment = entry.getValue();
+                assert fragment.isExtremity(entry.getKey());
+                if (!fragment.isEmpty()) {
+                    path = writeTo(path, fragment.toPolylines(), false);
+                    fragment.clear();       // Necessary because the same list appears twice in the map.
+                }
+                entry.setValue(null);       // Let the garbage collector do its work.
+            }
+            if (path != null) {
+                shape = path.build();
+                path  = null;
+            }
+        }
+
+        /**
+         * Appends the pixel coordinates of this level to the given path, for debugging purposes only.
+         * The {@link #gridToCRS} transform is <em>not</em> applied by this method.
+         * For avoiding confusing behavior, that transform should be null.
+         *
+         * @param  appendTo  where to append the coordinates.
+         *
+         * @see Isolines#toRawPath()
+         */
+        @Debug
+        final void toRawPath(final Map<PolylineStage,Path2D> appendTo) {
+            PolylineStage.FINAL.add(appendTo, (path != null) ? path.snapshot() : shape);
+            PolylineStage.FRAGMENT.add(appendTo, partialPaths);
+            polylineOnLeft.toRawPath(appendTo);
+            for (final PolylineBuffer p : polylinesOnTop) {
+                if (p != null) p.toRawPath(appendTo);
+            }
+        }
+    }
+
+    /**
+     * Writes all given polylines to the specified path builder. Null {@code PolylineBuffer} instances are ignored.
+     * {@code PolylineBuffer} instances at even index are written with their points in reverse order.
+     * All given polylines are cleared by this method.
+     *
+     * @param  path       where to write the polylines, or {@code null} if not yet created.
+     * @param  polylines  the polylines to write.
+     * @param  close      whether to close the polygon.
+     * @return the given path builder, or a newly created builder if the argument was null.
+     * @throws TransformException if the {@link #gridToCRS} transform can not be applied.
+     */
+    private Joiner writeTo(Joiner path, final PolylineBuffer[] polylines, final boolean close) throws TransformException {
+        for (int pi=0; pi < polylines.length; pi++) {
+            final PolylineBuffer p = polylines[pi];
+            if (p == null) {
+                continue;
+            }
+            final int size = p.size;
+            if (size == 0) {
+                assert p.isEmpty();
+                continue;
+            }
+            if (path == null) {
+                path = new Joiner(gridToCRS);
+            }
+            path.append(p.coordinates, size, (pi & 1) == 0);
+            p.clear();
+        }
+        if (path != null) {
+            path.createPolyline(close);
+        }
+        return path;
+    }
+}
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/package-info.java b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/package-info.java
new file mode 100644
index 0000000..8f337f5
--- /dev/null
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/processing/isoline/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.
+ */
+
+/**
+ * Isolines generation.
+ *
+ * <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.3
+ * @since   1.1
+ * @module
+ */
+package org.apache.sis.internal.processing.isoline;
diff --git a/core/sis-feature/src/test/java/org/apache/sis/image/ImageProcessorTest.java b/core/sis-feature/src/test/java/org/apache/sis/image/ImageProcessorTest.java
index 22dd7ae..6a123ba 100644
--- a/core/sis-feature/src/test/java/org/apache/sis/image/ImageProcessorTest.java
+++ b/core/sis-feature/src/test/java/org/apache/sis/image/ImageProcessorTest.java
@@ -20,7 +20,7 @@
 import java.awt.Shape;
 import java.awt.image.BufferedImage;
 import java.awt.image.RenderedImage;
-import org.apache.sis.internal.processing.image.IsolinesTest;
+import org.apache.sis.internal.processing.isoline.IsolinesTest;
 import org.apache.sis.test.DependsOn;
 import org.opengis.referencing.operation.MathTransform;
 import org.apache.sis.test.TestCase;
@@ -38,7 +38,7 @@
  * @since   1.1
  * @module
  */
-@DependsOn(org.apache.sis.internal.processing.image.IsolinesTest.class)
+@DependsOn(org.apache.sis.internal.processing.isoline.IsolinesTest.class)
 public final strictfp class ImageProcessorTest extends TestCase {
     /**
      * Tests {@link ImageProcessor#isolines(RenderedImage, double[][], MathTransform)}.
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/isoline/IsolinesTest.java
similarity index 99%
rename from core/sis-feature/src/test/java/org/apache/sis/internal/processing/image/IsolinesTest.java
rename to core/sis-feature/src/test/java/org/apache/sis/internal/processing/isoline/IsolinesTest.java
index 14568a0..e5fa451 100644
--- 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/isoline/IsolinesTest.java
@@ -14,7 +14,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package org.apache.sis.internal.processing.image;
+package org.apache.sis.internal.processing.isoline;
 
 import java.util.Map;
 import java.awt.Shape;
diff --git a/core/sis-feature/src/test/java/org/apache/sis/internal/processing/isoline/StepsViewer.java b/core/sis-feature/src/test/java/org/apache/sis/internal/processing/isoline/StepsViewer.java
new file mode 100644
index 0000000..9003ee3
--- /dev/null
+++ b/core/sis-feature/src/test/java/org/apache/sis/internal/processing/isoline/StepsViewer.java
@@ -0,0 +1,426 @@
+/*
+ * 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.isoline;
+
+import java.util.Map;
+import java.util.EnumMap;
+import java.awt.Shape;
+import java.awt.Color;
+import java.awt.Graphics;
+import java.awt.Graphics2D;
+import java.awt.BasicStroke;
+import java.awt.EventQueue;
+import java.awt.BorderLayout;
+import java.awt.Container;
+import java.awt.Dimension;
+import java.awt.Rectangle;
+import java.awt.event.ActionEvent;
+import java.awt.event.ActionListener;
+import java.awt.geom.AffineTransform;
+import java.awt.geom.Path2D;
+import java.awt.geom.PathIterator;
+import java.awt.image.BufferedImage;
+import java.awt.image.RenderedImage;
+import java.awt.image.WritableRaster;
+import javax.swing.Timer;
+import javax.swing.JFrame;
+import javax.swing.JPanel;
+import javax.swing.JLabel;
+import javax.swing.JButton;
+import javax.swing.JComponent;
+import javax.swing.ButtonModel;
+import javax.swing.event.ChangeEvent;
+import javax.swing.event.ChangeListener;
+import java.util.function.BiConsumer;
+import java.util.concurrent.CountDownLatch;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.internal.referencing.j2d.AffineTransform2D;
+
+import static org.junit.Assert.*;
+
+
+/**
+ * A viewer for showing isoline generation step-by-step.
+ * For enabling the use of this class, temporarily remove {@code private} and {@code final} keywords in
+ * {@link Isolines#LISTENER}, then uncomment the {@link #setListener(StepsViewer)} constructor body.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.3
+ * @since   1.3
+ * @module
+ */
+@SuppressWarnings("serial")
+public final class StepsViewer extends JComponent implements BiConsumer<String,Isolines>, ChangeListener, ActionListener {
+    /**
+     * Sets the component to be notified after each row of isolines generated from the rendered image.
+     * The body of this method is commented-out because {@link Isolines#LISTENER} is private and final.
+     * The body should be uncommented only temporarily during debugging phases.
+     */
+    private static void setListener(final StepsViewer listener) {
+        // Isolines.LISTENER = listener;
+    }
+
+    /**
+     * Entry point for debugging. Edit this method body as needed for loading an image to use as test data.
+     *
+     * @param  args  ignored.
+     * @throws Exception if an error occurred during I/O or isoline generation.
+     */
+    public static void main(final String[] args) throws Exception {
+        // showStepByStep(local.test.DebugIsoline.data(), 0);
+    }
+
+    /**
+     * Size of the window and spacing between borders and isolines. All values are in pixels.
+     */
+    private static final int CANVAS_WIDTH = 1600, CANVAS_HEIGHT = 1000, PADDING = 3;
+
+    /**
+     * Whether to flip X and/or Y axis.
+     */
+    private static final boolean FLIP_X = false, FLIP_Y = true;
+
+    /**
+     * Description of current step. This title is updated at each isoline generation step,
+     * when {@link #accept(String, Shape)} is invoked.
+     */
+    private final JLabel stepTitle;
+
+    /**
+     * The button for moving to the next step. When this button is enabled, the isoline process is blocked
+     * by {@link #blocker} until this button is pressed. When this button is pressed, the isoline process
+     * continue until {@link #accept(String, Shape)} is invoked again.
+     *
+     * @see #actionPerformed(ActionEvent)
+     */
+    private final JButton next;
+
+    /**
+     * Simulate a "next" action after some delay. This is used when users keep the "Next" button pressed.
+     */
+    private final Timer delayedNext;
+
+    /**
+     * Blocks the isoline computation thread until the user is ready to see the next step.
+     */
+    private CountDownLatch blocker;
+
+    /**
+     * The isolines to show.
+     */
+    private final Map<PolylineStage,Path2D> isolines;
+
+    /**
+     * The colors to associate to the isoline for each stage.
+     * Array indices are {@link PolylineStage#ordinal()} values.
+     */
+    private final Color[] stageColors;
+
+    /**
+     * Bounds of {@link #isolines}, slightly expanded for making easier to see.
+     */
+    private Rectangle bounds;
+
+    /**
+     * Conversion from pixel indices in the source image to pixel indices in the displayed window.
+     */
+    private final AffineTransform2D sourceToCanvas;
+
+    /**
+     * The image to show in background. This image has the size of the canvas.
+     */
+    private final BufferedImage background;
+
+    /**
+     * The final result returned by public API.
+     */
+    private Shape result;
+
+    /**
+     * Creates a new viewer.
+     *
+     * @param  data  the source of data for isolines.
+     * @param  pane  the container where to add components.
+     */
+    @SuppressWarnings("ThisEscapedInObjectConstruction")
+    private StepsViewer(final RenderedImage data, final Container pane) {
+        isolines    = new EnumMap<>(PolylineStage.class);
+        stageColors = new Color[] {Color.YELLOW, Color.PINK, Color.BLUE};
+        /*
+         * Computes a transform from indices in the data matrix to pixel coordinates in the canvas.
+         */
+        setMaximumSize(new Dimension(CANVAS_WIDTH, CANVAS_HEIGHT));
+        final double scaleX = (CANVAS_WIDTH  - 2*PADDING) / (double) data.getWidth();
+        final double scaleY = (CANVAS_HEIGHT - 2*PADDING) / (double) data.getHeight();
+        sourceToCanvas = new AffineTransform2D(
+                FLIP_X ? -scaleX : scaleX, 0, 0, FLIP_Y ? -scaleY : scaleY,
+                scaleX * (PADDING + data.getMinX() + (FLIP_X ? data.getWidth()  : 0)),
+                scaleY * (PADDING + data.getMinY() + (FLIP_Y ? data.getHeight() : 0)));
+        /*
+         * Creates a background image as a grayscale image with pixel values in the range 0 to 64.
+         * It produces a dark image for making the isolines (in bright colors) easier to see.
+         */
+        background = new BufferedImage(CANVAS_WIDTH, CANVAS_HEIGHT, BufferedImage.TYPE_BYTE_GRAY);
+        {   // For keeping variable locale.
+            final Graphics2D gh = background.createGraphics();
+            gh.drawRenderedImage(data, sourceToCanvas);
+            gh.dispose();
+            final WritableRaster r = background.getRaster();
+            int min = 256, max = 2;
+            for (int y=0; y<CANVAS_HEIGHT; y++) {
+                for (int x=0; x<CANVAS_WIDTH; x++) {
+                    final int value = r.getSample(x, y, 0);
+                    if (value != 0) {
+                        if (value < min) min = value;
+                        if (value > max) max = value;
+                    }
+                }
+            }
+            if (--min < max) {
+                final float scale = 128f / (max - min);
+                for (int y=0; y<CANVAS_HEIGHT; y++) {
+                    for (int x=0; x<CANVAS_WIDTH; x++) {
+                        final int value = r.getSample(x, y, 0);
+                        r.setSample(x, y, 0, Math.max(Math.round(scale * (value - min)), 0));
+                    }
+                }
+            }
+        }
+        /*
+         * Swing controls.
+         */
+        stepTitle = new JLabel();
+        next = new JButton("Next");
+        next.setEnabled(false);
+        next.addActionListener(this);
+        next.getModel().addChangeListener(this);
+        delayedNext = new Timer(1000, this::fastForward);       // 1 second delay before fast forward.
+        delayedNext.setRepeats(false);
+
+        final JPanel bar = new JPanel(new BorderLayout());
+        bar .add(stepTitle, BorderLayout.CENTER);
+        bar .add(next,      BorderLayout.EAST);
+        pane.add(bar,       BorderLayout.NORTH);
+        pane.add(this,      BorderLayout.CENTER);
+    }
+
+    /**
+     * Generates isolines for the given image and show the result step by step.
+     * The given image shall have only one band.
+     *
+     * @param  data    the source of data for isolines.
+     * @param  levels  levels of isolones to generate.
+     */
+    public static void showStepByStep(final RenderedImage data, final double... levels) {
+        assertEquals("Unsupported number of bands.", 1, data.getSampleModel().getNumBands());
+        final JFrame frame = new JFrame("Step-by-step isoline viewer");
+        frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
+        frame.setLayout(new BorderLayout());
+        final StepsViewer viewer = new StepsViewer(data, frame.getContentPane());
+        final Isolines iso;
+        try {
+            setListener(viewer);
+            frame.setVisible(true);
+            frame.setSize(CANVAS_WIDTH, CANVAS_HEIGHT);
+            iso = Isolines.generate(data, new double[][] {levels}, null)[0];
+        } catch (TransformException e) {
+            throw new AssertionError(e);        // Should not happen because we specified an identity transform.
+        } finally {
+            setListener(null);
+        }
+        final Path2D path = new Path2D.Float();
+        for (final Shape shape : iso.polylines().values()) {
+            path.append(shape, false);
+        }
+        path.transform(viewer.sourceToCanvas);
+        viewer.result = path;
+        viewer.repaint();
+    }
+
+    /**
+     * Invoked when the isolines need to be drawn.
+     */
+    @Override
+    protected void paintComponent(final Graphics g) {
+        super.paintComponent(g);
+        final Graphics2D gh = (Graphics2D) g;
+        gh.drawRenderedImage(background, new AffineTransform());
+        if (bounds != null) {
+            gh.setStroke(new BasicStroke(2));
+            gh.setColor(Color.ORANGE);
+            gh.draw(bounds);
+        }
+        for (final Map.Entry<PolylineStage,Path2D> entry : isolines.entrySet()) {
+            final PolylineStage stage = entry.getKey();
+            gh.setStroke(new BasicStroke((result == null) ? (stage == PolylineStage.BUFFER ? 2 : 0) : 5));
+            gh.setColor((result == null) ? stageColors[stage.ordinal()] : Color.YELLOW);
+            gh.draw(entry.getValue());
+        }
+        if (result != null) {
+            gh.setStroke(new BasicStroke(2));
+            gh.setColor(Color.BLUE);
+            gh.draw(result);
+        }
+    }
+
+    /**
+     * Returns {@code true} if the shapes described by given iterators are equal.
+     * This is used for deciding if it is worth to bother the user with a request
+     * for pressing the "Next" button.
+     */
+    private static boolean equal(final PathIterator it1, final PathIterator it2) {
+        final float[] a1 = new float[6];
+        final float[] a2 = new float[6];
+        while (!it1.isDone()) {
+            if (it2.isDone()) return false;
+            final int code = it1.currentSegment(a1);
+            if (code != it2.currentSegment(a2)) {
+                return false;
+            }
+            int n;
+            switch (code) {
+                case PathIterator.SEG_MOVETO:
+                case PathIterator.SEG_LINETO:  n = 2; break;
+                case PathIterator.SEG_QUADTO:  n = 4; break;
+                case PathIterator.SEG_CUBICTO: n = 6; break;
+                case PathIterator.SEG_CLOSE:   n = 0; break;
+                default: throw new AssertionError(code);
+            }
+            while (--n >= 0) {
+                if (Float.floatToIntBits(a1[n]) != Float.floatToIntBits(a2[n])) {
+                    return false;
+                }
+            }
+            it1.next();
+            it2.next();
+        }
+        return it2.isDone();
+    }
+
+    /**
+     * Invoked after a row has been processed during the isoline generation.
+     * This is invoked from the main thread (<strong>not</strong> the Swing thread).
+     *
+     * @param  title      description of current state.
+     * @param  generator  new generator of isolines.
+     */
+    @Override
+    public void accept(final String title, final Isolines generator) {
+        final Map<PolylineStage, Path2D> paths = generator.toRawPath();
+        for (final Map.Entry<PolylineStage,Path2D> entry : paths.entrySet()) {
+            entry.getValue().transform(sourceToCanvas);
+        }
+        try {
+            final CountDownLatch c = new CountDownLatch(1);
+            EventQueue.invokeLater(() -> {
+                Rectangle b = null;
+                boolean unchanged = true;
+                for (final PolylineStage stage : PolylineStage.values()) {
+                    final Path2D current = isolines.get(stage);
+                    final Path2D update  = paths.get(stage);
+                    if (unchanged && current != update && !(current != null && update != null &&
+                            equal(current.getPathIterator(null), update.getPathIterator(null))))
+                    {
+                        unchanged = false;
+                    }
+                    if (update == null) {
+                        isolines.remove(stage);
+                    } else {
+                        isolines.put(stage, update);
+                        if (stage == PolylineStage.BUFFER) {
+                            b = update.getBounds();
+                            if (b.isEmpty()) {
+                                b = null;
+                            } else {
+                                b.x      -= PADDING;
+                                b.y      -= PADDING;
+                                b.width  += PADDING * 2;
+                                b.height += PADDING * 2;
+                            }
+                        }
+                    }
+                }
+                bounds = b;
+                if (unchanged) {
+                    stepTitle.setText(title + " (no change)");
+                    c.countDown();
+                } else {
+                    stepTitle.setText(title);
+                    repaint();
+                    assertNull(blocker);
+                    if (next.getModel().isPressed()) {
+                        c.countDown();
+                    } else {
+                        blocker = c;
+                        next.setEnabled(true);
+                    }
+                }
+            });
+            c.await();
+        } catch (InterruptedException  e) {
+            throw new AssertionError(e);            // Stop the test.
+        }
+    }
+
+    /**
+     * Invoked by Swing when user presses the "Next" button.
+     * This method resumes isoline computation.
+     *
+     * @param  event  ignored.
+     */
+    @Override
+    public void actionPerformed(final ActionEvent event) {
+        next.setEnabled(false);
+        if (blocker != null) {
+            blocker.countDown();
+            blocker = null;
+        }
+    }
+
+    /**
+     * Invoked when the "Next" button is kept pressed.
+     * The effect is to start the "fast forward" mode.
+     * This method shall be invoked in Swing thread.
+     *
+     * @param  event  ignored.
+     */
+    private void fastForward(final ActionEvent event) {
+        if (next.getModel().isPressed()) {
+            if (blocker != null) {
+                blocker.countDown();
+                blocker = null;
+            }
+        }
+    }
+
+    /**
+     * Invoked by Swing when the state of the "Next" button (pressed or not) changed.
+     * If the button is pressed one second without being released, then we enter a
+     * "fast forward" mode until the button is released.
+     *
+     * @param  event  ignored.
+     */
+    @Override
+    public void stateChanged(final ChangeEvent event) {
+        final ButtonModel m = (ButtonModel) event.getSource();
+        if (m.isPressed()) {
+            delayedNext.restart();
+        } else {
+            delayedNext.stop();
+        }
+    }
+}
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 db88eca..78e4a38 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
@@ -88,7 +88,7 @@
     org.apache.sis.internal.coverage.j2d.ScaledColorSpaceTest.class,
     org.apache.sis.internal.coverage.j2d.ColorizerTest.class,
     org.apache.sis.internal.coverage.j2d.SampleModelFactoryTest.class,
-    org.apache.sis.internal.processing.image.IsolinesTest.class,
+    org.apache.sis.internal.processing.isoline.IsolinesTest.class,
     org.apache.sis.image.DataTypeTest.class,
     org.apache.sis.image.PlanarImageTest.class,
     org.apache.sis.image.ComputedImageTest.class,
diff --git a/core/sis-metadata/src/main/java/org/apache/sis/metadata/iso/DefaultMetadata.java b/core/sis-metadata/src/main/java/org/apache/sis/metadata/iso/DefaultMetadata.java
index 3a21772..ed127c7 100644
--- a/core/sis-metadata/src/main/java/org/apache/sis/metadata/iso/DefaultMetadata.java
+++ b/core/sis-metadata/src/main/java/org/apache/sis/metadata/iso/DefaultMetadata.java
@@ -144,6 +144,7 @@
  * @since 0.3
  * @module
  */
+@SuppressWarnings("serial")                             // Fields are not statically typed as Serializable.
 @XmlType(name = "MD_Metadata_Type", propOrder = {
     // Attributes new in ISO 19115:2014
     "metadataIdentifier",
@@ -1646,8 +1647,8 @@
     }
 
     /**
-     * Invoked by JAXB {@link javax.xml.bind.Marshaller} after this object has been marshalled to
-     * XML. This method restores the locale to be used for XML marshalling to its previous value.
+     * Invoked by JAXB {@link javax.xml.bind.Marshaller} after this object has been marshalled to XML.
+     * This method restores the locale to be used for XML marshalling to its previous value.
      */
     @SuppressWarnings("unused")
     private void afterMarshal(final Marshaller marshaller) {
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/AbstractIdentifiedObject.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/AbstractIdentifiedObject.java
index 92341bf..4fd2a6f 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/AbstractIdentifiedObject.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/AbstractIdentifiedObject.java
@@ -181,6 +181,7 @@
      * @see #getName()
      * @see #getNames()
      */
+    @SuppressWarnings("serial")         // Not statically typed as Serializable.
     private ReferenceIdentifier name;
 
     /**
@@ -191,6 +192,7 @@
      * <p><b>Consider this field as final!</b>
      * This field is modified only at unmarshalling time by {@code Names.add(Identifier)}.</p>
      */
+    @SuppressWarnings("serial")         // Not statically typed as Serializable.
     private Collection<GenericName> alias;
 
     /**
@@ -203,6 +205,7 @@
      * @see #getIdentifiers()
      * @see #getIdentifier()
      */
+    @SuppressWarnings("serial")         // Not statically typed as Serializable.
     private Set<ReferenceIdentifier> identifiers;
 
     /**
@@ -213,6 +216,7 @@
      *
      * @see #getRemarks()
      */
+    @SuppressWarnings("serial")         // Not statically typed as Serializable.
     private InternationalString remarks;
 
     /**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/AbstractReferenceSystem.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/AbstractReferenceSystem.java
index 1f3457a..915c24b 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/AbstractReferenceSystem.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/AbstractReferenceSystem.java
@@ -82,6 +82,7 @@
      *
      * @see #getDomainOfValidity()
      */
+    @SuppressWarnings("serial")         // Not statically typed as Serializable.
     private Extent domainOfValidity;
 
     /**
@@ -93,6 +94,7 @@
      *
      * @see #getScope()
      */
+    @SuppressWarnings("serial")         // Not statically typed as Serializable.
     private InternationalString scope;
 
     /**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/AbstractCRS.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/AbstractCRS.java
index 656f309..078bcf3 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/AbstractCRS.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/AbstractCRS.java
@@ -118,6 +118,7 @@
      *
      * @see #getCoordinateSystem()
      */
+    @SuppressWarnings("serial")         // Not statically typed as Serializable.
     private CoordinateSystem coordinateSystem;
 
     /**
@@ -281,7 +282,7 @@
     final <T extends CoordinateSystem> T getCoordinateSystem(final Class<T> type) {
         final CoordinateSystem cs = coordinateSystem;
         if (type.isInstance(cs)) {
-            // Special case for AfficeCS: must ensure that the cs is not the CartesianCS subtype.
+            // Special case for AffineCS: must ensure that the cs is not the CartesianCS subtype.
             if (type != AffineCS.class || !(cs instanceof CartesianCS)) {
                 return (T) cs;
             }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/DefaultEngineeringCRS.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/DefaultEngineeringCRS.java
index 8d8965e..35eb3b4 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/DefaultEngineeringCRS.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/DefaultEngineeringCRS.java
@@ -19,15 +19,16 @@
 import java.util.Map;
 import javax.xml.bind.annotation.XmlType;
 import javax.xml.bind.annotation.XmlElement;
-import javax.xml.bind.annotation.XmlElements;
 import javax.xml.bind.annotation.XmlRootElement;
-import org.opengis.referencing.cs.CoordinateSystem;
+import javax.xml.bind.annotation.adapters.XmlJavaTypeAdapter;
+import org.opengis.referencing.cs.*;
 import org.opengis.referencing.crs.EngineeringCRS;
 import org.opengis.referencing.datum.EngineeringDatum;
 import org.apache.sis.referencing.cs.*;
 import org.apache.sis.referencing.AbstractReferenceSystem;
 import org.apache.sis.internal.metadata.MetadataUtilities;
 import org.apache.sis.internal.referencing.WKTKeywords;
+import org.apache.sis.internal.jaxb.referencing.CS_CoordinateSystem;
 import org.apache.sis.io.wkt.Formatter;
 
 import static org.apache.sis.util.ArgumentChecks.ensureNonNull;
@@ -60,7 +61,7 @@
  * in the javadoc, this condition holds if all components were created using only SIS factories and static constants.
  *
  * @author  Martin Desruisseaux (IRD, Geomatys)
- * @version 0.7
+ * @version 1.3
  *
  * @see org.apache.sis.referencing.datum.DefaultEngineeringDatum
  * @see org.apache.sis.referencing.factory.GeodeticAuthorityFactory#createEngineeringCRS(String)
@@ -69,7 +70,14 @@
  * @module
  */
 @XmlType(name = "EngineeringCRSType", propOrder = {
-    "coordinateSystem",
+    "abstractCS",
+    "affineCS",
+    "cartesianCS",
+    "cylindricalCS",
+    "linearCS",
+    "polarCS",
+    "sphericalCS",
+    "userDefinedCS",
     "datum"
 })
 @XmlRootElement(name = "EngineeringCRS")
@@ -87,6 +95,7 @@
      *
      * @see #getDatum()
      */
+    @SuppressWarnings("serial")         // Not statically typed as Serializable.
     private EngineeringDatum datum;
 
     /**
@@ -208,25 +217,6 @@
     }
 
     /**
-     * Returns the coordinate system.
-     *
-     * @return the coordinate system.
-     */
-    @Override
-    @XmlElements({
-        @XmlElement(name = "cartesianCS",   type = DefaultCartesianCS.class),
-        @XmlElement(name = "affineCS",      type = DefaultAffineCS.class),
-        @XmlElement(name = "cylindricalCS", type = DefaultCylindricalCS.class),
-        @XmlElement(name = "linearCS",      type = DefaultLinearCS.class),
-        @XmlElement(name = "polarCS",       type = DefaultPolarCS.class),
-        @XmlElement(name = "sphericalCS",   type = DefaultSphericalCS.class),
-        @XmlElement(name = "userDefinedCS", type = DefaultUserDefinedCS.class)
-    })
-    public CoordinateSystem getCoordinateSystem() {
-        return super.getCoordinateSystem();
-    }
-
-    /**
      * {@inheritDoc}
      *
      * @return {@inheritDoc}
@@ -302,10 +292,77 @@
 
     /**
      * Used by JAXB only (invoked by reflection).
+     * Only one of {@code getFooCS()} methods can return a non-null value.
      *
-     * @see #getCoordinateSystem()
+     * <div class="note"><b>Implementation note:</b>
+     * The usual way to handle {@code <xs:choice>} with JAXB is to annotate a single method like below:
+     *
+     * {@preformat java
+     *     &#64;Override
+     *     &#64;XmlElements({
+     *       &#64;XmlElement(name = "cartesianCS",   type = DefaultCartesianCS.class),
+     *       &#64;XmlElement(name = "affineCS",      type = DefaultAffineCS.class),
+     *       &#64;XmlElement(name = "cylindricalCS", type = DefaultCylindricalCS.class),
+     *       &#64;XmlElement(name = "linearCS",      type = DefaultLinearCS.class),
+     *       &#64;XmlElement(name = "polarCS",       type = DefaultPolarCS.class),
+     *       &#64;XmlElement(name = "sphericalCS",   type = DefaultSphericalCS.class),
+     *       &#64;XmlElement(name = "userDefinedCS", type = DefaultUserDefinedCS.class)
+     *     })
+     *     public CoordinateSystem getCoordinateSystem() {
+     *         return super.getCoordinateSystem();
+     *     }
+     * }
+     *
+     * However our attempts to apply this approach worked for {@code DefaultParameterValue} but not for this class:
+     * for an unknown reason, the unmarshalled CS object is empty.</div>
+     *
+     * @see <a href="http://issues.apache.org/jira/browse/SIS-166">SIS-166</a>
      */
-    private void setCoordinateSystem(final CoordinateSystem cs) {
-        setCoordinateSystem(null, cs);  // 'null' here means to infer the XML property name from the cs type.
+    @XmlElement(name="affineCS")      private AffineCS      getAffineCS()      {return getCoordinateSystem(AffineCS     .class);}
+    @XmlElement(name="cartesianCS")   private CartesianCS   getCartesianCS()   {return getCoordinateSystem(CartesianCS  .class);}
+    @XmlElement(name="cylindricalCS") private CylindricalCS getCylindricalCS() {return getCoordinateSystem(CylindricalCS.class);}
+    @XmlElement(name="linearCS")      private LinearCS      getLinearCS()      {return getCoordinateSystem(LinearCS     .class);}
+    @XmlElement(name="polarCS")       private PolarCS       getPolarCS()       {return getCoordinateSystem(PolarCS      .class);}
+    @XmlElement(name="sphericalCS")   private SphericalCS   getSphericalCS()   {return getCoordinateSystem(SphericalCS  .class);}
+    @XmlElement(name="userDefinedCS") private UserDefinedCS getUserDefinedCS() {return getCoordinateSystem(UserDefinedCS.class);}
+
+    /**
+     * Invoked by JAXB at unmarshalling time.
+     */
+    private void setAffineCS     (final AffineCS      cs) {super.setCoordinateSystem("affineCS",      cs);}
+    private void setCartesianCS  (final CartesianCS   cs) {super.setCoordinateSystem("cartesianCS",   cs);}
+    private void setCylindricalCS(final CylindricalCS cs) {super.setCoordinateSystem("cylindricalCS", cs);}
+    private void setLinearCS     (final LinearCS      cs) {super.setCoordinateSystem("linearCS",      cs);}
+    private void setPolarCS      (final PolarCS       cs) {super.setCoordinateSystem("polarCS",       cs);}
+    private void setSphericalCS  (final SphericalCS   cs) {super.setCoordinateSystem("sphericalCS",   cs);}
+    private void setUserDefinedCS(final UserDefinedCS cs) {super.setCoordinateSystem("userDefinedCS", cs);}
+
+    /**
+     * The types for which a specialized method exists.
+     * Not including {@link CartesianCS}, because this case is already covered by {@link AffineCS}.
+     */
+    private static final Class<?>[] SPECIALIZED_TYPES = {
+        AffineCS.class, SphericalCS.class, CylindricalCS.class, PolarCS.class, LinearCS.class, UserDefinedCS.class
+    };
+
+    /**
+     * Returns the coordinate system if it is not an instance of any of the types handled by specialized methods.
+     * It is the case of {@link EllipsoidalCS}, {@link VerticalCS}, {@link TimeCS} and {@link ParametricCS}.
+     */
+    @XmlElement(name = "coordinateSystem", required = true)
+    @XmlJavaTypeAdapter(CS_CoordinateSystem.class)
+    private CoordinateSystem getAbstractCS() {
+        final CoordinateSystem cs = getCoordinateSystem();
+        for (final Class<?> t : SPECIALIZED_TYPES) {
+            if (t.isInstance(cs)) return null;
+        }
+        return cs;
+    }
+
+    /**
+     * Used by JAXB only (invoked by reflection).
+     */
+    private void setAbstractCS(final CoordinateSystem cs) {
+        setCoordinateSystem(null, cs);      // `null` here means to infer the XML property name from the cs type.
     }
 }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/DefaultGeodeticCRS.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/DefaultGeodeticCRS.java
index a36bc1a..9aaa2c9 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/DefaultGeodeticCRS.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/DefaultGeodeticCRS.java
@@ -324,10 +324,8 @@
      *     }
      * }
      *
-     * However our attempts to apply this approach worked for {@link DefaultEngineeringCRS} but not for this class:
-     * for an unknown reason, the unmarshalled CS object is empty. Maybe this is because the covariant return type
-     * in the {@link DefaultGeographicCRS} ({@code EllipsoidCS} instead of {@code CoordinateSystem} in above code)
-     * is causing confusion.</div>
+     * However our attempts to apply this approach worked for {@code DefaultParameterValue} but not for this class:
+     * for an unknown reason, the unmarshalled CS object is empty.</div>
      *
      * @see <a href="http://issues.apache.org/jira/browse/SIS-166">SIS-166</a>
      */
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/DefaultImageCRS.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/DefaultImageCRS.java
index d307556..265e9ea 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/DefaultImageCRS.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/crs/DefaultImageCRS.java
@@ -306,9 +306,8 @@
      *     }
      * }
      *
-     * However our attempts to apply this approach worked for {@link DefaultEngineeringCRS} but not for this class:
-     * for an unknown reason, the unmarshalled CS object is empty. Maybe this is because the covariant return type
-     * ({@code AffineCS} instead of {@code CoordinateSystem} in above code) is causing confusion.</div>
+     * However our attempts to apply this approach worked for {@code DefaultParameterValue} but not for this class:
+     * for an unknown reason, the unmarshalled CS object is empty.</div>
      *
      * @see <a href="http://issues.apache.org/jira/browse/SIS-166">SIS-166</a>
      */
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/cs/AbstractCS.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/cs/AbstractCS.java
index cd5fef4..31c2c7e 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/cs/AbstractCS.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/cs/AbstractCS.java
@@ -117,6 +117,7 @@
      *
      * @see #getAxis(int)
      */
+    @SuppressWarnings("serial")         // Not statically typed as Serializable.
     private CoordinateSystemAxis[] axes;
 
     /**
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/crs/DefaultEngineeringCRSTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/crs/DefaultEngineeringCRSTest.java
index e130975..0748ff2 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/crs/DefaultEngineeringCRSTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/crs/DefaultEngineeringCRSTest.java
@@ -36,7 +36,7 @@
  * Tests {@link DefaultEngineeringCRS}.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.3
  * @since   0.6
  * @module
  */
@@ -117,22 +117,24 @@
         assertXmlEquals(
                 "<gml:EngineeringCRS xmlns:gml=\"" + Namespaces.GML + "\">\n" +
                 "  <gml:name>A construction site CRS</gml:name>\n" +
-                "  <gml:cartesianCS gml:id=\"Cartesian2D\">\n" +
-                "    <gml:name>Cartesian 2D</gml:name>\n" +
-                "    <gml:axis>\n" +
-                "      <gml:CoordinateSystemAxis uom=\"urn:ogc:def:uom:EPSG::9001\" gml:id=\"x\">\n" +
-                "        <gml:name>x</gml:name>\n" +
-                "        <gml:axisAbbrev>x</gml:axisAbbrev>\n" +
-                "        <gml:axisDirection codeSpace=\"EPSG\">east</gml:axisDirection>\n" +
-                "      </gml:CoordinateSystemAxis>\n" +
-                "    </gml:axis>\n" +
-                "    <gml:axis>\n" +
-                "      <gml:CoordinateSystemAxis uom=\"urn:ogc:def:uom:EPSG::9001\" gml:id=\"y\">\n" +
-                "        <gml:name>y</gml:name>\n" +
-                "        <gml:axisAbbrev>y</gml:axisAbbrev>\n" +
-                "        <gml:axisDirection codeSpace=\"EPSG\">north</gml:axisDirection>\n" +
-                "      </gml:CoordinateSystemAxis>\n" +
-                "    </gml:axis>\n" +
+                "  <gml:cartesianCS>\n" +
+                "    <gml:CartesianCS gml:id=\"Cartesian2D\">\n" +
+                "      <gml:name>Cartesian 2D</gml:name>\n" +
+                "      <gml:axis>\n" +
+                "        <gml:CoordinateSystemAxis uom=\"urn:ogc:def:uom:EPSG::9001\" gml:id=\"x\">\n" +
+                "          <gml:name>x</gml:name>\n" +
+                "          <gml:axisAbbrev>x</gml:axisAbbrev>\n" +
+                "          <gml:axisDirection codeSpace=\"EPSG\">east</gml:axisDirection>\n" +
+                "        </gml:CoordinateSystemAxis>\n" +
+                "      </gml:axis>\n" +
+                "      <gml:axis>\n" +
+                "        <gml:CoordinateSystemAxis uom=\"urn:ogc:def:uom:EPSG::9001\" gml:id=\"y\">\n" +
+                "          <gml:name>y</gml:name>\n" +
+                "          <gml:axisAbbrev>y</gml:axisAbbrev>\n" +
+                "          <gml:axisDirection codeSpace=\"EPSG\">north</gml:axisDirection>\n" +
+                "        </gml:CoordinateSystemAxis>\n" +
+                "      </gml:axis>\n" +
+                "    </gml:CartesianCS>\n" +
                 "  </gml:cartesianCS>\n" +
                 "  <gml:engineeringDatum>\n" +
                 "    <gml:EngineeringDatum gml:id=\"P1\">\n" +
@@ -167,37 +169,39 @@
         assertXmlEquals(
                 "<gml:EngineeringCRS xmlns:gml=\"" + Namespaces.GML + "\">\n" +
                 "  <gml:name>A spherical CRS</gml:name>\n" +
-                "  <gml:sphericalCS gml:id=\"Spherical\">\n" +
-                "    <gml:name>Spherical</gml:name>\n" +
-                "    <gml:axis>\n" +
-                "      <gml:CoordinateSystemAxis uom=\"urn:ogc:def:uom:EPSG::9122\" gml:id=\"SphericalLatitude\">\n" +
-                "        <gml:name>Spherical latitude</gml:name>\n" +
-                "        <gml:axisAbbrev>Ω</gml:axisAbbrev>\n" +
-                "        <gml:axisDirection codeSpace=\"EPSG\">north</gml:axisDirection>\n" +
-                "        <gml:minimumValue>-90.0</gml:minimumValue>\n" +
-                "        <gml:maximumValue>90.0</gml:maximumValue>\n" +
-                "        <gml:rangeMeaning codeSpace=\"EPSG\">exact</gml:rangeMeaning>\n" +
-                "      </gml:CoordinateSystemAxis>\n" +
-                "    </gml:axis>\n" +
-                "    <gml:axis>\n" +
-                "      <gml:CoordinateSystemAxis uom=\"urn:ogc:def:uom:EPSG::9122\" gml:id=\"SphericalLongitude\">\n" +
-                "        <gml:name>Spherical longitude</gml:name>\n" +
-                "        <gml:axisAbbrev>θ</gml:axisAbbrev>\n" +
-                "        <gml:axisDirection codeSpace=\"EPSG\">east</gml:axisDirection>\n" +
-                "        <gml:minimumValue>-180.0</gml:minimumValue>\n" +
-                "        <gml:maximumValue>180.0</gml:maximumValue>\n" +
-                "        <gml:rangeMeaning codeSpace=\"EPSG\">wraparound</gml:rangeMeaning>\n" +
-                "      </gml:CoordinateSystemAxis>\n" +
-                "    </gml:axis>\n" +
-                "    <gml:axis>\n" +
-                "      <gml:CoordinateSystemAxis uom=\"urn:ogc:def:uom:EPSG::9001\" gml:id=\"GeocentricRadius\">\n" +
-                "        <gml:name>Geocentric radius</gml:name>\n" +
-                "        <gml:axisAbbrev>r</gml:axisAbbrev>\n" +
-                "        <gml:axisDirection codeSpace=\"EPSG\">up</gml:axisDirection>\n" +
-                "        <gml:minimumValue>0.0</gml:minimumValue>\n" +
-                "        <gml:rangeMeaning codeSpace=\"EPSG\">exact</gml:rangeMeaning>\n" +
-                "      </gml:CoordinateSystemAxis>\n" +
-                "    </gml:axis>\n" +
+                "  <gml:sphericalCS>\n" +
+                "    <gml:SphericalCS gml:id=\"Spherical\">\n" +
+                "      <gml:name>Spherical</gml:name>\n" +
+                "      <gml:axis>\n" +
+                "        <gml:CoordinateSystemAxis uom=\"urn:ogc:def:uom:EPSG::9122\" gml:id=\"SphericalLatitude\">\n" +
+                "          <gml:name>Spherical latitude</gml:name>\n" +
+                "          <gml:axisAbbrev>Ω</gml:axisAbbrev>\n" +
+                "          <gml:axisDirection codeSpace=\"EPSG\">north</gml:axisDirection>\n" +
+                "          <gml:minimumValue>-90.0</gml:minimumValue>\n" +
+                "          <gml:maximumValue>90.0</gml:maximumValue>\n" +
+                "          <gml:rangeMeaning codeSpace=\"EPSG\">exact</gml:rangeMeaning>\n" +
+                "        </gml:CoordinateSystemAxis>\n" +
+                "      </gml:axis>\n" +
+                "      <gml:axis>\n" +
+                "        <gml:CoordinateSystemAxis uom=\"urn:ogc:def:uom:EPSG::9122\" gml:id=\"SphericalLongitude\">\n" +
+                "          <gml:name>Spherical longitude</gml:name>\n" +
+                "          <gml:axisAbbrev>θ</gml:axisAbbrev>\n" +
+                "          <gml:axisDirection codeSpace=\"EPSG\">east</gml:axisDirection>\n" +
+                "          <gml:minimumValue>-180.0</gml:minimumValue>\n" +
+                "          <gml:maximumValue>180.0</gml:maximumValue>\n" +
+                "          <gml:rangeMeaning codeSpace=\"EPSG\">wraparound</gml:rangeMeaning>\n" +
+                "        </gml:CoordinateSystemAxis>\n" +
+                "      </gml:axis>\n" +
+                "      <gml:axis>\n" +
+                "        <gml:CoordinateSystemAxis uom=\"urn:ogc:def:uom:EPSG::9001\" gml:id=\"GeocentricRadius\">\n" +
+                "          <gml:name>Geocentric radius</gml:name>\n" +
+                "          <gml:axisAbbrev>r</gml:axisAbbrev>\n" +
+                "          <gml:axisDirection codeSpace=\"EPSG\">up</gml:axisDirection>\n" +
+                "          <gml:minimumValue>0.0</gml:minimumValue>\n" +
+                "          <gml:rangeMeaning codeSpace=\"EPSG\">exact</gml:rangeMeaning>\n" +
+                "        </gml:CoordinateSystemAxis>\n" +
+                "      </gml:axis>\n" +
+                "    </gml:SphericalCS>\n" +
                 "  </gml:sphericalCS>\n" +
                 "  <gml:engineeringDatum>\n" +
                 "    <gml:EngineeringDatum gml:id=\"Centre\">\n" +