--- a/scripts/image/private/__tiff_imread__.m	Sat Sep 03 20:36:33 2022 +0200
+++ b/scripts/image/private/__tiff_imread__.m	Sat Sep 03 23:51:06 2022 +0200
@@ -122,11 +122,23 @@
   alpha = [];
   for page_idx = 1:numel(pages)
     img.setDirectory (pages(page_idx));
-    ## FIXME: This an ineffecient way to read pixel regions because it
-    ## always reads the entire image first. A better way is to figure
-    ## out which strips/tiles are needed for the region and read only those.
-    data = img.read ();
-    data = data (region{1}, region{2}, :);
+
+    tiled = img.isTiled ();
+    if (numel (region{1}) == info.height
+        && (! tiled || numel (region{2}) == info.width))
+      ## If we need the entire tiled image or all the rows of a stripped
+      ## image then the best approach is to read the entire image directly
+      data = img.read ();
+      if (numel (region{2}) != info.width)
+        data = data (:, region{2}, :);
+      endif
+    elseif (tiled)
+      ## Otherwise it should be better to read the specific strips and tiles
+      ## that the requested region falls in
+      data = read_tiled_data (img, region);
+    else
+      data = read_stripped_data (img, region);
+    endif
     A = cat (4, A, data);
     if (info.photometric == Tiff.Photometric.Palette)
       cmap = cat (3, cmap, img.getTag (Tiff.TagID.ColorMap));
@@ -136,6 +148,64 @@
   img.close ();
 endfunction
 
+function data = read_tiled_data (img, region)
+  info = get_image_info (img);
+  tile = img.readEncodedTile (1);
+  data = zeros (numel (region{1}), numel (region{2}),
+                info.nchannels, class (tile));
+
+endfunction
+
+function data = read_stripped_data (img, region)
+  ## Get all necessary image info
+  info = get_image_info (img);
+  rows_per_strip = img.getTag (Tiff.TagID.RowsPerStrip);
+  planar = img.getTag (Tiff.TagID.PlanarConfiguration);
+  
+  ## Create array to hold the output data
+  ## (we need a strip to infer the data type)
+  strip = img.readEncodedStrip (1);
+  data = zeros (numel (region{1}), numel (region{2}),
+                info.nchannels, class (strip));
+  last_strip_index = 1;
+  if (planar == Tiff.PlanarConfiguration.Chunky)
+    ## For chunky planes, we need to iterate over rows only
+    data_row = 1;
+    for row = region{1}
+      strip_index = img.computeStrip (row);
+      ## The index of the row in the strip
+      strip_row = rem (row - 1, rows_per_strip) + 1;
+      if (strip_index != last_strip_index)
+        ## If the required strip isn't cached, grab it
+        strip = img.readEncodedStrip (strip_index);
+        last_strip_index = strip_index;
+      endif
+      data(data_row, :, :) = strip (strip_row, region{2}, :);
+      data_row += 1;
+    endfor
+  else
+    ## For separate planes we need to loop over channels as well as rows
+    ## It's better to make the outer loop for channels so we can make use
+    ## of strip caching, because otherwise, we would need to read a different
+    ## strip each iteration of the inner loop
+    for ch = 1:info.nchannels
+      data_row = 1;
+      for row = region{1}
+        strip_index = img.computeStrip (row, ch);
+        ## The index of the row in the strip
+        strip_row = rem (row - 1, rows_per_strip) + 1;
+        if (strip_index != last_strip_index)
+          ## If the needed strip is not cached, grab it
+          strip = img.readEncodedStrip (strip_index);
+          last_strip_index = strip_index;
+        endif
+        data(data_row, :, ch) = strip (strip_row, region{2});
+        data_row += 1;
+      endfor
+    endfor
+  endif
+endfunction
+
 function bool = is_valid_index_option (arg)
   bool = isvector (arg) && isnumeric (arg) && isreal (arg);
 endfunction