--- a/scripts/image/private/__tiff_imread__.m	Sat Sep 03 23:51:06 2022 +0200
+++ b/scripts/image/private/__tiff_imread__.m	Sun Sep 04 15:30:27 2022 +0200
@@ -149,11 +149,100 @@
 endfunction
 
 function data = read_tiled_data (img, region)
+  ## Get all necessary image info
   info = get_image_info (img);
+  planar = img.getTag (Tiff.TagID.PlanarConfiguration);
+  tile_width = img.getTag (Tiff.TagID.TileWidth);
+  tile_height = img.getTag (Tiff.TagID.TileLength);
+
+  ## Calculate the number of tiles in the two dimensions
+  tiles_across = (info.width + tile_width - 1) / tile_width;
+  tiles_down = (info.height + tile_height - 1) / tile_height;
+
+  ## Create array to hold output data (we need a tile to infer the datatype)
   tile = img.readEncodedTile (1);
   data = zeros (numel (region{1}), numel (region{2}),
                 info.nchannels, class (tile));
+  
+  ## Calculate the row stride
+  rows = region{1};
+  row_stride = uint32 (info.height);
+  if (numel (rows) > 1)
+    row_stride = uint32 (rows(2) - rows(1));
+  endif
 
+  ## Calculate the column stride
+  cols = region{2};
+  col_stride = uint32 (info.width);
+  if (numel (cols) > 1)
+    col_stride = uint32 (cols(2) - cols(1));
+  endif
+
+  tile_idx = 1;
+  if (planar == Tiff.PlanarConfiguration.Chunky)
+    ## For chunky planes all channels are stores in the same plane
+    for tj = 1:tiles_down
+      ## Calculate the index of the first row of the tile
+      t_rowstart = (tj - 1) * tile_height + 1;
+      ## Extract the rows that fall in the tile
+      rows_in_tile = rows (rows >= t_rowstart
+                          & rows < (t_rowstart + tile_height));
+      ## Calculate the indices of these rown in the tile array
+      tile_rows = rem (rows_in_tile - 1, tile_height) + 1;
+      ## Calculate the indices of the same rows in the output array
+      data_rows = idivide (rows_in_tile - rows (1), row_stride, "fix") + 1;
+      for ti = 1:tiles_across
+        ## Calculate the index of the first column of the tile
+        t_colstart = (ti - 1) * tile_width + 1;
+        ## Extract the columns that fall in the tile
+        cols_in_tile = cols (cols >= t_colstart
+                            & cols < (t_colstart + tile_width));
+        ## Calculate the indices of these columns in the tile array
+        tile_cols = rem (cols_in_tile - 1, tile_width) + 1;
+        ## Calculate the indices of the same columns in the output array
+        data_cols = idivide (cols_in_tile - cols (1), col_stride, "fix") + 1;
+
+        tile = img.readEncodedTile (tile_idx);
+        tile_idx += 1;
+
+        ## Copy the needed data from the tile array to the output array
+        data (data_rows, data_cols, :) = tile (tile_rows, tile_cols, :);
+      endfor
+    endfor
+  else
+    for ch = 1:info.nchannels
+      ## For separate planes all the tiles of the first channel are found
+      ## first then all the tiles of the second channel and so on
+      for tj = 1:tiles_down
+        ## Calculate the index of the first row of the tile
+        t_rowstart = (tj - 1) * tile_height + 1;
+        ## Extract the rows that fall in the tile
+        rows_in_tile = rows (rows >= t_rowstart
+                            & rows < (t_rowstart + tile_height));
+        ## Calculate the indices of these rown in the tile array
+        tile_rows = rem (rows_in_tile - 1, tile_height) + 1;
+        ## Calculate the indices of the same rows in the output array
+        data_rows = idivide (rows_in_tile - rows (1), row_stride, "fix") + 1;
+        for ti = 1:tiles_across
+          ## Calculate the index of the first column of the tile
+          t_colstart = (ti - 1) * tile_width + 1;
+          ## Extract the columns that fall in the tile
+          cols_in_tile = cols (cols >= t_colstart
+                              & cols < (t_colstart + tile_width));
+          ## Calculate the indices of these columns in the tile array
+          tile_cols = rem (cols_in_tile - 1, tile_width) + 1;
+          ## Calculate the indices of the same columns in the output array
+          data_cols = idivide (cols_in_tile - cols (1), col_stride, "fix") + 1;
+
+          tile = img.readEncodedTile (tile_idx);
+          tile_idx += 1;
+
+          ## Copy the needed data from the tile array to the output array
+          data (data_rows, data_cols, ch) = tile (tile_rows, tile_cols);
+        endfor
+      endfor
+    endfor
+  endif
 endfunction
 
 function data = read_stripped_data (img, region)