sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 07/07: Take in account the band subset specified by the `range` argument in a call to `GridCoverage.render(…)` on a GeoTIFF image.
Date Thu, 29 Jul 2021 13:07:22 GMT
This is an automated email from the ASF dual-hosted git repository.

desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git

commit d53b8a5db8e77371dfb32beebaabda80defe2f0d
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed Jul 28 18:26:49 2021 +0200

    Take in account the band subset specified by the `range` argument in a call to `GridCoverage.render(…)` on a GeoTIFF image.
---
 .../org/apache/sis/internal/geotiff/Inflater.java  |  52 +++---
 .../apache/sis/internal/geotiff/Uncompressed.java  |  53 +++---
 .../sis/storage/geotiff/CompressedSubset.java      | 106 ++++++++++--
 .../org/apache/sis/storage/geotiff/DataCube.java   |  12 +-
 .../org/apache/sis/storage/geotiff/DataSubset.java | 107 ++++++------
 .../sis/storage/geotiff/ImageFileDirectory.java    |  10 ++
 .../sis/internal/storage/AbstractGridResource.java |  18 +-
 .../sis/internal/storage/TiledGridCoverage.java    |  15 ++
 .../sis/internal/storage/TiledGridResource.java    |  14 +-
 .../internal/storage/io/HyperRectangleReader.java  |  11 ++
 .../sis/test/storage/CoverageReadConsistency.java  | 182 ++++++++++++++++-----
 11 files changed, 431 insertions(+), 149 deletions(-)

diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Inflater.java b/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Inflater.java
index e5fa40c..9cd502f 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Inflater.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Inflater.java
@@ -41,47 +41,59 @@ public abstract class Inflater {
 
     /**
      * Number of chunk per row, as a strictly positive integer.
-     * A chunk is a pixel, except if we can optimize by reading the whole row as a single chunk.
+     * See {@link #samplesPerChunk} for more details about what is a "chunk".
      */
     protected final int chunksPerRow;
 
     /**
      * Number of sample values per chunk, as a strictly positive integer.
-     * A chunk is a pixel, except if we can optimize by reading the whole row as a single chunk.
+     * A chunk can be:
+     * <ul>
+     *   <li>A sample value ({@code samplesPerChunk} = 1).</li>
+     *   <li>A pixel with one sample value per band ({@code samplesPerChunk} = number of bands).</li>
+     *   <li>A full row (optimization when it is possible to read the row in a single I/O method call).</li>
+     * </ul>
      */
     protected final int samplesPerChunk;
 
     /**
-     * Number of sample values to skip between pixels. Positive but often zero.
+     * Number of sample values to skip between chunks on the same row, or {@code null} if none.
+     * If non-null then after reading the chunk at zero-based index <var>x</var>, inflater shall skip
+     * {@code skipAfterChunks[x % skipAfterChunks.length]} sample values before to read the next chunk.
+     * The <var>x</var> index is reset to zero at the beginning of every new row.
      */
-    protected final int interpixels;
+    protected final int[] skipAfterChunks;
 
     /**
      * Creates a new instance.
      *
-     * @param  input            the source of data to decompress.
-     * @param  pixelsPerRow     number of pixels per row. Must be strictly positive.
-     * @param  samplesPerPixel  number of sample values per pixel. Must be strictly positive.
-     * @param  interpixels      number of sample values to skip between pixels. May be zero.
+     * @param  input              the source of data to decompress.
+     * @param  elementsPerRow     number of elements (usually pixels) per row. Must be strictly positive.
+     * @param  samplesPerElement  number of sample values per element (usually pixel). Must be strictly positive.
+     * @param  skipAfterElements  number of sample values to skip between elements (pixels). May be empty or null.
      */
-    protected Inflater(final ChannelDataInput input, final int pixelsPerRow, final int samplesPerPixel, final int interpixels) {
-        ensureStrictlyPositive("pixelsPerRow",    pixelsPerRow);
-        ensureStrictlyPositive("samplesPerPixel", samplesPerPixel);
-        ensurePositive        ("interpixels",     interpixels);
-        ensureNonNull         ("input",           input);
-        if (interpixels == 0) {
-            chunksPerRow    = 1;
-            samplesPerChunk = Math.multiplyExact(samplesPerPixel, pixelsPerRow);
+    protected Inflater(final ChannelDataInput input, final int elementsPerRow, final int samplesPerElement, final int[] skipAfterElements) {
+        ensureNonNull("input", input);
+        ensureStrictlyPositive("elementsPerRow",    elementsPerRow);
+        ensureStrictlyPositive("samplesPerElement", samplesPerElement);
+        this.input = input;
+        skipAfterChunks = skipAfterElements;
+        if (skipAfterElements != null) {
+            for (int i=0; i<skipAfterElements.length; i++) {
+                ensurePositive("skipAfterElements", skipAfterElements[i]);
+            }
+            chunksPerRow    = elementsPerRow;
+            samplesPerChunk = samplesPerElement;
         } else {
-            chunksPerRow    = pixelsPerRow;
-            samplesPerChunk = samplesPerPixel;
+            chunksPerRow    = 1;
+            samplesPerChunk = Math.multiplyExact(samplesPerElement, elementsPerRow);
         }
-        this.interpixels = interpixels;
-        this.input       = input;
     }
 
     /**
      * Reads a row of sample values and stores them in the target buffer.
+     * This is not a complete row; caller may invoke {@link #skip(long)}
+     * for ignoring leading and trailing sample values.
      *
      * @throws IOException if an error occurred while reading the input channel.
      */
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Uncompressed.java b/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Uncompressed.java
index 3cb3ce5..9b695a6 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Uncompressed.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Uncompressed.java
@@ -26,6 +26,9 @@ import org.apache.sis.util.Classes;
 
 /**
  * A pseudo-inflater which copy values unchanged.
+ * This implementation is useful for handling more complex subsampling
+ * than what {@link org.apache.sis.internal.storage.io.HyperRectangleReader} can handle.
+ * It is also useful for testing purpose.
  *
  * @author  Martin Desruisseaux (Geomatys)
  * @version 1.1
@@ -52,8 +55,10 @@ public abstract class Uncompressed extends Inflater {
     /**
      * For constructors in inner classes.
      */
-    private Uncompressed(ChannelDataInput input, long start, int pixelsPerRow, int samplesPerPixel, int interpixels, final int sampleSize) {
-        super(input, pixelsPerRow, samplesPerPixel, interpixels);
+    private Uncompressed(final ChannelDataInput input, final long start, final int elementsPerRow,
+                         final int samplesPerElement, final int[] skipAfterElements, final int sampleSize)
+    {
+        super(input, elementsPerRow, samplesPerElement, skipAfterElements);
         this.streamPosition = start;
         this.sampleSize = sampleSize;
     }
@@ -61,25 +66,26 @@ public abstract class Uncompressed extends Inflater {
     /**
      * Creates a new instance.
      *
-     * @param  input            the source of data to decompress.
-     * @param  start            stream position where to start reading.
-     * @param  pixelsPerRow     number of pixels per row. Must be strictly positive.
-     * @param  samplesPerPixel  number of sample values per pixel. Must be strictly positive.
-     * @param  interpixels      number of sample values to skip between pixels. May be zero.
-     * @param  target           where to store sample values.
+     * @param  input   the source of data to decompress.
+     * @param  start   stream position where to start reading.
+     * @param  count   number of elements (usually pixels) per row. Must be strictly positive.
+     * @param  size    number of sample values per element (usually pixel). Must be strictly positive.
+     * @param  skips   number of sample values to skip between elements (pixels). May be empty or null.
+     * @param  target  where to store sample values.
      * @return the inflater for the given targe type.
      * @throws IllegalArgumentException if the buffer type is not recognized.
      */
     public static Uncompressed create(final ChannelDataInput input, final long start,
-            final int pixelsPerRow, final int samplesPerPixel, final int interpixels, final Buffer target)
+            final int count, final int size, final int[] skips, final Buffer target)
     {
-        if (target instanceof ByteBuffer) return new Bytes(input, start, pixelsPerRow, samplesPerPixel, interpixels, (ByteBuffer) target);
+        if (target instanceof ByteBuffer) return new Bytes(input, start, count, size, skips, (ByteBuffer) target);
         throw new IllegalArgumentException(Errors.format(Errors.Keys.UnsupportedType_1, Classes.getClass(target)));
     }
 
     /**
      * Reads a row of sample values and stores them in the target buffer.
-     * Subclasses must override and invoke {@code super.uncompress()} before to do the actual reading.
+     * Subclasses must override this method and invoke {@code super.uncompress()}
+     * before to do the actual reading.
      */
     @Override
     public void uncompressRow() throws IOException {
@@ -108,23 +114,22 @@ public abstract class Uncompressed extends Inflater {
     }
 
     /**
-     * Inflater when the values to read and store are bytes.
+     * Inflater for sample values stored as bytes.
      */
     private static final class Bytes extends Uncompressed {
         /** Where to copy the values that we will read. */
         private final ByteBuffer target;
 
         /** Creates a new inflater which will write in the given buffer. */
-        Bytes(ChannelDataInput input, long start, int pixelsPerRow, int samplesPerPixel, int interpixels, ByteBuffer target) {
-            super(input, start, pixelsPerRow, samplesPerPixel, interpixels, Byte.BYTES);
+        Bytes(ChannelDataInput input, long start, int count, int size, int[] skips, ByteBuffer target) {
+            super(input, start, count, size, skips, Byte.BYTES);
             this.target = target;
         }
 
-        /**
-         * Reads and decompress a row of sample values.
-         */
+        /** Reads and decompress a row of sample values. */
         @Override public void uncompressRow() throws IOException {
             super.uncompressRow();
+            int ip = 0;
             for (int i = chunksPerRow; --i > 0;) {      // (chunksPerRow - 1) iterations.
                 int n = samplesPerChunk;
                 do target.put(input.readByte());
@@ -134,14 +139,18 @@ public abstract class Uncompressed extends Inflater {
                  * We invoke `readByte()` in a loop instead of invoking `skip` because if
                  * the number of bytes to skip is small, this is more efficient.
                  */
-                for (n = interpixels; --n >= 0;) {
-                    input.readByte();
+                if (skipAfterChunks != null) {
+                    for (n = skipAfterChunks[ip]; --n >= 0;) {
+                        input.readByte();
+                    }
+                    if (++ip >= skipAfterChunks.length) ip = 0;
                 }
             }
             /*
-             * Read the last element that was not read in first above `for` loop, but without
-             * skipping `interpixels` values after it. This is necessary for avoiding EOF if
-             * the last pixel to read is in the last column of the tile.
+             * Read the last element that was not read in above `for` loop,
+             * but without skipping `skipAfterElements` sample values after.
+             * This is necessary for avoiding EOF if the last pixel to read
+             * is in the last column of the tile.
              */
             int n = samplesPerChunk;
             do target.put(input.readByte());
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/CompressedSubset.java b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/CompressedSubset.java
index d071da6..1b36eb1 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/CompressedSubset.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/CompressedSubset.java
@@ -41,11 +41,57 @@ import static org.apache.sis.internal.jdk9.JDK9.multiplyFull;
  */
 final class CompressedSubset extends DataSubset {
     /**
-     * Number of sample values to skip for moving to the next row of a tile.
+     * Number of sample values to skip for moving to the next row of a tile in the GeoTIFF file.
+     * This is not necessarily the same scanline stride than for the tiles created by this class.
      */
     private final long scanlineStride;
 
     /**
+     * Number of sample values to skip before to read the first value of the first pixel in a row.
+     * The first pixel is at column index 0; subsampling offset is not included in this calculation.
+     */
+    private final int beforeFirstBand;
+
+    /**
+     * Number of sample values to skip for reaching end-of-row after reading the last value of the
+     * <em>first</em> pixel in a row. For computing the actual number of sample values to skip,
+     * the number of sample values read of skipped before the last pixel must be subtracted.
+     */
+    private final int afterLastBand;
+
+    /**
+     * Number of sample values to skip after an element has been read, or {@code null} if none.
+     * An element is a sample value or a complete pixel, depending on {@link #samplesPerElement}.
+     * The {@code skipAfterElements} array is used in a cyclic way:
+     *
+     * <ul>
+     *   <li>Skip {@code skipAfterElements[0]} sample values between element 0 and element 1.</li>
+     *   <li>Skip {@code skipAfterElements[1]} sample values between element 1 and element 2.</li>
+     *   <li><i>etc</i>. When we reach the array end, continue at {@code skipAfterElements[0]}.</li>
+     *   <li>When we start a new row, unconditionally restart at {@code skipAfterElements[0]}.</li>
+     * </ul>
+     *
+     * More generally, skip {@code skipAfterElements[x % skipAfterElements.length]} sample values
+     * after element at the zero-based column index <var>x</var>.
+     */
+    private final int[] skipAfterElements;
+
+    /**
+     * Number of sample values that compose an element (pixel or sample) in the GeoTIFF file.
+     * The value of this field can be:
+     *
+     * <ul>
+     *   <li>1 if an element is a sample value.</li>
+     *   <li>{@link #sourcePixelStride} if an element is a full pixel.</li>
+     *   <li>Any intermediate value if some optimizations have been applied,
+     *       for example for taking advantage of consecutive indices in {@link #selectedBands}.</li>
+     * </ul>
+     *
+     * This value shall always be a divisor of {@link #targetPixelStride}.
+     */
+    private final int samplesPerElement;
+
+    /**
      * Creates a new data subset. All parameters should have been validated
      * by {@link ImageFileDirectory#validateMandatoryTags()} before this call.
      * This constructor should be invoked inside a synchronized block.
@@ -57,7 +103,33 @@ final class CompressedSubset extends DataSubset {
      */
     CompressedSubset(final DataCube source, final TiledGridResource.Subset subset) throws DataStoreException {
         super(source, subset);
-        scanlineStride = multiplyFull(getTileSize(0), numInterleaved);
+        scanlineStride = multiplyFull(getTileSize(0), sourcePixelStride);
+        final int between = sourcePixelStride * (getSubsampling(0) - 1);
+        int afterLastBand = sourcePixelStride * (getTileSize(0) - 1);
+        if (selectedBands != null && sourcePixelStride > 1) {
+            final int n = selectedBands.length;
+            skipAfterElements = new int[n];
+            int b = sourcePixelStride;
+            for (int i = n; --i >= 0;) {
+                // Number of sample values to skip after each band.
+                skipAfterElements[i] = b - (b = selectedBands[i]) - 1;
+            }
+            beforeFirstBand         = b;
+            afterLastBand          += skipAfterElements[n-1];       // Add trailing bands that were left unread.
+            skipAfterElements[n-1] += between + beforeFirstBand;    // Add pixels skipped by subsampling and move to first band.
+            samplesPerElement       = 1;
+            // TODO: we could optimize if we find that all sample values to read are consecutive.
+        } else {
+            skipAfterElements = (between != 0) ? new int[] {between} : null;
+            samplesPerElement = sourcePixelStride;
+            beforeFirstBand   = 0;
+        }
+        this.afterLastBand = afterLastBand;
+        /*
+         * Invariant documented in Javadoc.
+         * Calculation correctness depends on this condition.
+         */
+        assert targetPixelStride % samplesPerElement == 0 : samplesPerElement;
     }
 
     /**
@@ -79,39 +151,43 @@ final class CompressedSubset extends DataSubset {
      * @param  upper        (<var>x</var>, <var>y</var>) coordinates after the last pixel to read relative to the tile.
      * @param  subsampling  (<var>sx</var>, <var>sy</var>) subsampling factors.
      * @param  location     pixel coordinates in the upper-left corner of the tile to return.
-     * @return image decoded from the GeoTIFF file.
+     * @return a single tile decoded from the GeoTIFF file.
      */
     @Override
     WritableRaster readSlice(final long[] offsets, final long[] byteCounts, final long[] lower, final long[] upper,
                              final int[] subsampling, final Point location) throws IOException, DataStoreException
     {
-        final DataType type   = getDataType();
-        final int      width  = pixelCount(lower, upper, subsampling, 0);
-        final int      height = pixelCount(lower, upper, subsampling, 1);
-        final int      skipY  = subsampling[1] - 1;
-        final int      skipX  = numInterleaved * (subsampling[0] - 1);
-        final long     head   = numInterleaved * lower[0];
-        final long     tail   = numInterleaved * (getTileSize(0) - width*subsampling[0]) + skipX - head;
+        final DataType type       = getDataType();
+        final int  width          = pixelCount(lower, upper, subsampling, 0);
+        final int  height         = pixelCount(lower, upper, subsampling, 1);
+        final int  elementsPerRow = width * (targetPixelStride / samplesPerElement);
+        final int  betweenRows    = subsampling[1] - 1;
+        final long head           = beforeFirstBand + sourcePixelStride * (lower[0]);
+        final long tail           = afterLastBand   - sourcePixelStride * (lower[0] + (width-1)*subsampling[0]);
         /*
          * `head` and `tail` are the number of sample values to skip at the beginning and end of each row.
-         * Then `skipX` is the number of sample values to skip between each pixel.
+         * `betweenPixels` is the number of sample values to skip between each pixel, but the actual skips
+         * are more complicated if only a subset of the bands are read. The actual number of sample values
+         * to skip between "elements" is detailed by `skipAfterElements`.
          */
-        final Buffer[] banks  = new Buffer[numBanks];
+        final Buffer[] banks = new Buffer[numBanks];
         for (int b=0; b<numBanks; b++) {
             final Buffer   bank = RasterFactory.createBuffer(type, capacity);
-            final Inflater algo = Uncompressed.create(reader().input, offsets[b], width, numInterleaved, skipX, bank);
+            final Inflater algo = Uncompressed.create(reader().input, offsets[b], elementsPerRow, samplesPerElement, skipAfterElements, bank);
             for (long y = lower[1]; --y >= 0;) {
                 algo.skip(scanlineStride);
+                // TODO: after we finished to implement decompression algorithms,
+                // revisit if we can safely replace the loop by a multiplication.
             }
             for (int y = height; --y > 0;) {        // (height - 1) iterations.
                 algo.skip(head);
                 algo.uncompressRow();
                 algo.skip(tail);
-                for (int j=skipY; --j>=0;) {
+                for (int j=betweenRows; --j>=0;) {
                     algo.skip(scanlineStride);
                 }
             }
-            algo.skip(head);                        // Last iteration without last `skip(…)` calls.
+            algo.skip(head);                        // Last iteration without the trailing `skip(…)` calls.
             algo.uncompressRow();
             fillRemainingRows(bank.flip());
             banks[b] = bank;
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataCube.java b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataCube.java
index 8fd4870..23f2e66 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataCube.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataCube.java
@@ -75,6 +75,16 @@ abstract class DataCube extends TiledGridResource implements ResourceOnFileSyste
     }
 
     /**
+     * Returns the number of components per pixel in the image stored in GeoTIFF file.
+     * This the same value than the one returned by {@code getSampleModel().getNumBands()},
+     * and is also the size of the collection returned by {@link #getSampleDimensions()}.
+     *
+     * @see #getSampleModel()
+     * @see SampleModel#getNumBands()
+     */
+    abstract int getNumBands();
+
+    /**
      * Returns the total number of tiles. This is used for computing the stride between a
      * band and the next band in {@link #tileOffsets} and {@link #tileByteCounts} vectors.
      */
@@ -123,7 +133,7 @@ abstract class DataCube extends TiledGridResource implements ResourceOnFileSyste
                 }
                 switch (compression) {
                     case NONE: {
-                        if (subset.hasSubsampling(0) && isInterleaved()) {
+                        if (subset.hasBandSubset() || (subset.hasSubsampling(0) && isInterleaved())) {
                             coverage = new CompressedSubset(this, subset);
                         } else {
                             coverage = new DataSubset(this, subset);
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataSubset.java b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataSubset.java
index 562af30..313b262 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataSubset.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataSubset.java
@@ -100,21 +100,32 @@ class DataSubset extends TiledGridCoverage implements Localized {
     private final int numTiles;
 
     /**
-     * Number of banks in the image data buffer.
+     * Number of banks retained for the target image data buffer.
      * This is equal to the number of bands only for planar images, and 1 in all other cases.
+     * If the user asked to read only a subset of the bands, then "number of bands" in above
+     * sentence is the {@linkplain #selectedBands subset} size.
      */
     protected final int numBanks;
 
     /**
-     * Number of interleaved sample values in a pixel. For planar images, this is equal to 1.
-     * For interleaved sample model, this is equal to the number of bands. This value is often
-     * equal to the {@linkplain java.awt.image.ComponentSampleModel#getPixelStride() pixel stride}.
+     * Number of interleaved sample values in a pixel in the GeoTIFF file (ignoring band subset).
+     * For planar images (banded sample model), this is equal to 1. For pixel interleaved image,
+     * this is equal to the number of bands in the original image.
+     *
+     * @see java.awt.image.ComponentSampleModel#getPixelStride()
      */
-    protected final int numInterleaved;
+    protected final int sourcePixelStride;
+
+    /**
+     * Number of interleaved sample values in a pixel of the image to load in memory.
+     * This is similar to {@link #sourcePixelStride}, but taking in account the number
+     * of bands requested by user at reading time.
+     */
+    protected final int targetPixelStride;
 
     /**
      * Number of sample values in a bank (not necessarily a band).
-     * This is tile width × height × {@link #numInterleaved}.
+     * This is tile width × height × {@code targetPixelStride}.
      */
     protected final int capacity;
 
@@ -139,19 +150,24 @@ class DataSubset extends TiledGridCoverage implements Localized {
          * "Banks" (in `java.awt.image.DataBuffer` sense) are synonymous to "bands" for planar image only.
          * Otherwise there is only one bank no matter the amount of bands. Each bank will be read separately.
          */
+        final int maxBank;
         if (model instanceof BandedSampleModel) {
             numBanks = model.getNumBands();
-            numInterleaved = 1;
+            sourcePixelStride = targetPixelStride = 1;
+            maxBank = (selectedBands != null) ? selectedBands[selectedBands.length - 1] : numBanks - 1;
+            // Note: `selectedBands` is in strictly increasing order, so taking the last value is okay.
         } else {
+            maxBank  = 0;
             numBanks = 1;
-            numInterleaved = model.getNumBands();
+            sourcePixelStride = source.getNumBands();
+            targetPixelStride = model .getNumBands();
         }
         final int n = tileOffsets.size();
-        if (numBanks > n / numTiles) {
+        if (maxBank >= n / numTiles) {
             throw new DataStoreContentException(source.reader.errors().getString(
-                    Errors.Keys.TooFewCollectionElements_3, "tileOffsets", numBanks * numTiles, n));
+                    Errors.Keys.TooFewCollectionElements_3, "tileOffsets", (maxBank + 1) * numTiles, n));
         }
-        capacity = multiplyExact(multiplyExact(model.getWidth(), model.getHeight()), numInterleaved);
+        capacity = multiplyExact(multiplyExact(model.getWidth(), model.getHeight()), targetPixelStride);
     }
 
     /**
@@ -192,8 +208,8 @@ class DataSubset extends TiledGridCoverage implements Localized {
         /**
          * Value of {@link DataSubset#tileOffsets} at index {@link #indexInTileVector}.
          * If pixel data are stored in different planes ("banks" in Java2D terminology),
-         * this is the smallest value of all banks.
-         * This is the value that we want in increasing order.
+         * then current implementation takes only the offset of the first bank to read.
+         * This field contains the value that we want in increasing order.
          *
          * @see #compareTo(Tile)
          */
@@ -204,16 +220,13 @@ class DataSubset extends TiledGridCoverage implements Localized {
          *
          * @param iterator  the iterator for which to create a snapshot of its current position.
          */
-        Tile(final AOI domain, final Vector tileOffsets, final int numTiles) {
+        Tile(final AOI domain, final Vector tileOffsets, final int[] selectedBanks, final int numTiles) {
             super(domain);
             int p = indexInTileVector;
-            long offset = tileOffsets.longValue(p);
-            final int limit = tileOffsets.size() - numTiles;
-            while (p < limit) {
-                // TODO: should take in account only the bands selected by user.
-                offset = Math.min(offset, tileOffsets.longValue(p += numTiles));
+            if (selectedBanks != null) {
+                p += selectedBanks[0] * numTiles;
             }
-            byteOffset = offset;
+            byteOffset = tileOffsets.longValue(p);
         }
 
         /**
@@ -225,11 +238,10 @@ class DataSubset extends TiledGridCoverage implements Localized {
          * @param  target    the array where to copy vector values.
          * @param  numTiles  value of {@link DataSubset#numTiles}.
          */
-        final void copyTileInfo(final Vector source, final long[] target, final int numTiles) {
-            int i = indexInTileVector;
+        final void copyTileInfo(final Vector source, final long[] target, final int[] selectedBanks, final int numTiles) {
             for (int j=0; j<target.length; j++) {
+                final int i = indexInTileVector + numTiles * (selectedBanks != null ? selectedBanks[j] : j);
                 target[j] = source.longValue(i);
-                i += numTiles;
             }
         }
 
@@ -263,10 +275,14 @@ class DataSubset extends TiledGridCoverage implements Localized {
         /*
          * Prepare an array for all tiles to be returned. Tiles that are already in memory will be stored
          * in this array directly. Other tiles will be declared in the `missings` array and loaded later.
+         * Each tile will either store all sample values in an interleaved fashion inside a single bank
+         * (`sourcePixelStride` > 1) or use one separated bank per band (`sourcePixelStride` == 1).
          */
+        final int[] selectedBanks = (sourcePixelStride == 1) ? selectedBands : null;
         final WritableRaster[] result = new WritableRaster[iterator.tileCountInQuery];
         final Tile[] missings = new Tile[iterator.tileCountInQuery];
         int numMissings = 0;
+        boolean needsCompaction = false;
         synchronized (reader().store) {
             do {
                 final WritableRaster tile = iterator.getCachedTile();
@@ -274,7 +290,7 @@ class DataSubset extends TiledGridCoverage implements Localized {
                     result[iterator.getIndexInResultArray()] = tile;
                 } else {
                     // Tile not yet loaded. Add to a queue of tiles to load later.
-                    missings[numMissings++] = new Tile(iterator, tileOffsets, numTiles);
+                    missings[numMissings++] = new Tile(iterator, tileOffsets, selectedBanks, numTiles);
                 }
             } while (iterator.next());
             Arrays.sort(missings, 0, numMissings);
@@ -292,14 +308,13 @@ class DataSubset extends TiledGridCoverage implements Localized {
             final Point  origin      = new Point();
             final long[] offsets     = new long[numBanks];
             final long[] byteCounts  = new long[numBanks];
-            boolean needsCompaction  = false;
             for (int i=0; i<numMissings; i++) {
                 final Tile tile = missings[i];
                 if (tile.getRegionInsideTile(lower, upper, subsampling, BIDIMENSIONAL)) {
                     origin.x = tile.originX;
                     origin.y = tile.originY;
-                    tile.copyTileInfo(tileOffsets,    offsets,    numTiles);
-                    tile.copyTileInfo(tileByteCounts, byteCounts, numTiles);
+                    tile.copyTileInfo(tileOffsets,    offsets,    selectedBanks, numTiles);
+                    tile.copyTileInfo(tileByteCounts, byteCounts, selectedBanks, numTiles);
                     for (int b=0; b<offsets.length; b++) {
                         offsets[b] = addExact(offsets[b], reader().origin);
                     }
@@ -309,18 +324,18 @@ class DataSubset extends TiledGridCoverage implements Localized {
                     needsCompaction = true;
                 }
             }
-            /*
-             * If the subsampling is larger than tile size, some tiles were empty and excluded.
-             * The corresponding elements in the `result` array were left to the null value.
-             * We need to compact the array by removing those null elements.
-             */
-            if (needsCompaction) {
-                int n = 0;
-                for (final WritableRaster tile : result) {
-                    if (tile != null) result[n++] = tile;
-                }
-                return Arrays.copyOf(result, n);
+        }
+        /*
+         * If the subsampling is larger than tile size, some tiles were empty and excluded.
+         * The corresponding elements in the `result` array were left to the null value.
+         * We need to compact the array by removing those null elements.
+         */
+        if (needsCompaction) {
+            int n = 0;
+            for (final WritableRaster tile : result) {
+                if (tile != null) result[n++] = tile;
             }
+            return Arrays.copyOf(result, n);
         }
         return result;
     }
@@ -342,7 +357,7 @@ class DataSubset extends TiledGridCoverage implements Localized {
      * @param  upper        (<var>x</var>, <var>y</var>) coordinates after the last pixel to read relative to the tile.
      * @param  subsampling  (<var>sx</var>, <var>sy</var>) subsampling factors.
      * @param  location     pixel coordinates in the upper-left corner of the tile to return.
-     * @return image decoded from the GeoTIFF file.
+     * @return a single tile decoded from the GeoTIFF file.
      * @throws IOException if an I/O error occurred.
      * @throws DataStoreException if a logical error occurred.
      * @throws RuntimeException if the Java2D image can not be created for another reason
@@ -360,16 +375,16 @@ class DataSubset extends TiledGridCoverage implements Localized {
          * is smaller than the "standard" tile size of the image. It happens often when reading the last strip.
          */
         final long length  = multiplyExact(type.size() / Byte.SIZE,
-                             multiplyExact(multiplyExact(width, height), numInterleaved));
-        final long[] size = new long[] {multiplyFull(numInterleaved, getTileSize(0)), getTileSize(1)};
+                             multiplyExact(multiplyExact(width, height), sourcePixelStride));
+        final long[] size = new long[] {multiplyFull(sourcePixelStride, getTileSize(0)), getTileSize(1)};
         /*
-         * If we use an interleaved sample model, each "element" from `HyperRectangleReader` perspective is actually
-         * a group of `numInterleaved` values. Note that in such case, we can not handle subsampling on the first axis.
+         * If we use an interleaved sample model, each "element" from `HyperRectangleReader` perspective is actually a
+         * group of `sourcePixelStride` values. Note that in such case, we can not handle subsampling on the first axis.
          * Such case should be handled by the `CompressedSubset` subclass instead, even if there is no compression.
          */
-        assert numInterleaved == 1 || subsampling[0] == 1;
-        lower[0] *= numInterleaved;
-        upper[0] *= numInterleaved;
+        assert sourcePixelStride == 1 || subsampling[0] == 1;
+        lower[0] *= sourcePixelStride;
+        upper[0] *= sourcePixelStride;
         /*
          * Read each plane ("banks" in Java2D terminology). Note that a single bank contain all bands
          * in the interleaved sample model case.
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/ImageFileDirectory.java b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/ImageFileDirectory.java
index cb6540b..778c3f5 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/ImageFileDirectory.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/ImageFileDirectory.java
@@ -221,6 +221,8 @@ final class ImageFileDirectory extends DataCube {
      * The {@code samplesPerPixel} value is usually 1 for bilevel, grayscale and palette-color images,
      * and 3 for RGB images. If this value is higher, then the {@code ExtraSamples} TIFF tag should
      * give an indication of the meaning of the additional channels.
+     *
+     * @see #getNumBands()
      */
     private short samplesPerPixel;
 
@@ -1360,6 +1362,14 @@ final class ImageFileDirectory extends DataCube {
     }
 
     /**
+     * Returns the number of components per pixel.
+     */
+    @Override
+    protected int getNumBands() {
+        return samplesPerPixel;
+    }
+
+    /**
      * Returns the size of tiles. This is also the size of the image sample model.
      * The number of dimensions is always 2 for {@code ImageFileDirectory}.
      *
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/AbstractGridResource.java b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/AbstractGridResource.java
index c2804a3..d502e8c 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/AbstractGridResource.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/AbstractGridResource.java
@@ -168,7 +168,7 @@ public abstract class AbstractGridResource extends AbstractResource implements G
                 previous = r;
             }
         }
-        return new RangeArgument(packed);
+        return new RangeArgument(packed, packed.length == numSampleDimensions);
     }
 
     /**
@@ -184,6 +184,11 @@ public abstract class AbstractGridResource extends AbstractResource implements G
         private final long[] packed;
 
         /**
+         * Whether the selection contains all bands of the resource, not necessarily in order.
+         */
+        public final boolean hasAllBands;
+
+        /**
          * If a {@linkplain #insertSubsampling subsampling} has been applied, indices of the first and last band
          * to read, together with the interval (stride) between bands.  Those information are computed only when
          * the {@code insertFoo(…)} methods are invoked.
@@ -201,9 +206,10 @@ public abstract class AbstractGridResource extends AbstractResource implements G
         /**
          * Encapsulates the given {@code range} argument packed in high bits.
          */
-        private RangeArgument(final long[] packed) {
-            this.packed = packed;
-            interval = 1;
+        private RangeArgument(final long[] packed, final boolean hasAllBands) {
+            this.packed      = packed;
+            this.hasAllBands = hasAllBands;
+            this.interval    = 1;
         }
 
         /**
@@ -213,7 +219,9 @@ public abstract class AbstractGridResource extends AbstractResource implements G
          * @return whether user specified all bands in increasing order without subsampling inserted.
          */
         public boolean isIdentity() {
-            if (interval != 1) return false;
+            if (!hasAllBands || interval != 1) {
+                return false;
+            }
             for (int i=0; i<packed.length; i++) {
                 if (packed[i] != ((((long) i) << Integer.SIZE) | i)) {
                     return false;
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridCoverage.java b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridCoverage.java
index f9abba0..296e6a9 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridCoverage.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridCoverage.java
@@ -119,6 +119,7 @@ public abstract class TiledGridCoverage extends GridCoverage {
      *
      * This transform maps {@linkplain org.opengis.referencing.datum.PixelInCell#CELL_CORNER pixel corners}.
      *
+     * @see #getSubsampling(int)
      * @see #toFullResolution(long, int)
      */
     private final int[] subsampling, subsamplingOffsets;
@@ -235,6 +236,16 @@ public abstract class TiledGridCoverage extends GridCoverage {
     }
 
     /**
+     * Returns the subsampling in the given dimension.
+     *
+     * @param  dimension  dimension for which to get subsampling.
+     * @return subsampling as a value ≥ 1.
+     */
+    protected final int getSubsampling(final int dimension) {
+        return subsampling[dimension];
+    }
+
+    /**
      * Converts a cell coordinate from this {@code TiledGridCoverage} coordinate space to full resolution.
      * This method removes the subsampling effect. Note that since this {@code TiledGridCoverage} uses the
      * same coordinate space than {@link TiledGridResource}, the converted coordinates should be valid in
@@ -581,6 +592,10 @@ public abstract class TiledGridCoverage extends GridCoverage {
          * if the enclosing {@link TiledGridCoverage} contains only one (potentially big) tile.
          * In that case, the reading process is more like untiled image reading.
          *
+         * <p>The {@link TiledGridCoverage} subsampling is provided for convenience,
+         * but is constant for all tiles regardless the subregion to read.
+         * The same values can be obtained by {@link #getSubsampling(int)}.</p>
+         *
          * @param  lower        a pre-allocated array where to store relative coordinates of the first pixel.
          * @param  upper        a pre-allocated array where to store relative coordinates after the last pixel.
          * @param  subsampling  a pre-allocated array where to store subsampling.
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridResource.java b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridResource.java
index cdd4e59..73a662e 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridResource.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridResource.java
@@ -259,7 +259,8 @@ public abstract class TiledGridResource extends AbstractGridResource {
                     selectedBands[i] = rangeIndices.getSourceIndex(i);
                 }
                 assert ArraysExt.isSorted(selectedBands, true);
-                if (ArraysExt.isRange(0, selectedBands)) {
+                if (rangeIndices.hasAllBands) {
+                    assert ArraysExt.isRange(0, selectedBands);
                     selectedBands = null;
                 }
             }
@@ -285,5 +286,16 @@ public abstract class TiledGridResource extends AbstractGridResource {
         public boolean hasSubsampling(final int dimension) {
             return subsampling[dimension] != 1;
         }
+
+        /**
+         * Returns {@code true} if the user requested a subset of bands, or {@code false}
+         * if all bands shall be read in order. For this method, a change of band order
+         * is considered as a "subset".
+         *
+         * @return {@code true} if only a subset of bands will be read or if bands will be read out of order.
+         */
+        public boolean hasBandSubset() {
+            return selectedBands != null;
+        }
     }
 }
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/io/HyperRectangleReader.java b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/io/HyperRectangleReader.java
index f2a93d2..122336c 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/io/HyperRectangleReader.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/io/HyperRectangleReader.java
@@ -28,6 +28,17 @@ import org.apache.sis.storage.DataStoreContentException;
  * Helper methods for reading a rectangular area, a cube or a hyper-cube from a channel.
  * The data can be stored in an existing array, or a new array can be created.
  * This class does not handle compression; it is rather designed for efficient reading of uncompressed data.
+ * It tries to read the largest possible contiguous blocks of data with single
+ * {@link java.nio.channels.ReadableByteChannel#read(ByteBuffer)} and
+ * {@link ByteBuffer#get(byte[], int, int)} method calls.
+ *
+ * <p>This reader supports subsampling in any dimension. However subsampling in the first dimension
+ * (the one with fastest varying index) is generally not efficient because it forces a large amount
+ * of seek operations. This class makes no special case for making that specific subsampling faster.
+ * It is generally not worth because subsampling in the first dimension is a special case anyway.
+ * It is the "dimension" of bands in an image using the pixel interleaved sample model, so the caller
+ * often needs to process subsampling in the first dimension in a different way than other dimensions
+ * anyway.</p>
  *
  * @author  Johann Sorel (Geomatys)
  * @author  Martin Desruisseaux (Geomatys)
diff --git a/storage/sis-storage/src/test/java/org/apache/sis/test/storage/CoverageReadConsistency.java b/storage/sis-storage/src/test/java/org/apache/sis/test/storage/CoverageReadConsistency.java
index a3f5286..b9665e2 100644
--- a/storage/sis-storage/src/test/java/org/apache/sis/test/storage/CoverageReadConsistency.java
+++ b/storage/sis-storage/src/test/java/org/apache/sis/test/storage/CoverageReadConsistency.java
@@ -28,8 +28,10 @@ import org.apache.sis.coverage.grid.GridExtent;
 import org.apache.sis.storage.DataStoreException;
 import org.apache.sis.storage.GridCoverageResource;
 import org.apache.sis.internal.util.StandardDateFormat;
+import org.apache.sis.internal.util.Numerics;
 import org.apache.sis.image.PixelIterator;
 import org.apache.sis.math.Statistics;
+import org.apache.sis.util.ArraysExt;
 import org.apache.sis.test.DependsOnMethod;
 import org.apache.sis.test.TestUtilities;
 import org.apache.sis.test.TestCase;
@@ -73,6 +75,23 @@ public strictfp class CoverageReadConsistency extends TestCase {
     private final GridCoverage full;
 
     /**
+     * Whether to allow random sub-regions to start elsewhere than (0,0).
+     */
+    private boolean allowOffsets;
+
+    /**
+     * Whether to use random subsamplings.
+     */
+    private boolean allowSubsamplings;
+
+    /**
+     * Whether to use random selection of bands.
+     *
+     * @see #randomRange()
+     */
+    private boolean allowBandSubset;
+
+    /**
      * The random number generator to use.
      */
     private final Random random;
@@ -133,7 +152,7 @@ public strictfp class CoverageReadConsistency extends TestCase {
      */
     @Test
     public void testSubRegionAtOrigin() throws DataStoreException {
-        readAndCompareRandomRegions(false, false);
+        readAndCompareRandomRegions();
     }
 
     /**
@@ -145,7 +164,8 @@ public strictfp class CoverageReadConsistency extends TestCase {
     @Test
     @DependsOnMethod("testSubRegionAtOrigin")
     public void testSubRegionsAnywhere() throws DataStoreException {
-        readAndCompareRandomRegions(true, false);
+        allowOffsets = true;
+        readAndCompareRandomRegions();
     }
 
     /**
@@ -157,7 +177,8 @@ public strictfp class CoverageReadConsistency extends TestCase {
     @Test
     @DependsOnMethod("testSubRegionAtOrigin")
     public void testSubsamplingAtOrigin() throws DataStoreException {
-        readAndCompareRandomRegions(false, true);
+        allowSubsamplings = true;
+        readAndCompareRandomRegions();
     }
 
     /**
@@ -169,26 +190,114 @@ public strictfp class CoverageReadConsistency extends TestCase {
     @Test
     @DependsOnMethod({"testSubsamplingAtOrigin", "testSubRegionsAnywhere"})
     public void testSubsamplingAnywhere() throws DataStoreException {
-        readAndCompareRandomRegions(true, true);
+        allowOffsets      = true;
+        allowSubsamplings = true;
+        readAndCompareRandomRegions();
+    }
+
+    /**
+     * Tests reading random subset of bands in random sub-region starting at coordinates (0,0).
+     *
+     * @throws DataStoreException if an error occurred while using the resource.
+     */
+    @Test
+    @DependsOnMethod("testSubRegionAtOrigin")
+    public void testBandSubsetAtOrigin() throws DataStoreException {
+        allowBandSubset = true;
+        readAndCompareRandomRegions();
+    }
+
+    /**
+     * Tests reading random subset of bands in random sub-regions starting at random offsets.
+     *
+     * @throws DataStoreException if an error occurred while using the resource.
+     */
+    @Test
+    @DependsOnMethod({"testBandSubsetAtOrigin", "testSubRegionsAnywhere"})
+    public void testBandSubsetAnywhere() throws DataStoreException {
+        allowOffsets    = true;
+        allowBandSubset = true;
+        readAndCompareRandomRegions();
+    }
+
+    /**
+     * Tests reading random subset of bands in random sub-regions starting at random offsets
+     * with random subsampling applied.
+     *
+     * @throws DataStoreException if an error occurred while using the resource.
+     */
+    @Test
+    @DependsOnMethod({"testBandSubsetAnywhere", "testSubsamplingAnywhere"})
+    public void testAllAnywhere() throws DataStoreException {
+        allowOffsets      = true;
+        allowBandSubset   = true;
+        allowSubsamplings = true;
+        readAndCompareRandomRegions();
+    }
+
+    /**
+     * Creates a random domain to be used as a query on the {@link #resource} to test.
+     * All arrays given to this method will have their values overwritten.
+     *
+     * @param  gg           value of {@link GridCoverage#getGridGeometry()} on the resource to test.
+     * @param  low          pre-allocated array where to write the lower grid coordinates, inclusive.
+     * @param  high         pre-allocated array where to write the upper grid coordinates, inclusive.
+     * @param  subsampling  pre-allocated array where to write the subsampling.
+     */
+    private GridGeometry randomDomain(final GridGeometry gg, final long[] low, final long[] high, final int[] subsampling) {
+        final GridExtent fullExtent = gg.getExtent();
+        final int dimension = fullExtent.getDimension();
+        for (int d=0; d<dimension; d++) {
+            final int span = StrictMath.toIntExact(fullExtent.getSize(d));
+            final int rs = random.nextInt(span);                            // Span of the sub-region - 1.
+            if (allowOffsets) {
+                low[d] = random.nextInt(span - rs);                         // Note: (span - rs) > 0.
+            }
+            high[d] = low[d] + rs;
+            subsampling[d] = 1;
+            if (allowSubsamplings) {
+                subsampling[d] += random.nextInt(StrictMath.max(rs / 16, 1));
+            }
+        }
+        return gg.derive().subgrid(new GridExtent(null, low, high, true), subsampling).build();
+    }
+
+    /**
+     * Returns the subset of bands to use for testing. This method never return {@code null},
+     * but the set of bands is random only if {@link #allowBandSubset} is {@code true}.
+     */
+    private int[] randomRange(final int numBands) {
+        if (!allowBandSubset) {
+            return ArraysExt.range(0, numBands);
+        }
+        final int[] selectedBands = new int[numBands];
+        for (int i=0; i<numBands; i++) {
+            selectedBands[i] = random.nextInt(numBands);
+        }
+        // Remove duplicated elements.
+        long included = 0;
+        int count = 0;
+        for (final int b : selectedBands) {
+            if (included != (included |= Numerics.bitmask(b))) {
+                selectedBands[count++] = b;
+            }
+        }
+        return ArraysExt.resize(selectedBands, count);
     }
 
     /**
      * Implementation of methods testing reading in random sub-regions with random sub-samplings.
      *
-     * @param  allowOffsets       whether to allow sub-regions to start elsewhere than (0,0).
-     * @param  allowSubsamplings  whether to use random subsamplings.
      * @throws DataStoreException if an error occurred while using the resource.
      */
-    private void readAndCompareRandomRegions(final boolean allowOffsets, final boolean allowSubsamplings)
-            throws DataStoreException
-    {
+    private void readAndCompareRandomRegions() throws DataStoreException {
         final GridGeometry gg = resource.getGridGeometry();
-        final GridExtent fullExtent = gg.getExtent();
-        final int    dimension   = fullExtent.getDimension();
+        final int    dimension   = gg.getDimension();
         final long[] low         = new long[dimension];
         final long[] high        = new long[dimension];
         final int [] subsampling = new int [dimension];
         final int [] subOffsets  = new int [dimension];
+        final int    numBands    = resource.getSampleDimensions().size();
         /*
          * We will collect statistics on execution time only if the
          * test is executed in a more verbose mode than the default.
@@ -196,30 +305,15 @@ public strictfp class CoverageReadConsistency extends TestCase {
         final Statistics durations = (VERBOSE || !failOnMismatch) ? new Statistics("time (ms)") : null;
         int failuresCount = 0;
         for (int it=0; it < numIterations; it++) {
-            /*
-             * Create a random domain to be used as a query on the `GridCoverageResource`.
-             */
-            for (int d=0; d<dimension; d++) {
-                final int span = StrictMath.toIntExact(fullExtent.getSize(d));
-                final int rs = random.nextInt(span);                    // Span of the sub-region - 1.
-                if (allowOffsets) {
-                    low[d] = random.nextInt(span - rs);                 // Note: (span - rs) > 0.
-                }
-                high[d] = low[d] + rs;
-                subsampling[d] = 1;
-                if (allowSubsamplings) {
-                    subsampling[d] += random.nextInt(StrictMath.max(rs / 16, 1));
-                }
-            }
-            final GridGeometry domain = gg.derive()
-                    .subgrid(new GridExtent(null, low, high, true), subsampling).build();
+            final GridGeometry domain = randomDomain(gg, low, high, subsampling);
+            final int[] selectedBands = randomRange(numBands);
             /*
              * Read a coverage containing the requested sub-domain. Note that the reader is free to read
              * more data than requested. The extent actually read is `actualReadExtent`. It shall contain
              * fully the requested `domain`.
              */
             final long startTime = System.nanoTime();
-            final GridCoverage subset = resource.read(domain, null);
+            final GridCoverage subset = resource.read(domain, selectedBands);
             final GridExtent actualReadExtent = subset.getGridGeometry().getExtent();
             if (failOnMismatch) {
                 assertEquals("Unexpected number of dimensions.", dimension, actualReadExtent.getDimension());
@@ -251,16 +345,17 @@ nextSlice:  for (;;) {
                 System.arraycopy(sliceMin, BIDIMENSIONAL, sliceMax, BIDIMENSIONAL, dimension - BIDIMENSIONAL);
                 final PixelIterator itr = iterator(full,   sliceMin, sliceMax, subsampling, subOffsets, allowSubsamplings);
                 final PixelIterator itc = iterator(subset, sliceMin, sliceMax, subsampling, subOffsets, false);
-                if (durations != null) {
-                    durations.accept((System.nanoTime() - startTime) / (double) StandardDateFormat.NANOS_PER_MILLISECOND);
-                }
                 if (itr != null) {
                     assertEquals(itr.getDomain().getSize(), itc.getDomain().getSize());
-                    double[] expected = null, actual = null;
+                    final double[] expected = new double[selectedBands.length];
+                    double[] reference = null, actual = null;
                     while (itr.next()) {
                         assertTrue(itc.next());
-                        expected = itr.getPixel(expected);
-                        actual   = itc.getPixel(actual);
+                        reference = itr.getPixel(reference);
+                        actual    = itc.getPixel(actual);
+                        for (int i=0; i<selectedBands.length; i++) {
+                            expected[i] = reference[selectedBands[i]];
+                        }
                         if (!Arrays.equals(expected, actual)) {
                             failuresCount++;
                             if (!failOnMismatch) break;
@@ -269,7 +364,7 @@ nextSlice:  for (;;) {
                             final StringBuilder message = new StringBuilder(100).append("Mismatch at position (")
                                     .append(pr.x).append(", ").append(pr.y).append(") in full image and (")
                                     .append(pc.x).append(", ").append(pc.y).append(") in tested sub-image");
-                            findMatchPosition(itr, pr, expected, message);
+                            findMatchPosition(itr, pr, selectedBands, actual, message);
                             assertArrayEquals(message.toString(), expected, actual, STRICT);
                         }
                     }
@@ -291,6 +386,9 @@ nextSlice:  for (;;) {
                 }
                 break;
             }
+            if (durations != null) {
+                durations.accept((System.nanoTime() - startTime) / (double) StandardDateFormat.NANOS_PER_MILLISECOND);
+            }
         }
         /*
          * Show statistics only if the test are executed with the `VERBOSE` flag set,
@@ -373,8 +471,11 @@ nextSlice:  for (;;) {
      * Explores pixel values around the given position in search for a pixel having the expected values.
      * If a match is found, the error message is completed with information about the match position.
      */
-    private static void findMatchPosition(final PixelIterator ir, final Point pr, final double[] expected, final StringBuilder message) {
-        double[] actual = null;
+    private static void findMatchPosition(final PixelIterator ir, final Point pr, final int[] selectedBands,
+                                          final double[] actual, final StringBuilder message)
+    {
+        final double[] expected = new double[actual.length];
+        double[] reference = null;
         for (int dy=0; dy<10; dy++) {
             for (int dx=0; dx<10; dx++) {
                 if ((dx | dy) != 0) {
@@ -386,7 +487,10 @@ nextSlice:  for (;;) {
                         } catch (IndexOutOfBoundsException e) {
                             continue;
                         }
-                        actual = ir.getPixel(actual);
+                        reference = ir.getPixel(reference);
+                        for (int i=0; i<selectedBands.length; i++) {
+                            expected[i] = reference[selectedBands[i]];
+                        }
                         if (Arrays.equals(expected, actual)) {
                             message.append(" (note: found a match at offset (").append(x).append(", ").append(y)
                                    .append(") in full image)");

Mime
View raw message