--- a/libinterp/dldfcn/__tiff__.cc	Tue Jul 19 19:45:06 2022 +0200
+++ b/libinterp/dldfcn/__tiff__.cc	Wed Jul 20 02:48:22 2022 +0200
@@ -51,8 +51,6 @@
   {
     typedef typename T::element_type P;
 
-    T img;
-
     // For 1-bit and 4-bit images, each row must be byte aligned and padding
     // is added to the end of the row to reach the byte mark.
     // To facilitate reading the data, the matrix is defined with the padded
@@ -69,6 +67,7 @@
     // The matrix dimensions are defined in the order that corresponds to
     // the order of strip data read from LibTIFF.
     // At the end, the matrix is permutated to the order expected by Octave
+    T img;
     if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
       img = T (dim_vector (image_data->samples_per_pixel, padded_width,
                            image_data->height));
@@ -193,46 +192,168 @@
   {
     typedef typename T::element_type P;
 
-    T img = T (dim_vector (image_data->height, image_data->width,
-                            image_data->samples_per_pixel));
-    P *img_fvec = img.fortran_vec ();
-
     // Obtain the necessary data for handling the tiles
     uint32_t tile_width, tile_height;
     validate_tiff_get_field (TIFFGetField (tif, TIFFTAG_TILEWIDTH,
                                            &tile_width));
     validate_tiff_get_field (TIFFGetField (tif, TIFFTAG_TILELENGTH,
                                            &tile_height));
+    uint32_t tile_count = TIFFNumberOfTiles (tif);
+    uint32_t tiles_across = (image_data->width + tile_width - 1)
+                            / tile_width;
+    uint32_t tiles_down = (image_data->height + tile_height - 1)
+                          / tile_height;
 
-    tdata_t buf = _TIFFmalloc (TIFFTileSize (tif));
-    if (! buf)
-      error ("Failed to allocate buffer for tile data");
-    
-    // TODO(maged): implement support for ImageDepth?
-    uint32_t tile_size;
-    for (uint32_t base_row = 0; base_row < image_data->height; base_row += tile_height)
-      for (uint32_t base_col = 0; base_col < image_data->width; base_col += tile_width)
-        if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+    T img;
+    // The matrix dimensions are defined in the order that corresponds to
+    // the order of the tile data read from LibTIFF.
+    // At the end, the matrix is permutated, reshaped and resized to be in the
+    // shape expected by Octave
+    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+      img = T (dim_vector (image_data->samples_per_pixel, tile_width,
+                           tile_height, tiles_across, tiles_down));
+    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      img = T (dim_vector (tile_width, tile_height, tiles_across, tiles_down,
+                           image_data->samples_per_pixel));
+    else
+      error ("Unsupported Planar Configuration");
+
+    P *img_fvec = img.fortran_vec ();
+
+    // image_size is calculated from the tile data and not from the
+    // image parameters to account for the additional padding in the
+    // boundary tiles
+    uint64_t image_size = tile_width * tile_height * tile_count * sizeof(P);
+    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+      image_size *= image_data->samples_per_pixel;
+
+    // Can't rely on TileByteCounts because compressed images will report
+    // the number of bytes present in the file as opposed to the actual
+    // number of bytes of uncompressed data that is needed here
+    int64_t tile_size;    
+    uint32_t written_size = 0;
+    for (uint32_t tile = 0; tile < tile_count; tile++)
+      {
+        tile_size = TIFFReadEncodedTile(tif, tile, img_fvec, -1);
+        
+        if (tile_size == -1)
+          error ("Failed to read tile data");
+
+        // Checking if the read bytes would exceed the size of the matrix
+        if (tile_size + written_size > image_size)
+          error ("Tile data is larger than image size");
+        
+        if (image_data->bits_per_sample == 1)
           {
-            if ((tile_size = TIFFReadTile(tif, buf, base_col, base_row, 0, 0)) == -1)
-              error ("Failed to read tile data");
+            if (image_data->samples_per_pixel != 1)
+              error ("Bi-Level images must have one channel only");
+            
+            // The tile size is multiplied by 8 to reflect tha actual
+            // number of bytes written to the matrix since each byte
+            // in the original tile contains 8 pixels of data
+            tile_size *= 8;
+
+            // Checking bounds again with the new size
+            if (written_size + tile_size > image_size)
+              error ("Tile data is larger than the image size");
+
+            // Iterate over the memory region backwards to expand the bits
+            // to their respective bytes without overwriting the read data
+            for (int64_t pixel = tile_size - 1; pixel >= 0; pixel--)
+              {
+                // TODO(maged): is it necessary to check FillOrder?
+                uint8_t bit_number = 7 - pixel % 8;
+                img_fvec[pixel] = (img_fvec[pixel / 8] >> bit_number) & 0x01;
+              }
+            break;
+          }
+        else if (image_data->bits_per_sample == 4)
+          {
+            if (image_data->samples_per_pixel != 1)
+              error ("4-bit images are only supported for grayscale");
             
-            // Boundary tiles are zero padded so we might need to stop earlier than the end of the tile
-            for (uint32_t sub_row = 0; sub_row < tile_height && base_row + sub_row < image_data->height; sub_row++)
-              for (uint32_t sub_col = 0; sub_col < tile_width && base_col + sub_col < image_data->width; sub_col++)
-                for (uint16_t sample = 0; sample < image_data->samples_per_pixel; sample++)
-                  img_fvec[sample * image_data->width * image_data->height
-                           + (base_col + sub_col) * image_data->height
-                           + base_row + sub_row]
-                    = ((P *)buf)[sub_row * tile_width * image_data->samples_per_pixel
-                                 + sub_col * image_data->samples_per_pixel
-                                 + sample];
+            // tile size is multplied by as each byte contains 2 pixels
+            // and each pixel is later expanded into its own byte
+            tile_size *= 2;
+
+            // Checking bounds again with the ne size
+            if (written_size + tile_size > image_size)
+              error ("Tile data is larger than the image size");
+            
+            // Iterate over the memory region backwards to expand the nibbles
+            // to their respective bytes without overwriting the read data
+            for (int64_t pixel = tile_size - 1; pixel >= 0; pixel--)
+              {
+                uint8_t shift = pixel % 2 == 0? 4: 0;
+                img_fvec[pixel] = (img_fvec[pixel / 2] >> shift) & 0x0F;
+              }
+            break;
           }
-        else
-          // TODO(maged): Implement for separated planes
-          error ("Tiled images with separate planar configuration not supported");
+        else if (image_data->bits_per_sample != 8 &&
+                 image_data->bits_per_sample != 16 &&
+                 image_data->bits_per_sample != 32 &&
+                 image_data->bits_per_sample != 64)
+          error ("Unsupported bit depth");
+        
+        img_fvec = (P*)((uint8_t *)img_fvec + tile_size);
+        written_size += tile_size;
+      }
+    
+    // The data is now in the matrix but in a different order than expected
+    // by Octave and with additional padding in boundary tiles.
+    // To get it to the right order, the dimensions are permutated to
+    // align tiles to their correct grid, then reshaped to remove the
+    // extra dimensions (tiles_across, tiles_down), then resized to
+    // remove any extra padding, and finally permutated to the correct
+    // order that is: height x width x channels
+    Array<octave_idx_type> perm1 (dim_vector (5, 1));
+    Array<octave_idx_type> perm2 (dim_vector (3, 1));
+    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+      {
+        perm1(0) = 0;
+        perm1(1) = 1;
+        perm1(2) = 3;
+        perm1(3) = 2;
+        perm1(4) = 4;
+        img = img.permute (perm1);
 
-    _TIFFfree (buf);
+        img = img.reshape (dim_vector (image_data->samples_per_pixel,
+                                       tile_width * tiles_across,
+                                       tile_height * tiles_down));
+        
+        if (tile_width * tiles_across != image_data->width
+            || tile_height * tiles_down != image_data->height)
+          img.resize (dim_vector (image_data->samples_per_pixel,
+                                  image_data->width, image_data->height));
+        
+        perm2(0) = 2;
+        perm2(1) = 1;
+        perm2(2) = 0;
+        img = img.permute (perm2);
+      }
+    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      {
+        perm1(0) = 0;
+        perm1(1) = 2;
+        perm1(2) = 1;
+        perm1(3) = 3;
+        perm1(4) = 4;
+        img = img.permute (perm1);
+
+        img = img.reshape (dim_vector (tile_width * tiles_across,
+                                       tile_height * tiles_down,
+                                       image_data->samples_per_pixel));
+        
+        if (tile_width * tiles_across != image_data->width
+            || tile_height * tiles_down != image_data->height)
+          img.resize (dim_vector (image_data->width, image_data->height,
+                                  image_data->samples_per_pixel));
+        
+        perm2(0) = 1;
+        perm2(1) = 0;
+        perm2(2) = 2;
+        img = img.permute (perm2);
+      }
 
     return octave_value (img);
   }