--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/__tiff__.cc	Wed Aug 10 00:58:12 2022 +0200
@@ -0,0 +1,2369 @@
+#if defined (HAVE_CONFIG_H)
+#  include "config.h"
+#endif
+
+#include <string>
+#include <iostream>
+
+// Workaround to not have to include fcntl which doesn't exist on windows
+// TODO(maged): see what octave does for this (fcntl)
+#ifndef O_RDONLY
+#define O_RDONLY 0
+#endif
+
+#include "defun.h"
+#include "ov.h"
+#include "ovl.h"
+#include "error.h"
+
+#include "errwarn.h"
+
+#if defined (HAVE_TIFF)
+#  include <tiffio.h>
+#endif
+
+namespace octave
+{
+#if defined (HAVE_TIFF)
+
+  struct tiff_image_data
+  {
+  public:
+    uint32_t width;
+    uint32_t height;
+    uint16_t samples_per_pixel;
+    uint16_t bits_per_sample;
+    uint16_t planar_configuration;
+    uint16_t is_tiled;
+
+    tiff_image_data (TIFF *tif)
+    {
+      if (! TIFFGetField (tif, TIFFTAG_IMAGEWIDTH, &width))
+        error ("Failed to read image width");
+
+      if (! TIFFGetField (tif, TIFFTAG_IMAGELENGTH, &height))
+        error ("Failed to read image height");
+      
+      if (! TIFFGetFieldDefaulted (tif, TIFFTAG_SAMPLESPERPIXEL,
+                                   &samples_per_pixel))
+        error ("Failed to read the SamplesPerPixel tag");
+
+      if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
+                                   &bits_per_sample))
+        error ("Failed to read the BitsPerSample tag");
+      
+      // TODO(maged): this doesn't really work for writing as LibTIFF will
+      // refuse to write unless the tag is set
+      if (! TIFFGetField (tif, TIFFTAG_PLANARCONFIG,
+                          &planar_configuration))
+        // LibTIFF has a bug where it incorrectly returns 0 as a default
+        // value for PlanarConfiguration although the default value is 1
+        planar_configuration = 1;
+      
+      is_tiled = TIFFIsTiled(tif);
+    }
+  };
+
+  // Error if status is not 1 (success status for TIFFGetField)
+  void
+  validate_tiff_get_field (bool status, void *p_to_free=NULL)
+  {
+      if (status != 1)
+        {
+          if (p_to_free != NULL)
+            _TIFFfree (p_to_free);
+          error ("Failed to read tag");
+        }
+  }
+
+  uint32_t get_rows_in_strip (uint32_t strip_no, uint32_t strip_count,
+                              uint32_t rows_per_strip,
+                              tiff_image_data *image_data)
+  {
+    // Calculate the expected number of elements in the strip data array
+    // All strips have equal number of rows except strips at the bottom
+    // of the image can have less number of rows
+    uint32_t rows_in_strip = rows_per_strip;
+    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+      {
+        // All strips have equal number of rows excpet strips at the bottom
+        // of the image can have less number of rows
+        if (strip_no == strip_count - 1)
+          rows_in_strip = image_data->height - rows_in_strip * strip_no;
+      }
+    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      {
+        // For separate planes, we should check the last strip of each channel
+        uint32_t strips_per_channel
+          = strip_count / image_data->samples_per_pixel;
+        for (uint32_t boundary_strip = strips_per_channel - 1;
+             boundary_strip <= strip_count;
+             boundary_strip += strips_per_channel)
+          if (strip_no == boundary_strip)
+            rows_in_strip = image_data->height
+                            - rows_in_strip * (strip_no % strips_per_channel);
+      }
+    else
+      error ("Planar Configuration not supported");
+    
+    return rows_in_strip;
+  }
+
+  template <typename T>
+  octave_value
+  read_strip (TIFF *tif, uint32_t strip_no, tiff_image_data * image_data)
+  {
+    // ASSUMES stripped image and strip_no is a valid zero-based strip
+    // index for the tif image
+    
+    uint32_t rows_in_strip;
+    if (! TIFFGetFieldDefaulted (tif, TIFFTAG_ROWSPERSTRIP, &rows_in_strip))
+      error ("Failed to obtain a value for RowsPerStrip");
+    
+    if (rows_in_strip > image_data->height)
+      rows_in_strip = image_data->height;
+    
+    uint32_t strip_count = TIFFNumberOfStrips (tif);
+    rows_in_strip = get_rows_in_strip (strip_no, strip_count,
+                                       rows_in_strip, image_data);
+    dim_vector strip_dims;
+    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+      strip_dims = dim_vector (image_data->samples_per_pixel,
+                               image_data->width, rows_in_strip);
+    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      strip_dims = dim_vector (image_data->width, rows_in_strip, 1);
+    else
+      error ("Unsupported bit depth");
+    
+    T strip_data (strip_dims);
+    uint8_t *strip_fvec
+      = reinterpret_cast<uint8_t *> (strip_data.fortran_vec ());
+
+    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)
+      {
+        if (TIFFReadEncodedStrip (tif, strip_no, strip_fvec, -1) == -1)
+          error ("Failed to read strip data");
+      }
+    else if (image_data->bits_per_sample == 1)
+      {
+        if (image_data->samples_per_pixel != 1)
+          error ("Bi-Level images must have one channel only");
+        
+        // Create a buffer to hold the packed strip data
+        // Unique pointers are faster than vectors for constant size buffers
+        std::unique_ptr<uint8_t []> strip_ptr
+          = std::make_unique<uint8_t []> (TIFFStripSize (tif));
+        uint8_t *strip_buf = strip_ptr.get ();
+        if (TIFFReadEncodedStrip (tif, strip_no, strip_buf, -1) == -1)
+          error ("Failed to read strip data");
+        
+        // According to the format specification, the row should be byte
+        // aligned so the number of bytes is rounded up to the nearest byte
+        uint32_t padded_width = (image_data->width + 7) / 8;
+        // Packing the pixel data into bits
+        for (uint32_t row = 0; row < rows_in_strip; row++)
+          {
+            for (uint32_t col = 0; col < image_data->width; col++)
+            {
+              uint8_t shift = 7 - col % 8;
+              strip_fvec[col] = (strip_buf[col / 8] >> shift) & 0x1;
+            }
+            strip_fvec += image_data->width;
+            strip_buf += padded_width;
+          }
+      }
+    else
+      error ("Unsupported bit depth");
+    
+    Array<octave_idx_type> perm (dim_vector (3, 1));
+    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+      {
+        perm(0) = 2;
+        perm(1) = 1;
+        perm(2) = 0;
+      }
+    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      {
+        perm(0) = 1;
+        perm(1) = 0;
+        perm(2) = 2;
+      }
+    
+    strip_data = strip_data.permute (perm);
+    return octave_value (strip_data);
+  }
+
+  template <typename T>
+  octave_value
+  read_tile (TIFF *tif, uint32_t tile_no, tiff_image_data *image_data)
+  {
+    // ASSUMES tiled image and tile_no is a valid zero-based tile
+    // index for the tif image
+
+    // TODO(maged): refactor into a function?
+    uint32_t tile_width, tile_height;
+    if (! TIFFGetField (tif, TIFFTAG_TILELENGTH, &tile_height))
+      error ("Filed to read tile length");
+    
+    if (! TIFFGetField (tif, TIFFTAG_TILEWIDTH, &tile_width))
+      error ("Filed to read tile length");
+
+    if (tile_height == 0 || tile_height % 16 != 0
+        || tile_width == 0 || tile_width % 16 != 0)
+      error ("Tile dimesion tags are invalid");
+
+    dim_vector tile_dims;
+    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+      {
+        tile_dims = dim_vector (image_data->samples_per_pixel, tile_width,
+                                tile_height);
+      }
+    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      {
+        tile_dims = dim_vector (tile_width, tile_height, 1);
+      }
+    else
+      error ("Unsupported planar configuration");
+    
+    T tile_data (tile_dims);
+    uint8_t *tile_fvec
+      = reinterpret_cast<uint8_t *> (tile_data.fortran_vec ());
+
+    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)
+        {      
+          if (TIFFReadEncodedTile (tif, tile_no, tile_fvec, -1) == -1)
+            error ("Failed to read tile data");
+        }
+    else if (image_data->bits_per_sample == 1)
+      {
+        if (image_data->samples_per_pixel != 1)
+          error ("Bi-Level images must have one channel only");
+        
+        // Create a buffer to hold the packed tile data
+        // Unique pointers are faster than vectors for constant size buffers
+        std::unique_ptr<uint8_t []> tile_ptr
+          = std::make_unique<uint8_t []> (TIFFTileSize (tif));
+        uint8_t *tile_buf = tile_ptr.get ();
+        if (TIFFReadEncodedTile (tif, tile_no, tile_buf, -1) == -1)
+            error ("Failed to read tile data");
+        
+        // unpack tile bits into output matrix cells
+        for (uint32_t row = 0; row < tile_height; row++)
+          {
+            for (uint32_t col = 0; col < tile_width; col++)
+              {
+                uint8_t shift = 7 - col % 8;
+                tile_fvec[col] = (tile_buf [col / 8] >> shift) & 0x1;
+              }  
+            tile_fvec += tile_width;
+            tile_buf += tile_width / 8;
+          }
+      }
+    else
+      error ("Unsupported bit depth");
+    
+    Array<octave_idx_type> perm (dim_vector (3, 1));
+    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+      {
+        perm(0) = 2;
+        perm(1) = 1;
+        perm(2) = 0;
+      }
+    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      {
+        perm(0) = 1;
+        perm(1) = 0;
+        perm(2) = 2;
+      }
+    
+    tile_data = tile_data.permute (perm);
+
+    // Get the actual tile dimensions
+    uint32_t tiles_across = (image_data->width + tile_width - 1)
+                            / tile_width;
+    uint32_t tiles_down = (image_data->height + tile_height - 1)
+                          / tile_height;
+    uint32_t corrected_width = tile_width;
+    uint32_t corrected_height = tile_height;
+    if (tile_no % tiles_across == tiles_across - 1)
+      corrected_width = image_data->width - tile_width * (tiles_across - 1);
+    if ((tile_no / tiles_across) % tiles_down == tiles_down - 1)
+      corrected_height = image_data->height - tile_height * (tiles_down - 1);
+    
+    dim_vector corrected_dims;
+    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+      corrected_dims = dim_vector (corrected_height, corrected_width,
+                                   image_data->samples_per_pixel);
+    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      corrected_dims = dim_vector (corrected_height, corrected_width);
+    
+    tile_data.resize (corrected_dims);
+
+    return octave_value (tile_data);
+  }
+
+  template <typename T>
+  octave_value
+  read_strip_or_tile (TIFF *tif, uint32_t strip_tile_no,
+                      tiff_image_data *image_data)
+  {
+    if (image_data->is_tiled)
+      return read_tile<T> (tif, strip_tile_no, image_data);
+    else
+      return read_strip<T> (tif, strip_tile_no, image_data);
+  }
+
+  octave_value
+  handle_read_strip_or_tile (TIFF *tif, uint32_t strip_tile_no)
+  {
+    // Obtain all necessary tags
+    tiff_image_data image_data (tif);
+
+    uint16_t sample_format;
+    if (! TIFFGetFieldDefaulted(tif, TIFFTAG_SAMPLEFORMAT, &sample_format))
+      error ("Failed to obtain a value for sample format");
+
+    if (sample_format == 3)
+      {
+        if (image_data.bits_per_sample != 32
+            && image_data.bits_per_sample != 64)
+          error ("Floating point images are only supported for bit depths of 32 and 64");
+      }
+    else if (sample_format != 1 && sample_format != 4)
+      error ("Unsupported sample format");
+    
+    switch (image_data.bits_per_sample)
+      {
+      case 1:
+        return read_strip_or_tile<boolNDArray> (tif, strip_tile_no,
+                                                &image_data);
+        break;
+      case 8:
+        return read_strip_or_tile<uint8NDArray> (tif, strip_tile_no,
+                                                 &image_data);
+        break;
+      case 16:
+        return read_strip_or_tile<uint16NDArray> (tif, strip_tile_no,
+                                                  &image_data);
+        break;
+      case 32:
+        if (sample_format == 3)
+          return read_strip_or_tile<FloatNDArray> (tif, strip_tile_no,
+                                                   &image_data);
+        else
+          return read_strip_or_tile<uint32NDArray> (tif, strip_tile_no,
+                                                    &image_data);
+        break;
+      case 64:
+        if (sample_format == 3)
+          return read_strip_or_tile<NDArray> (tif, strip_tile_no,
+                                              &image_data);
+        else
+          return read_strip_or_tile<uint64NDArray> (tif, strip_tile_no,
+                                                    &image_data);
+        break;
+      default:
+        error ("Unsupported bit depth");
+      }
+  }
+
+  template <typename T>
+  octave_value
+  read_stripped_image (TIFF *tif, tiff_image_data *image_data)
+  {
+    typedef typename T::element_type P;
+
+    // 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
+    // size and the padding is removed at the end.
+    uint32_t padded_width = image_data->width;
+    uint8_t remove_padding = 0;
+    if ((image_data->bits_per_sample == 1 || image_data->bits_per_sample == 4)
+        && padded_width % 8 != 0)
+      {
+        padded_width += (8 - padded_width % 8) % 8;
+        remove_padding = 1;
+      }
+    
+    // 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));
+    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      img = T (dim_vector (padded_width, image_data->height,
+                           image_data->samples_per_pixel));
+    else
+      error ("Unsupported Planar Configuration");
+
+    P *img_fvec = img.fortran_vec ();
+
+    // Obtain the necessary data for handling the strips
+    uint32_t strip_count = TIFFNumberOfStrips (tif);
+    
+    // Can't rely on StripByteCounts because in compressed images
+    // the byte count reflect the actual number of bytes stored
+    // in the file not the size of the uncompressed strip
+    int64_t strip_size;
+    uint64_t written_size = 0;
+    uint64_t image_size = padded_width * image_data->height
+                          * image_data->samples_per_pixel
+                          * sizeof (P);
+    for (uint32_t strip = 0; strip < strip_count; strip++)
+      {
+        // Read the strip data into the matrix directly
+        strip_size = TIFFReadEncodedStrip (tif, strip, img_fvec, -1);
+
+        // Check if the strip read failed.
+        if (strip_size == -1)
+          error ("Failed to read strip data");
+        
+        // Check if the size being read exceeds the bounds of the matrix
+        // In case of a corrupt image with more data than needed
+        // This is probably redundant as LibTIFF checks sizes internally
+        if (written_size + strip_size > image_size)
+          error ("Strip data is larger than the image size");
+
+        if (image_data->bits_per_sample == 1)
+          {
+            if (image_data->samples_per_pixel != 1)
+              error ("Bi-Level images must have one channel only");
+            
+            // The strip size is multiplied by 8 to reflect tha actual
+            // number of bytes written to the matrix since each byte
+            // in the original strip contains 8 pixels of data
+            strip_size *= 8;
+
+            // Checking bounds again with the new size
+            if (written_size + strip_size > image_size)
+              error ("Strip 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 = strip_size - 1; pixel >= 0; pixel--)
+              {
+                // TODO(maged): is it necessary to check FillOrder?
+                uint8_t bit_number = 7 - pixel % 8;
+                uint8_t * img_u8 = reinterpret_cast<uint8_t *> (img_fvec);
+                img_fvec[pixel] = (img_u8[pixel / 8] >> bit_number) & 0x01;
+              }
+          }
+        else if (image_data->bits_per_sample == 4)
+          {
+            if (image_data->samples_per_pixel != 1)
+              error ("4-bit images are only supported for grayscale");
+            
+            // Strip size is multplied by as each byte contains 2 pixels
+            // and each pixel is later expanded into its own byte
+            strip_size *= 2;
+
+            // Checking bounds again with the ne size
+            if (written_size + strip_size > image_size)
+              error ("Strip 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 = strip_size - 1; pixel >= 0; pixel--)
+              {
+                uint8_t shift = pixel % 2 == 0? 4: 0;
+                uint8_t * img_u8 = reinterpret_cast<uint8_t *> (img_fvec);
+                img_fvec[pixel] = (img_u8[pixel / 2] >> shift) & 0x0F;
+              }
+          }
+        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");
+
+        // Advance the pointer by the amount of bytes read
+        img_fvec 
+          = reinterpret_cast<P *> (reinterpret_cast <uint8_t *> (img_fvec)
+                                   + strip_size);
+        written_size += strip_size;
+      }
+
+    // The matrices are permutated back to the shape expected by Octave
+    // which is height x width x channels
+    Array<octave_idx_type> perm (dim_vector (3, 1));
+    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+      {
+        perm(0) = 2;
+        perm(1) = 1;
+        perm(2) = 0;
+      }
+    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      {
+        perm(0) = 1;
+        perm(1) = 0;
+        perm(2) = 2;
+      }
+    
+    img = img.permute (perm);
+
+    if (remove_padding)
+      img.resize (dim_vector (image_data->height, image_data->width));
+
+    return octave_value (img);
+  }
+
+  template <typename T>
+  octave_value
+  read_tiled_image (TIFF *tif, tiff_image_data *image_data)
+  {
+    typedef typename T::element_type P;
+
+    // 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;
+
+    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;    
+    uint64_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 (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;
+                uint8_t * img_u8 = reinterpret_cast<uint8_t *> (img_fvec);
+                img_fvec[pixel]= (img_u8[pixel / 8] >> bit_number) & 0x01;
+              }
+          }
+        else if (image_data->bits_per_sample == 4)
+          {
+            if (image_data->samples_per_pixel != 1)
+              error ("4-bit images are only supported for grayscale");
+            
+            // 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;
+                uint8_t * img_u8 = reinterpret_cast<uint8_t *> (img_fvec);
+                img_fvec[pixel] = (img_u8[pixel / 2] >> shift) & 0x0F;
+              }
+          }
+        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
+          = reinterpret_cast<P *> (reinterpret_cast<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);
+
+        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);
+  }
+
+  template <typename T>
+  octave_value
+  read_image (TIFF *tif, tiff_image_data *image_data)
+  {
+    if (image_data->is_tiled)
+      return read_tiled_image<T> (tif, image_data);
+    else
+      return read_stripped_image<T> (tif, image_data);
+  }
+
+  // Convert tag value to double
+  octave_value
+  interpret_scalar_tag_data (void *data, TIFFDataType tag_datatype)
+  {
+    double retval;
+
+    switch (tag_datatype)
+      {
+      case TIFF_BYTE:
+      case TIFF_UNDEFINED:
+        {
+          retval = static_cast<double> (*(reinterpret_cast<uint8_t *> (data)));
+          break;
+        }
+      case TIFF_SHORT:
+        {
+          retval = static_cast<double> (*(reinterpret_cast<uint16_t *> (data)));
+          break;
+        }
+      case TIFF_LONG:
+        {
+          retval = static_cast<double> (*(reinterpret_cast<uint32_t *> (data)));
+          break;
+        }
+      case TIFF_LONG8:
+        {
+          retval = static_cast<double> (*(reinterpret_cast<uint64_t *> (data)));
+          break;
+        }
+      case TIFF_RATIONAL:
+        {
+          error ("TIFF_RATIONAL should have at least 2 elements but got only 1");
+          break;
+        }
+      case TIFF_SBYTE:
+        {
+          retval = static_cast<double> (*(reinterpret_cast<int8_t *> (data)));
+          break;
+        }
+      case TIFF_SSHORT:
+        {
+          retval = static_cast<double> (*(reinterpret_cast<int16_t *> (data)));
+          break;
+        }
+      case TIFF_SLONG:
+        {
+          retval = static_cast<double> (*(reinterpret_cast<int32_t *> (data)));
+          break;
+        }
+      case TIFF_SLONG8:
+        {
+          retval = static_cast<double> (*(reinterpret_cast<int64_t *> (data)));
+          break;
+        }
+      case TIFF_FLOAT:
+        {
+          retval = *(reinterpret_cast<float *> (data));
+          break;
+        }
+      case TIFF_DOUBLE:
+        {
+          retval = *(reinterpret_cast<double *> (data));
+          break;
+        }
+      case TIFF_SRATIONAL:
+        {
+          error ("TIFF_SRATIONAL should have at least 2 elements but got only 1");
+          break;
+        }
+      case TIFF_IFD:
+      case TIFF_IFD8:
+        error ("Unimplemented IFFD data type");
+        break;
+      default:
+        error ("Unsupported tag data type");
+      }
+
+    return octave_value (retval);
+  }
+
+  // Convert memory buffer into a suitable octave value
+  // depending on tag_datatype
+  octave_value
+  interpret_tag_data (void *data, uint32_t count, TIFFDataType tag_datatype)
+  {
+    octave_value retval;
+    // Apparently matlab converts scalar numerical values into double
+    // but doesn't do the same for arrays
+    if (count == 1 && tag_datatype != TIFF_ASCII)
+      {
+        retval = interpret_scalar_tag_data (data, tag_datatype);
+      }
+    else
+      {
+        dim_vector arr_dims (1, count);
+
+        switch (tag_datatype)
+          {
+          case TIFF_BYTE:
+          case TIFF_UNDEFINED:
+            {
+              uint8NDArray arr (arr_dims);
+              for (uint32_t i = 0; i < count; i++)
+                {
+                  arr(i) = (reinterpret_cast<uint8_t *> (data))[i];
+                }
+              retval = octave_value (arr);
+              break;
+            }
+          case TIFF_ASCII:
+            {
+              retval = octave_value (*(reinterpret_cast<char **> (data)));
+              break;
+            }
+          case TIFF_SHORT:
+            {
+              uint16NDArray arr (arr_dims);
+              for (uint32_t i = 0; i < count; i++)
+                {
+                  arr(i) = (reinterpret_cast<uint16_t *> (data))[i];
+                }
+              retval = octave_value (arr);
+              break;
+            }
+          case TIFF_LONG:
+            {
+              uint32NDArray arr (arr_dims);
+              for (uint32_t i = 0; i < count; i++)
+                {
+                  arr(i) = (reinterpret_cast<uint32_t *> (data))[i];
+                }
+              retval = octave_value (arr);
+              break;
+            }
+          case TIFF_LONG8:
+            {
+              uint64NDArray arr (arr_dims);
+              for (uint32_t i = 0; i < count; i++)
+                {
+                  arr(i) = (reinterpret_cast<uint64_t *> (data))[i];
+                }
+              retval = octave_value (arr);
+              break;
+            }
+          case TIFF_RATIONAL:
+            {
+              NDArray arr (arr_dims);
+              for (uint32_t i = 0; i < count; i+=2)
+                {
+                  arr(i / 2) = static_cast<float> ((reinterpret_cast<uint32_t *> (data))[i])
+                               / static_cast<float> ((reinterpret_cast<uint32_t *> (data))[i+1]);
+                }
+              retval = octave_value (arr);
+              break;
+            }
+          case TIFF_SBYTE:
+            {
+              int8NDArray arr (arr_dims);
+              for (uint32_t i = 0; i < count; i++)
+                {
+                  arr(i) = (reinterpret_cast<int8_t *> (data))[i];
+                }
+              retval = octave_value (arr);
+              break;
+            }
+          case TIFF_SSHORT:
+            {
+              int16NDArray arr (arr_dims);
+              for (uint32_t i = 0; i < count; i++)
+                {
+                  arr(i) = (reinterpret_cast<int16_t *> (data))[i];
+                }
+              retval = octave_value (arr);
+              break;
+            }
+          case TIFF_SLONG:
+            {
+              int32NDArray arr (arr_dims);
+              for (uint32_t i = 0; i < count; i++)
+                {
+                  arr(i) = (reinterpret_cast<int32_t *> (data))[i];
+                }
+              retval = octave_value (arr);
+              break;
+            }
+          case TIFF_SLONG8:
+            {
+              int64NDArray arr (arr_dims);
+              for (uint32_t i = 0; i < count; i++)
+                {
+                  arr(i) = (reinterpret_cast<int64_t *> (data))[i];
+                }
+              retval = octave_value (arr);
+              break;
+            }
+          case TIFF_FLOAT:
+            {
+              NDArray arr (arr_dims);
+              for (uint32_t i = 0; i < count; i++)
+                {
+                  arr(i) = (reinterpret_cast<float *> (data))[i];
+                }
+              retval = octave_value (arr);
+              break;
+            }
+          case TIFF_DOUBLE:
+            {
+              NDArray arr (arr_dims);
+              for (uint32_t i = 0; i < count; i++)
+              {
+                  arr(i) = (reinterpret_cast<double *> (data))[i];
+              }
+              retval = octave_value (arr);
+              break;
+            }
+          case TIFF_SRATIONAL:
+            {
+              NDArray arr (arr_dims);
+              for (uint32_t i = 0; i < count; i+=2)
+                {
+                  arr(i / 2) = static_cast<float> ((reinterpret_cast<int32_t *> (data))[i])
+                               / static_cast<float> ((reinterpret_cast<int32_t *> (data))[i+1]);
+                }
+              retval = octave_value (arr);
+              break;
+            }
+          case TIFF_IFD:
+          case TIFF_IFD8:
+            // TODO(maged): implement IFD datatype?
+            error ("Unimplemented IFD data type");
+            break;
+          default:
+            error ("Unsupported tag data type");
+          }
+      }
+
+    return retval;
+  }
+
+  octave_value
+  get_scalar_field_data (TIFF *tif, const TIFFField *fip)
+  {
+    uint32_t tag_id = TIFFFieldTag (fip);
+
+    // TIFFFieldReadCount returns VARIABLE for some scalar tags
+    // (e.g. Compression) But TIFFFieldPassCount seems consistent
+    // Since scalar tags are the last to be handled, any tag that
+    // require a count to be passed is an unsupported tag.
+    if (TIFFFieldPassCount (fip))
+      error ("Unsupported tag");
+    
+    int type_size = TIFFDataWidth (TIFFFieldDataType (fip));
+    
+    void *data = _TIFFmalloc (type_size);
+    if (tag_id == TIFFTAG_PLANARCONFIG)
+      {
+        // Workaround for a bug in LibTIFF where it incorrectly returns
+        // zero as the default value for PlanaConfiguration
+        if (! TIFFGetField(tif, tag_id, data))
+          *reinterpret_cast<uint16_t *> (data) = 1;
+      }
+    else
+      validate_tiff_get_field (TIFFGetFieldDefaulted (tif, tag_id, data), data);
+    
+    octave_value tag_data_ov = interpret_tag_data (data, 1,
+                                                   TIFFFieldDataType (fip));
+    _TIFFfree (data);
+
+    return tag_data_ov;
+  }
+
+  octave_value
+  get_array_field_data (TIFF *tif, const TIFFField *fip, uint32_t array_size)
+  {
+    void *data;
+    validate_tiff_get_field (TIFFGetField (tif, TIFFFieldTag (fip), &data));
+      
+    return interpret_tag_data (data, array_size, TIFFFieldDataType (fip));
+  }
+
+  octave_value
+  get_field_data (TIFF *tif, const TIFFField *fip)
+  {
+    octave_value tag_data_ov;
+    uint32_t tag_id = TIFFFieldTag (fip);
+
+    // TODO(maged): find/create images to test the special tags
+    switch (tag_id)
+      {
+      case TIFFTAG_STRIPBYTECOUNTS:
+      case TIFFTAG_STRIPOFFSETS:
+        tag_data_ov = get_array_field_data (tif, fip, 
+                                            TIFFNumberOfStrips (tif));
+        break;
+      case TIFFTAG_TILEBYTECOUNTS:
+      case TIFFTAG_TILEOFFSETS:
+        tag_data_ov = get_array_field_data (tif, fip, TIFFNumberOfTiles (tif));
+        break;
+      case TIFFTAG_YCBCRCOEFFICIENTS:
+        tag_data_ov = get_array_field_data (tif, fip, 3);
+        break;
+      case TIFFTAG_REFERENCEBLACKWHITE:
+        tag_data_ov = get_array_field_data (tif, fip, 6);
+        break;
+      case TIFFTAG_GRAYRESPONSECURVE:
+        {
+          uint16_t bits_per_sample;
+          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
+                                            &bits_per_sample))
+            error ("Failed to obtain the bit depth");
+          
+          tag_data_ov = get_array_field_data (tif, fip, 1<<bits_per_sample);
+          break;
+        }
+      case TIFFTAG_COLORMAP:
+        {
+          uint16_t bits_per_sample;
+          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
+                                       &bits_per_sample))
+            error ("Failed to obtain the bit depth");
+          
+          // According to the format specification, this field should
+          // be 8 or 16 only.
+          if (bits_per_sample > 16)
+            error ("Too high bit depth for a palette image");
+
+          uint32_t count = 1 << bits_per_sample;
+          uint16_t *red, *green, *blue;
+          validate_tiff_get_field (TIFFGetField (tif, TIFFTAG_COLORMAP,
+                                                 &red, &green, &blue));
+          
+          // Retrieving the data of the three channels and concatenating
+          // them together
+          OCTAVE_LOCAL_BUFFER (NDArray, array_list, 3);
+          dim_vector col_dims(count, 1);
+          array_list[0] = NDArray (interpret_tag_data (red,
+                                                       count,
+                                                       TIFFFieldDataType(fip))
+                                                       .uint16_array_value ()
+                                                       .reshape (col_dims));
+          array_list[1] = NDArray (interpret_tag_data (green,
+                                                       count,
+                                                       TIFFFieldDataType(fip))
+                                                       .uint16_array_value ()
+                                                       .reshape (col_dims));
+          array_list[2] = NDArray (interpret_tag_data (blue,
+                                                       count,
+                                                       TIFFFieldDataType(fip))
+                                                       .uint16_array_value ()
+                                                       .reshape (col_dims));
+
+          NDArray mat_out = NDArray::cat(1, 3, array_list);
+          // Normalize the range to be between 0 and 1
+          mat_out /= UINT16_MAX;
+
+          tag_data_ov = octave_value (mat_out);
+          break;
+        }
+      case TIFFTAG_TRANSFERFUNCTION:
+        {
+          uint16_t samples_per_pixel;
+          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_SAMPLESPERPIXEL,
+                                       &samples_per_pixel))
+            error ("Failed to obtain the number of samples per pixel");
+
+          uint16_t bits_per_sample;
+          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
+                                       &bits_per_sample))
+            error ("Failed to obtain the number of samples per pixel");
+
+          uint32_t count = 1 << bits_per_sample;
+          uint16_t *ch1, *ch2, *ch3;
+          if (samples_per_pixel == 1)
+            {
+              validate_tiff_get_field (TIFFGetField (tif, TIFFTAG_COLORMAP, &ch1));
+              tag_data_ov = interpret_tag_data (ch1, count,
+                                                TIFFFieldDataType (fip));
+            }
+          else
+            {
+              validate_tiff_get_field (TIFFGetField (tif, TIFFTAG_COLORMAP,
+                                                     &ch1, &ch2, &ch3));
+              
+              OCTAVE_LOCAL_BUFFER (uint16NDArray, array_list, 3);
+              dim_vector col_dims(count, 1);
+              array_list[0] = interpret_tag_data (ch1,
+                                                  count,
+                                                  TIFFFieldDataType (fip))
+                                                  .uint16_array_value ()
+                                                  .reshape (col_dims);
+              array_list[1] = interpret_tag_data (ch2,
+                                                  count,
+                                                  TIFFFieldDataType (fip))
+                                                  .uint16_array_value ()
+                                                  .reshape (col_dims);
+              array_list[2] = interpret_tag_data (ch3,
+                                                  count,
+                                                  TIFFFieldDataType (fip))
+                                                  .uint16_array_value ()
+                                                  .reshape (col_dims);
+              
+              tag_data_ov = octave_value (uint16NDArray::cat (1, 3, array_list));
+            }
+          break;
+        }
+      case TIFFTAG_PAGENUMBER:
+      case TIFFTAG_HALFTONEHINTS:
+      case TIFFTAG_DOTRANGE:
+      case TIFFTAG_YCBCRSUBSAMPLING:
+        {
+          uint16_t tag_part1, tag_part2;
+          validate_tiff_get_field (TIFFGetField (tif, tag_id,
+                                                 &tag_part1, &tag_part2));
+          
+          NDArray mat_out (dim_vector (1, 2));
+          mat_out(0)
+            = interpret_tag_data (&tag_part1, 1,
+                                  TIFFFieldDataType (fip)).double_value ();
+          mat_out(1)
+            = interpret_tag_data (&tag_part2, 1,
+                                  TIFFFieldDataType (fip)).double_value ();
+          
+          tag_data_ov = octave_value (mat_out);
+          break;
+        }
+      case TIFFTAG_SUBIFD:
+        {
+          uint16_t count;
+          uint64_t *offsets;
+          validate_tiff_get_field (TIFFGetField (tif, tag_id, &count, &offsets));
+          tag_data_ov = interpret_tag_data (offsets, count,
+                                            TIFFFieldDataType (fip));
+          break;
+        }
+      case TIFFTAG_EXTRASAMPLES:
+        {
+          uint16_t count;
+          uint16_t *types;
+          validate_tiff_get_field (TIFFGetField (tif, tag_id, &count, &types));
+          tag_data_ov = interpret_tag_data (types, count,
+                                            TIFFFieldDataType (fip));
+          break;
+        }
+      // TODO(maged): These tags are more complex to implement
+      //  will be implemented and tested later.
+      case TIFFTAG_XMLPACKET:
+      case TIFFTAG_RICHTIFFIPTC:
+      case TIFFTAG_PHOTOSHOP:
+      case TIFFTAG_ICCPROFILE:
+        {
+          error ("Complex Tags not implemented");
+          break;
+        }
+      // These tags are not mentioned in the LibTIFF documentation
+      // but are handled correctly by the library
+      case TIFFTAG_ZIPQUALITY:
+      case TIFFTAG_SGILOGDATAFMT:
+      case TIFFTAG_GRAYRESPONSEUNIT:
+        {
+          tag_data_ov = get_scalar_field_data (tif, fip);
+          break;
+        }
+      default:
+        tag_data_ov = get_scalar_field_data (tif, fip);
+    }
+    
+    return tag_data_ov;
+  }
+
+  void
+  set_field_data (TIFF *tif, const TIFFField *fip, octave_value tag_ov)
+  {
+    // TODO(maged): prevent calling here for read-only images
+    
+    // TODO(maged): complete the implementation of this function
+    uint32_t tag_id = TIFFFieldTag (fip);
+    uint32_t tag_data = tag_ov.double_value ();
+
+    if (! TIFFSetField(tif, tag_id, tag_data))
+      error ("Failed to set tag value");
+  }
+  
+  template <typename T>
+  void
+  write_strip (TIFF *tif, uint32_t strip_no, T strip_data,
+               tiff_image_data *image_data)
+  { 
+    uint32_t rows_in_strip;
+    if (! TIFFGetFieldDefaulted (tif, TIFFTAG_ROWSPERSTRIP, &rows_in_strip))
+      error ("Failed to obtain the RowsPerStrip tag");
+    
+    // The tag has a default value of UINT32_MAX which means the entire
+    // image, but we need to cap it to the height for later calculations
+    if (rows_in_strip > image_data->height)
+      rows_in_strip = image_data->height;
+    
+    // LibTIFF uses zero-based indexing as opposed to Octave's 1-based
+    strip_no--;
+
+    uint32_t strip_count = TIFFNumberOfStrips (tif);
+    rows_in_strip = get_rows_in_strip (strip_no, strip_count,
+                                       rows_in_strip, image_data);
+
+    dim_vector strip_dimensions;
+    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+      strip_dimensions = dim_vector (rows_in_strip, image_data->width,
+                                     image_data->samples_per_pixel);
+    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      strip_dimensions = dim_vector (rows_in_strip, image_data->width, 1);
+    else
+      error ("Planar configuration not supported");
+
+    if (strip_data.dim1 () > rows_in_strip)
+      warning ("The strip is composed of %u rows but the input has %ld rows.",
+               rows_in_strip,
+               strip_data.dim1 ());
+    
+    if (strip_data.dim2 () > image_data->width)
+      warning ("The image width is %u but the input has %ld columns.",
+               image_data->width,
+               strip_data.dim2 ());
+    
+    if (strip_data.ndims () > 2)
+      {
+        if (image_data->planar_configuration == PLANARCONFIG_CONTIG
+            && strip_data.dim3 () > image_data->samples_per_pixel)
+          warning ("The strip is composed of %u channels but the input has %ld channels.",
+                   image_data->samples_per_pixel,
+                   strip_data.dim3 ());
+        else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE
+                && strip_data.dim3 () > 1)
+          warning ("The strip is composed of %u channel but the input has %ld channels.",
+                   1, strip_data.dim3 ());
+      }
+
+    strip_data.resize (strip_dimensions);
+
+    // Permute the dimesions of the strip to match the expected memory
+    // arrangement of LibTIFF (channel x width x height)
+    Array<octave_idx_type> perm (dim_vector (3, 1));
+    perm(0) = 2;
+    perm(1) = 1;
+    perm(2) = 0;
+    strip_data = strip_data.permute (perm);
+
+    void *data_vec = strip_data.fortran_vec ();
+    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)
+      {
+        // Can't rely on LibTIFF's TIFFStripSize because boundary strips
+        // can be smaller in size
+        tsize_t strip_size = strip_data.numel ()
+                             * image_data->bits_per_sample / 8;
+        if (TIFFWriteEncodedStrip (tif, strip_no, data_vec, strip_size) == -1)
+          error ("Failed to write strip data to image");
+        
+      }
+    else if (image_data->bits_per_sample == 1)
+      {
+        if (image_data->samples_per_pixel != 1)
+          error ("Bi-Level images must have one channel only");
+        
+        // Create a buffer to hold the packed strip data
+        // Unique pointers are faster than vectors for constant size buffers
+        std::unique_ptr<uint8_t []> strip_ptr
+          = std::make_unique<uint8_t []> (TIFFStripSize (tif));
+        uint8_t *strip_buf = strip_ptr.get ();
+        uint8_t *data_u8 = reinterpret_cast<uint8_t *> (data_vec);
+        // According to the format specification, the row should be byte
+        // aligned so the number of bytes is rounded up to the nearest byte
+        uint32_t padded_width = (image_data->width + 7) / 8;
+        // Packing the pixel data into bits
+        for (uint32_t row = 0; row < rows_in_strip; row++)
+          {
+            for (uint32_t col = 0; col < image_data->width; col++)
+            {
+              uint8_t shift = 7 - col % 8;
+              strip_buf[row * padded_width + col / 8] |= data_u8[col] << shift;
+            }
+            data_u8 += image_data->width;
+          }
+        tsize_t strip_size = padded_width * rows_in_strip;
+        if (TIFFWriteEncodedStrip (tif, strip_no, strip_buf, strip_size) == -1)
+          error ("Failed to write strip data to image");
+      }
+    else
+      {
+        error ("Unsupported bit depth");
+      }
+  }
+
+  template <typename T>
+  void
+  write_tile (TIFF *tif, uint32_t tile_no, T tile_data,
+              tiff_image_data *image_data)
+  {
+    uint32_t tile_width, tile_height;
+    if (! TIFFGetField (tif, TIFFTAG_TILEWIDTH, &tile_width))
+      error ("Failed to get the tile width");
+    if (! TIFFGetField (tif, TIFFTAG_TILELENGTH, &tile_height))
+      error ("Failed to get the tile length");
+
+    if (tile_height == 0 || tile_height % 16 != 0
+        || tile_width == 0 || tile_width % 16 != 0)
+      error ("Tile dimesion tags are invalid");
+    
+    if (tile_no < 1 || tile_no > TIFFNumberOfTiles (tif))
+      error ("Tile number out of bounds");
+
+    if (tile_data.dim1 () > tile_height)
+      warning ("The tile is composed of %u rows but input has %ld rows",
+               tile_height, tile_data.dim1 ());
+    if (tile_data.dim2 () > tile_width)
+      warning ("The tile is composed of %u columns but input has %ld columns",
+               tile_width, tile_data.dim2 ());
+    if (tile_data.ndims () > 2)
+      {
+        if (image_data->planar_configuration == PLANARCONFIG_CONTIG
+            && tile_data.dim3 () > image_data->samples_per_pixel)
+          warning ("The tile is composed of %u channels but input has %ld channels",
+                  image_data->samples_per_pixel, tile_data.dim3 ());
+        else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE
+                 && tile_data.dim3 () > 1)
+          warning ("The tile is composed of %u channels but input has %ld channels",
+                   1, tile_data.dim3 ());
+      }
+    
+    dim_vector tile_dimensions;
+    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+      tile_dimensions = dim_vector (tile_height, tile_width,
+                                    image_data->samples_per_pixel);
+    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      tile_dimensions = dim_vector (tile_height, tile_width, 1);
+    else
+      error ("Planar configuration not supported");
+    
+    tile_data.resize (tile_dimensions);
+    Array<octave_idx_type> perm (dim_vector (3, 1));
+    perm(0) = 2;
+    perm(1) = 1;
+    perm(2) = 0;
+    tile_data = tile_data.permute (perm);
+    
+    // Octave indexing is 1-based while LibTIFF is zero-based
+    tile_no--;
+    void *data_vec = tile_data.fortran_vec ();
+    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)
+      {
+        if (TIFFWriteEncodedTile (tif, tile_no, data_vec,
+                                  TIFFTileSize (tif)) == -1)
+          error ("Failed to write tile data to image");
+        
+      }
+    else if (image_data->bits_per_sample == 1)
+      {
+        if (image_data->samples_per_pixel != 1)
+          error ("Bi-Level images must have one channel only");
+        
+        // Create a buffer to hold the packed tile data
+        // Unique pointers are faster than vectors for constant size buffers
+        std::unique_ptr<uint8_t []> tile_ptr
+          = std::make_unique<uint8_t []> (TIFFTileSize (tif));
+        uint8_t *tile_buf = tile_ptr.get ();
+        uint8_t *data_u8 = reinterpret_cast<uint8_t *> (data_vec);
+        // Packing the pixel data into bits
+        for (uint32_t row = 0; row < tile_height; row++)
+          {
+            for (uint32_t col = 0; col < tile_width; col++)
+            {
+              uint8_t shift = 7 - col % 8;
+              tile_buf[row * tile_width/8 + col/8] |= data_u8[col] << shift;
+            }
+            data_u8 += tile_width;
+          }
+        if (TIFFWriteEncodedTile (tif, tile_no, tile_buf,
+                                  TIFFTileSize (tif)) == -1)
+          error ("Failed to write tile data to image");
+      }
+    else
+      {
+        error ("Unsupported bit depth");
+      }
+  }
+
+  template <typename T>
+  void
+  write_strip_or_tile (TIFF *tif, uint32_t strip_tile_no, T strip_data,
+                       tiff_image_data *image_data)
+  {
+    if (image_data->is_tiled)
+      write_tile<T> (tif, strip_tile_no, strip_data, image_data);
+    else
+      write_strip<T> (tif, strip_tile_no, strip_data, image_data);
+  }
+
+  void
+  handle_write_strip_or_tile (TIFF *tif, uint32_t strip_tile_no,
+                              octave_value data_ov,
+                              tiff_image_data *image_data)
+  {
+
+    // SampleFormat tag is not a required field and has a default value of 1
+    // So we need to use TIFFGetFieldDefaulted in case it is not present in
+    // the file
+    uint16_t sample_format;
+    if (! TIFFGetFieldDefaulted(tif, TIFFTAG_SAMPLEFORMAT, &sample_format))
+      error ("Failed to obtain a value for sample format");
+
+    // TODO(maged): add support for signed integer images
+    if (sample_format == 3)
+      {
+        if (image_data->bits_per_sample != 32
+            && image_data->bits_per_sample != 64)
+          error ("Floating point images are only supported for bit depths of 32 and 64");
+      }
+
+    // The standard specifies that a SampleFormat of 4 should be treated
+    // the same as 1 (unsigned integer)
+    else if (sample_format != 1 && sample_format != 4)
+      error ("Unsupported sample format");
+    
+    switch (image_data->bits_per_sample)
+      {
+      case 1:
+        // We need to check for both scalar and matrix types to handle single
+        // element strip
+        if (data_ov.is_bool_scalar () || data_ov.is_bool_matrix ())
+          write_strip_or_tile<boolNDArray> (tif, strip_tile_no,
+                                            data_ov.bool_array_value (),
+                                            image_data);
+        else
+          error ("Expected logical matrix for BiLevel image");
+        break;
+      case 8:
+        if (data_ov.is_uint8_type ())
+          write_strip_or_tile<uint8NDArray> (tif, strip_tile_no,
+                                             data_ov.uint8_array_value (),
+                                             image_data);
+        else
+          error ("Only uint8 data is allowed for uint images with bit depth of 8");
+        break;
+      case 16:
+        if (data_ov.is_uint16_type ())
+          write_strip_or_tile<uint16NDArray> (tif, strip_tile_no,
+                                              data_ov.uint16_array_value (),
+                                              image_data);
+        else
+          error ("Only uint16 data is allowed for uint images with bit depth of 16");
+        break;
+      case 32:
+        if (sample_format == 3)
+          if (data_ov.is_single_type () || data_ov.is_double_type ())
+            write_strip_or_tile<FloatNDArray> (tif, strip_tile_no,
+                                               data_ov.float_array_value (),
+                                               image_data);
+          else
+            error ("Only single and double data are allowed for floating-point images");
+        else
+          if (data_ov.is_uint32_type ())
+            write_strip_or_tile<uint32NDArray> (tif, strip_tile_no,
+                                                data_ov.uint32_array_value (),
+                                                image_data);
+          else
+            error ("Only uint32 data is allowed for uint images with bit depth of 32");
+        break;
+      case 64:
+        if (sample_format == 3)
+          if (data_ov.is_single_type () || data_ov.is_double_type ())
+            write_strip_or_tile<NDArray> (tif, strip_tile_no,
+                                          data_ov.array_value (),
+                                          image_data);
+          else
+            error ("Only single and double data are allowed for floating-point images");
+        else  
+          if (data_ov.is_uint64_type ())
+            write_strip_or_tile<uint64NDArray> (tif, strip_tile_no,
+                                                data_ov.uint64_array_value (),
+                                                image_data);
+          else
+            error ("Only uint64 data is allowed for uint images with bit depth of 64");
+        break;
+      default:
+        error ("Unsupported bit depth");
+      }
+  }
+
+  template <typename T>
+  void
+  write_stripped_image (TIFF *tif, T pixel_data, tiff_image_data *image_data)
+  {
+    // TODO(maged): remove this? ASSUMES pixel data dimensions are already validated
+
+    typedef typename T::element_type P;
+
+    // Permute pixel data to be aligned in memory to the way LibTIFF
+    // expects the data to be (i.e. channel x width x height for chunky
+    // and width x height x channel for separate planes)
+    Array<octave_idx_type> perm (dim_vector (3, 1));
+    if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      {
+        perm(0) = 1;
+        perm(1) = 0;
+        perm(2) = 2;
+      }
+    else
+      {
+        perm(0) = 2;
+        perm(1) = 1;
+        perm(2) = 0;
+      }
+    pixel_data = pixel_data.permute (perm);
+
+    uint32_t row_per_strip;
+    if (! TIFFGetFieldDefaulted (tif, TIFFTAG_ROWSPERSTRIP, &row_per_strip))
+      error ("Failed to obtain the RowPerStrip tag");
+    
+    // The default value is INT_MAX so we need to cap it to the image height
+    if (row_per_strip > image_data->height)
+      row_per_strip = image_data->height;
+
+    uint8_t *pixel_fvec = reinterpret_cast<uint8_t *> (pixel_data.fortran_vec ());
+    uint32_t strip_count = TIFFNumberOfStrips (tif);
+    tsize_t strip_size;
+    uint32_t rows_in_strip;
+    for (uint32_t strip = 0; strip < strip_count; strip++)
+      {
+        rows_in_strip = get_rows_in_strip (strip, strip_count,
+                                           row_per_strip, image_data);
+        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)
+          {
+            strip_size = rows_in_strip * image_data->width * sizeof (P);
+            if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
+              strip_size *= image_data->samples_per_pixel;
+            if (! TIFFWriteEncodedStrip (tif, strip, pixel_fvec, strip_size))
+              error ("Failed to rite strip data");
+            pixel_fvec += strip_size;
+          }
+        else if (image_data->bits_per_sample == 1)
+          {
+            if (image_data->samples_per_pixel != 1)
+              error ("Bi-Level images must have one channel only");
+            
+            // Create a buffer to hold the packed strip data
+            // Unique pointers are faster than vectors for constant size buffers
+            std::unique_ptr<uint8_t []> strip_ptr
+              = std::make_unique<uint8_t []> (TIFFStripSize (tif));
+            uint8_t *strip_buf = strip_ptr.get ();
+            // According to the format specification, the row should be byte
+            // aligned so the number of bytes is rounded up to the nearest byte
+            uint32_t padded_width = (image_data->width + 7) / 8;
+            // Packing the pixel data into bits
+            for (uint32_t row = 0; row < rows_in_strip; row++)
+              {
+                for (uint32_t col = 0; col < image_data->width; col++)
+                {
+                  uint8_t shift = 7 - col % 8;
+                  strip_buf[row * padded_width + col / 8]
+                    |= pixel_fvec[col] << shift;
+                }
+                pixel_fvec += image_data->width;
+              }
+            strip_size = padded_width * rows_in_strip;
+            if (TIFFWriteEncodedStrip (tif, strip, strip_buf, strip_size) == -1)
+              error ("Failed to write strip data to image");
+          }
+        else
+          error ("Unsupported bit depth");
+      }
+  }
+
+  template <typename T>
+  void
+  write_tiled_image (TIFF *tif, T pixel_data, tiff_image_data *image_data)
+  {
+    // TODO(maged): remove this? ASSUMES pixel data dimensions are already validated
+
+    uint32_t tile_width, tile_height;
+    if (! TIFFGetField (tif, TIFFTAG_TILEWIDTH, &tile_width))
+      error ("Failed to get the tile width");
+    if (! TIFFGetField (tif, TIFFTAG_TILELENGTH, &tile_height))
+      error ("Failed to get the tile length");
+
+    if (tile_height == 0 || tile_height % 16 != 0
+        || tile_width == 0 || tile_width % 16 != 0)
+      error ("Tile dimesion tags are invalid");
+    
+    uint32_t tiles_across = (image_data->width + tile_width - 1)
+                            / tile_width;
+    uint32_t tiles_down = (image_data->height + tile_height - 1)
+                          / tile_height;
+    
+    // Resize the image data to add tile padding
+    dim_vector padded_dims (tiles_down * tile_height,
+                            tiles_across * tile_width,
+                            image_data->samples_per_pixel);
+    pixel_data.resize (padded_dims);
+
+    // Reshape the image to separate tiles into their own dimension
+    // so we can permute them to the right order
+    dim_vector tiled_dims (tile_height, tiles_down, tile_width, tiles_across,
+                           image_data->samples_per_pixel);
+    pixel_data = pixel_data.reshape (tiled_dims);
+
+    // Permute the dimesnions to get the memory alignment to match LibTIFF
+    Array<octave_idx_type> perm (dim_vector (5, 1));
+    if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
+      {
+        // For separate planes, the data coming from libtiff will have all
+        // tiles of the first sample then all tiles of the second sample
+        // and so on. Tiles of each sample will be ordered from left to right
+        // and from top to bottom. And data inside each tile is organised as
+        // rows and each row contains columns.
+        // So the order for LibTIFF is:
+        //   samples x tiles_down x tiles_across x tile_height x tile_width
+        // But memory orientation of Octave arrays is reversed so we set it to
+        //   tile_width x tile_height x tiles_across x tiles_down x samples
+        perm(0) = 2;
+        perm(1) = 0;
+        perm(2) = 3;
+        perm(3) = 1;
+        perm(4) = 4;
+      }
+    else
+      {
+        // For chunky planes, the data coming from libtiff will be ordered
+        // from left to right and from top to bottom. And data inside each
+        // tile is organised as rows and each row contains columns and each
+        // column contains samples.
+        // So the order for LibTIFF is:
+        //   tiles_down x tiles_across x tile_height x tile_width x samples
+        // But memory orientation of Octave arrays is reversed so we set it to
+        //   samples x tile_width x tile_height x tiles_across x tiles_down
+        perm(0) = 4;
+        perm(1) = 2;
+        perm(2) = 0;
+        perm(3) = 3;
+        perm(4) = 1;
+      }
+    pixel_data = pixel_data.permute (perm);
+
+    uint8_t *pixel_fvec
+      = reinterpret_cast<uint8_t *> (pixel_data.fortran_vec ());
+    uint32_t tile_count = TIFFNumberOfTiles (tif);
+    tsize_t tile_size = TIFFTileSize (tif);
+
+    for (uint32_t tile = 0; tile <tile_count; tile++)
+      {
+        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)
+          {
+            if (! TIFFWriteEncodedTile (tif, tile, pixel_fvec, tile_size))
+              error ("Failed to write tile data");
+            pixel_fvec += tile_size;
+          }
+        else if (image_data->bits_per_sample == 1)
+          {
+            if (image_data->samples_per_pixel != 1)
+              error ("Bi-Level images must have one channel only");
+            
+            // Create a buffer to hold the packed tile data
+            // Unique pointers are faster than vectors for
+            // constant size buffers
+            std::unique_ptr<uint8_t []> tile_ptr
+              = std::make_unique<uint8_t []> (tile_size);
+            uint8_t *tile_buf = tile_ptr.get ();
+            // Packing the pixel data into bits
+            for (uint32_t row = 0; row < tile_height; row++)
+              {
+                for (uint32_t col = 0; col < tile_width; col++)
+                {
+                  uint8_t shift = 7 - col % 8;
+                  tile_buf[row * tile_width / 8 + col / 8]
+                    |= pixel_fvec[col] << shift;
+                }
+                pixel_fvec += tile_width;
+              }
+            if (TIFFWriteEncodedTile (tif, tile, tile_buf,
+                                      TIFFTileSize (tif)) == -1)
+              error ("Failed to write tile data to image");
+          }
+        else
+          error ("Unsupported bit depth");
+      }
+
+  }
+
+  template <typename T>
+  void
+  write_image (TIFF *tif, T pixel_data, tiff_image_data *image_data)
+  {
+    // TODO(maged): matlab sets the remaining to zero for less width and height
+    // and issues a warning for larger widths but not for larger height? (we should error)
+    // and errors for less or more number of channels
+    if (image_data->height != pixel_data.dim1 ()
+        || image_data->width != pixel_data.dim2 ()
+        || pixel_data.ndims () > 3
+        || (image_data->samples_per_pixel > 1 && pixel_data.ndims () < 3)
+        || (pixel_data.ndims () > 2
+            && image_data->samples_per_pixel != pixel_data.dim3 ()))
+      error ("Dimensions of the input don't match image dimenions");
+    
+    if (image_data->is_tiled)
+      write_tiled_image<T> (tif, pixel_data, image_data);
+    else
+      write_stripped_image<T> (tif, pixel_data, image_data);
+
+  }
+
+
+#endif
+
+  DEFUN (__open_tiff__, args, ,
+             "Open a Tiff file and return its handle")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    if (nargin == 0 || nargin > 2)
+      {
+        // TODO(maged): return invalid object instead??
+        error ("No filename supplied\n");
+      }
+
+    std::string filename = args(0).string_value ();
+    std::string mode = "r";
+
+    if (nargin == 2)
+      mode = args(1).string_value ();
+
+    const std::vector<std::string> supported_modes {"r", "w", "w8", "a"};
+      
+    if (std::find (supported_modes.cbegin (), supported_modes.cend (), mode)
+          == supported_modes.cend ())
+      {
+        if (mode == "r+")
+          error ("Openning files in r+ mode is not yet supported");
+        else
+          error ("Invalid mode for openning Tiff file: %s", mode.c_str ());
+      }
+    
+    TIFF *tif = TIFFOpen (filename.c_str (), mode.c_str ());
+    
+    if (! tif)
+      error ("Failed to open Tiff file\n");
+
+    // TODO(maged): use inheritance of octave_base_value instead
+    octave_value tiff_ov = octave_value ((uint64_t)tif);
+    return octave_value_list (tiff_ov);
+#else
+    octave_unused_parameter (args);
+    err_disabled_feature ("Tiff", "Tiff");
+#endif
+  }
+
+
+  DEFUN (__close_tiff__, args, ,
+             "Close a tiff file")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    if (nargin == 0)
+      error ("No handle provided\n");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+    TIFFClose (tif);
+
+    return octave_value_list ();
+#else
+    err_disabled_feature ("close", "Tiff");
+#endif
+  }
+
+
+  DEFUN (__tiff_get_tag__, args, ,
+             "Get the value of a tag from a tiff image")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    if (nargin == 0)
+      error ("No handle provided\n");
+    
+    if (nargin < 2)
+      error ("No tag name provided\n");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+
+    uint32_t tag_id;
+    const TIFFField *fip;
+    
+    if (args(1).is_string ())
+      {
+        std::string tagName = args(1).string_value ();
+        fip = TIFFFieldWithName (tif, tagName.c_str ());
+        if (! fip)
+          error ("Tiff tag not found");
+        
+        tag_id = TIFFFieldTag (fip);
+      }
+    else
+      {
+        tag_id = args(1).int_value ();
+        fip = TIFFFieldWithTag (tif, tag_id);
+        
+        if (! fip)
+          error ("Tiff tag not found");
+      }
+
+    return octave_value_list (get_field_data (tif, fip));
+#else
+    err_disabled_feature ("getTag", "Tiff");
+#endif
+  }
+
+
+  DEFUN (__tiff_set_tag__, args, ,
+             "Set the value of a tag in a tiff image")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+    
+    if (nargin < 2)
+      error ("Too few arguments provided\n");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+    
+    if (args(1).isstruct ())
+      {
+        octave_scalar_map tags = args(1).xscalar_map_value ("__tiff_set_tag__: struct argument must be a scalar structure");
+        string_vector keys = tags.fieldnames ();
+        // Using iterators instead of this loop method seems to process
+        // the elements of the struct in a different order than they were
+        // assigned which makes the behavior here incompatible with matlab
+        // in case of errors.
+        for (octave_idx_type i = 0; i < keys.numel (); i++)
+          {
+            std::string key = keys[i];
+            const TIFFField *fip =  TIFFFieldWithName (tif, key.c_str ());
+            if (! fip)
+              error ("Tag %s not found", key.c_str ());
+            set_field_data (tif, fip, tags.contents (key));
+          }
+      }
+    else
+      {
+        if (nargin < 3)
+          error ("Too few arguments provided");
+
+        const TIFFField *fip;
+        // TODO(maged): matlab actually checks for its own strings not LibTIFF's
+        // Make it work in both ways
+        if (args(1).is_string ())
+          {
+            std::string tagName = args(1).string_value ();
+            fip = TIFFFieldWithName (tif, tagName.c_str ());
+            if (! fip)
+              error ("Tiff tag not found");
+          }
+        else
+          {
+            uint32_t tag_id = args(1).int_value ();
+            fip = TIFFFieldWithTag (tif, tag_id);
+            if (! fip)
+              error ("Tiff tag not found");
+          }
+        
+        set_field_data (tif, fip, args(2));
+      }
+
+    return octave_value_list ();
+#else
+    err_disabled_feature ("setTag", "Tiff");
+#endif
+  }
+
+  DEFUN (__tiff_read__, args, nargout,
+             "Read the image in the current IFD")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    if (nargin == 0)
+      error ("No handle provided\n");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+
+    // TODO(maged): nargout and ycbcr
+    octave_unused_parameter (nargout);
+
+    // Obtain all necessary tags
+    // The main purpose of this struct is to hold all the necessary tags that
+    // will be common among all the different cases of handling the the image
+    // data to avoid repeating the same calls to TIFFGetField in the different
+    // functions as each call is a possible point of failure
+    tiff_image_data image_data (tif);
+
+    uint16_t sample_format;
+    if (! TIFFGetFieldDefaulted(tif, TIFFTAG_SAMPLEFORMAT, &sample_format))
+      error ("Failed to obtain a value for sample format");
+
+    if (sample_format == 3)
+      {
+        if (image_data.bits_per_sample != 32 && image_data.bits_per_sample != 64)
+          error ("Floating point images are only supported for bit depths of 32 and 64");
+      }
+    else if (sample_format != 1 && sample_format != 4)
+      error ("Unsupported sample format");
+    
+    octave_value_list retval;
+    switch (image_data.bits_per_sample)
+      {
+      case 1:
+        retval(0) = read_image<boolNDArray> (tif, &image_data);
+        break;
+      case 4:
+      case 8:
+        retval(0) = read_image<uint8NDArray> (tif, &image_data);
+        break;
+      case 16:
+        retval(0) = read_image<uint16NDArray> (tif, &image_data);
+        break;
+      case 32:
+        if (sample_format == 3)
+          retval(0) = read_image<FloatNDArray> (tif, &image_data);
+        else
+          retval(0) = read_image<uint32NDArray> (tif, &image_data);
+        break;
+      case 64:
+        if (sample_format == 3)
+          retval(0) = read_image<NDArray> (tif, &image_data);
+        else
+          retval(0) = read_image<uint64NDArray> (tif, &image_data);
+        break;
+      default:
+        error ("Unsupported bit depth");
+      }
+    
+    return retval;
+#else
+    err_disabled_feature ("read", "Tiff");
+#endif
+  }
+
+  DEFUN (__tiff_read_encoded_strip__, args, nargout,
+             "Read a strip from the image in the current IFD")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    if (nargin != 2)
+      error ("rong number of arguments");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+
+    // TODO(maged): nargout and ycbcr
+    octave_unused_parameter (nargout);
+
+    if (TIFFIsTiled (tif))
+      error ("The image is tiled not stripped");
+    
+    uint32_t strip_no;
+    if (args(1).is_scalar_type ())
+      strip_no = args(1).uint32_scalar_value ();
+    else
+      error ("Expected scalar for strip number");
+    
+    if (strip_no < 1 || strip_no > TIFFNumberOfStrips (tif))
+      error ("Strip number out of bounds");
+    
+    // Convert from Octave's 1-based indexing to zero-based indexing
+    strip_no--;
+
+    return octave_value_list (handle_read_strip_or_tile (tif, strip_no));
+#else
+    err_disabled_feature ("readEncodedStrip", "Tiff");
+#endif
+  }
+
+  DEFUN (__tiff_read_encoded_tile__, args, nargout,
+             "Read a tile from the image in the current IFD")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    if (nargin != 2)
+      error ("rong number of arguments");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+
+    // TODO(maged): nargout and ycbcr
+    octave_unused_parameter (nargout);
+
+    if (! TIFFIsTiled (tif))
+      error ("The image is stripped not tiled");
+    
+    uint32_t tile_no;
+    if (args(1).is_scalar_type ())
+      tile_no = args(1).uint32_scalar_value ();
+    else
+      error ("Expected scalar for tile number");
+    
+    if (tile_no < 1 || tile_no > TIFFNumberOfTiles (tif))
+      error ("Tile number out of bounds");
+    
+    // Convert from Octave's 1-based indexing to zero-based indexing
+    tile_no--;
+
+    return octave_value_list (handle_read_strip_or_tile (tif, tile_no));
+#else
+    err_disabled_feature ("readEncodedTile", "Tiff");
+#endif
+  }
+
+  DEFUN (__tiff_write__, args, ,
+             "Write the image data to the current IFD")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    if (nargin < 2)
+      error ("Wrong number of arguments\n");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+
+    // TODO(maged): check on windows
+    if (TIFFGetMode (tif) == O_RDONLY)
+      error ("Can't write data to a file opened in read-only mode");
+
+    // Obtain all necessary tags
+    tiff_image_data image_data (tif);
+
+    uint16_t sample_format;
+    if (! TIFFGetFieldDefaulted(tif, TIFFTAG_SAMPLEFORMAT, &sample_format))
+      error ("Failed to obtain a value for sample format");
+
+    if (sample_format == 3)
+      {
+        if (image_data.bits_per_sample != 32
+            && image_data.bits_per_sample != 64)
+          error ("Floating point images are only supported for bit depths of 32 and 64");
+      }
+    else if (sample_format != 1 && sample_format != 4)
+      error ("Unsupported sample format");
+    
+    switch (image_data.bits_per_sample)
+      {
+      case 1:
+        // We need to check for both scalar and matrix types to handle single
+        // pixel image
+        if (args (1).is_bool_scalar () || args (1).is_bool_matrix ())
+          write_image<boolNDArray> (tif, args (1).bool_array_value (),
+                                    &image_data);
+        else
+          error ("Expected logical matrix for BiLevel image");
+        break;
+      case 8:
+        if (args (1).is_uint8_type ())
+          write_image<uint8NDArray> (tif, args (1).uint8_array_value (),
+                                     &image_data);
+        else
+          error ("Only uint8 data is allowed for uint images with bit depth of 8");
+        break;
+      case 16:
+        if (args (1).is_uint16_type ())
+          write_image<uint16NDArray> (tif, args (1).uint16_array_value (),
+                                      &image_data);
+        else
+          error ("Only uint16 data is allowed for uint images with bit depth of 16");
+        break;
+      case 32:
+        if (sample_format == 3)
+          if (args (1).is_single_type () || args (1).is_double_type ())
+            write_image<FloatNDArray> (tif, args (1).float_array_value (),
+                                       &image_data);
+          else
+            error ("Only single and double data are allowed for floating-point images");
+        else
+          if (args (1).is_uint32_type ())
+            write_image<uint32NDArray> (tif, args (1).uint32_array_value (),
+                                        &image_data);
+          else
+            error ("Only uint32 data is allowed for uint images with bit depth of 32");
+        break;
+      case 64:
+        if (sample_format == 3)
+          if (args (1).is_single_type () || args (1).is_double_type ())
+            write_image<NDArray> (tif, args (1).array_value (), &image_data);
+          else
+            error ("Only single and double data are allowed for floating-point images");
+        else  
+          if (args (1).is_uint64_type ())
+            write_image<uint64NDArray> (tif, args (1).uint64_array_value (),
+                                        &image_data);
+          else
+            error ("Only uint64 data is allowed for uint images with bit depth of 64");
+        break;
+      default:
+        error ("Unsupported bit depth");
+      }
+    
+    return octave_value_list ();
+#else
+    err_disabled_feature ("write", "Tiff");
+#endif
+  }
+
+  DEFUN (__tiff_write_encoded_strip__, args, ,
+             "Write an encoded strip to the image")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    // TODO(maged): add support for YCbCr data
+    if (nargin < 3)
+      error ("Too few arguments provided\n");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+
+    // TODO(maged): check on windows
+    if (TIFFGetMode (tif) == O_RDONLY)
+      error ("Can't write data to a file opened in read-only mode");
+
+    // Obtain all necessary tags
+    tiff_image_data image_data (tif);
+
+    if (image_data.is_tiled)
+      error ("Can't write strips to a tiled image");
+
+    uint32_t strip_no = args(1).uint32_scalar_value ();
+    if (strip_no < 1 || strip_no > TIFFNumberOfStrips (tif))
+      error ("Strip number out of range");
+    
+    handle_write_strip_or_tile (tif, strip_no, args(2), &image_data);
+    
+    return octave_value_list ();
+#else
+    err_disabled_feature ("writeEncodedStrip", "Tiff");
+#endif
+  }
+
+  DEFUN (__tiff_write_encoded_tile__, args, ,
+             "Write an encoded tile to the image")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    // TODO(maged): add support for YCbCr data
+    if (nargin < 3)
+      error ("Too few arguments provided\n");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+
+    // TODO(maged): check on windows
+    if (TIFFGetMode (tif) == O_RDONLY)
+      error ("Can't write data to a file opened in read-only mode");
+
+    // Obtain all necessary tags
+    tiff_image_data image_data (tif);
+
+    if (! image_data.is_tiled)
+      error ("Can't write tiles to a stripped image");
+
+    uint32_t tile_no = args(1).uint32_scalar_value ();
+    if (tile_no < 1 || tile_no > TIFFNumberOfTiles (tif))
+      error ("Tile number out of range");
+
+    handle_write_strip_or_tile (tif, tile_no, args(2), &image_data);
+    
+    return octave_value_list ();
+#else
+    err_disabled_feature ("writeEncodedTile", "Tiff");
+#endif
+  }
+  
+  DEFUN (__tiff_is_tiled__, args, ,
+             "Get whether the image is tiled")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    if (nargin == 0)
+      error ("No handle provided\n");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+    bool is_tiled = static_cast<bool> (TIFFIsTiled (tif));
+    return octave_value_list (octave_value (is_tiled));
+#else
+    err_disabled_feature ("isTiled", "Tiff");
+#endif
+  }
+
+  DEFUN (__tiff_number_of_strips__, args, ,
+             "Get the number of strips in the image")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    if (nargin == 0)
+      error ("No handle provided\n");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+    
+    if (TIFFIsTiled (tif))
+      error ("The image is tiled not stripped");
+    
+    double strip_count = static_cast<double> (TIFFNumberOfStrips (tif));
+    return octave_value_list (octave_value (strip_count));
+#else
+    err_disabled_feature ("numberOfStrips", "Tiff");
+#endif
+  }
+
+  DEFUN (__tiff_number_of_tiles__, args, ,
+             "Get the number of tiles in the image")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    if (nargin == 0)
+      error ("No handle provided\n");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+    
+    if (! TIFFIsTiled (tif))
+      error ("The image is stripped not tiled");
+    
+    double tile_count = static_cast<double> (TIFFNumberOfTiles (tif));
+    return octave_value_list (octave_value (tile_count));
+#else
+    err_disabled_feature ("numberOfTiles", "Tiff");
+#endif
+  }
+
+  DEFUN (__tiff_compute_strip__, args, ,
+             "Get the strip index containing the given row")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    if (nargin < 2 || nargin > 3)
+      error ("Wrong number of arguments\n");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+    
+    if (TIFFIsTiled (tif))
+      error ("The image is tiled not stripped");
+    
+    tiff_image_data image_data (tif);
+
+    uint32_t row = args(1).uint32_scalar_value ();
+    if (row > image_data.height)
+      row = image_data.height;
+
+    // Convert from 1-based to zero-based indexing but avoid underflow
+    if (row > 0)
+      row--;
+    
+    uint16_t plane;
+    if (nargin > 2)
+      {
+        if (image_data.planar_configuration == PLANARCONFIG_CONTIG)
+          error ("Can't use plane argument for images with chunky PlanarConfiguration");
+        plane = args(2).uint16_scalar_value ();
+        if (plane > image_data.samples_per_pixel)
+          plane = image_data.samples_per_pixel;
+        if (plane > 0)
+          plane--;
+      }
+    else
+      {
+        if (image_data.planar_configuration == PLANARCONFIG_SEPARATE)
+          error ("The plane argument is required for images with separate PlanarConfiguration");
+        plane = 0;
+      }
+
+    double strip_number = TIFFComputeStrip (tif, row, plane) + 1;
+    if (strip_number > TIFFNumberOfStrips (tif))
+      strip_number = TIFFNumberOfStrips (tif);
+    return octave_value_list (octave_value (strip_number));
+#else
+    err_disabled_feature ("computeStrip", "Tiff");
+#endif
+  }
+
+  DEFUN (__tiff_compute_tile__, args, ,
+             "Get the tile index containing the given row and column")
+  {
+#if defined (HAVE_TIFF)
+    int nargin = args.length ();
+
+    if (nargin < 2 || nargin > 3)
+      error ("Wrong number of arguments\n");
+    
+    TIFF *tif = (TIFF *)(args(0).uint64_value ());
+    
+    if (! TIFFIsTiled (tif))
+      error ("The image is stripped not tiled");
+    
+    uint32NDArray coords = args(1).uint32_array_value ();
+    if (coords.dim2() < 2)
+      error ("Coordinates must be in the shape [row, col]");
+    uint32_t row = coords(0, 0);
+    uint32_t col = coords(0, 1);
+
+    tiff_image_data image_data (tif);
+
+    if (col > image_data.width)
+      col = image_data.width;
+    if (row > image_data.height)
+      row = image_data.height;
+
+    // Convert from 1-based to zero-based indexing but avoid underflow
+    if (row > 0)
+      row--;
+    if (col > 0)
+      col--;
+    
+    uint16_t plane;
+    if (nargin > 2)
+      {
+        if (image_data.planar_configuration == PLANARCONFIG_CONTIG)
+          error ("Can't use plane argument for images with chunky PlanarConfiguration");
+        plane = args(2).uint16_scalar_value ();
+        if (plane > image_data.samples_per_pixel)
+          plane = image_data.samples_per_pixel;
+        if (plane > 0)
+          plane--;
+      }
+    else
+      {
+        if (image_data.planar_configuration == PLANARCONFIG_SEPARATE)
+          error ("The plane argument is required for images with separate PlanarConfiguration");
+        plane = 0;
+      }
+
+    double tile_number = TIFFComputeTile (tif, col, row, 0, plane) + 1;
+    if (tile_number > TIFFNumberOfTiles (tif))
+      tile_number = TIFFNumberOfTiles (tif);
+    return octave_value_list (octave_value (tile_number));
+#else
+    err_disabled_feature ("computeTile", "Tiff");
+#endif
+  }
+
+  DEFUN (__tiff_version__, , ,
+             "Get the version stamp of LibTIFF")
+  {
+#if defined (HAVE_TIFF)
+    std::string version = std::string (TIFFGetVersion ());
+    return octave_value_list (octave_value (version));
+#else
+    err_disabled_feature ("getVersion", "Tiff");
+#endif
+  }
+
+}

--- a/libinterp/corefcn/module.mk	Mon Aug 08 00:06:19 2022 +0200
+++ b/libinterp/corefcn/module.mk	Wed Aug 10 00:58:12 2022 +0200
@@ -126,6 +126,7 @@
   %reldir%/__magick_read__.cc \
   %reldir%/__pchip_deriv__.cc \
   %reldir%/__qp__.cc \
+  %reldir%/__tiff__.cc \
   %reldir%/amd.cc \
   %reldir%/auto-shlib.cc \
   %reldir%/balance.cc \

--- a/libinterp/dldfcn/__tiff__.cc	Mon Aug 08 00:06:19 2022 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,2368 +0,0 @@
-#if defined (HAVE_CONFIG_H)
-#  include "config.h"
-#endif
-
-#include <string>
-#include <iostream>
-
-// Workaround to not have to include fcntl which doesn't exist on windows
-// TODO(maged): see what octave does for this (fcntl)
-#ifndef O_RDONLY
-#define O_RDONLY 0
-#endif
-
-#include "defun-dld.h"
-#include "ov.h"
-#include "ovl.h"
-#include "error.h"
-
-#include "errwarn.h"
-
-#if defined (HAVE_TIFF)
-#  include <tiffio.h>
-#endif
-
-namespace octve
-{
-#if defined (HAVE_TIFF)
-
-  struct tiff_image_data
-  {
-  public:
-    uint32_t width;
-    uint32_t height;
-    uint16_t samples_per_pixel;
-    uint16_t bits_per_sample;
-    uint16_t planar_configuration;
-    uint16_t is_tiled;
-
-    tiff_image_data (TIFF *tif)
-    {
-      if (! TIFFGetField (tif, TIFFTAG_IMAGEWIDTH, &width))
-        error ("Failed to read image width");
-
-      if (! TIFFGetField (tif, TIFFTAG_IMAGELENGTH, &height))
-        error ("Failed to read image height");
-      
-      if (! TIFFGetFieldDefaulted (tif, TIFFTAG_SAMPLESPERPIXEL,
-                                   &samples_per_pixel))
-        error ("Failed to read the SamplesPerPixel tag");
-
-      if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
-                                   &bits_per_sample))
-        error ("Failed to read the BitsPerSample tag");
-      
-      // TODO(maged): this doesn't really work for writing as LibTIFF will
-      // refuse to write unless the tag is set
-      if (! TIFFGetField (tif, TIFFTAG_PLANARCONFIG,
-                          &planar_configuration))
-        // LibTIFF has a bug where it incorrectly returns 0 as a default
-        // value for PlanarConfiguration although the default value is 1
-        planar_configuration = 1;
-      
-      is_tiled = TIFFIsTiled(tif);
-    }
-  };
-
-  // Error if status is not 1 (success status for TIFFGetField)
-  void
-  validate_tiff_get_field (bool status, void *p_to_free=NULL)
-  {
-      if (status != 1)
-        {
-          if (p_to_free != NULL)
-            _TIFFfree (p_to_free);
-          error ("Failed to read tag");
-        }
-  }
-
-  uint32_t get_rows_in_strip (uint32_t strip_no, uint32_t strip_count,
-                              uint32_t rows_per_strip,
-                              tiff_image_data *image_data)
-  {
-    // Calculate the expected number of elements in the strip data array
-    // All strips have equal number of rows except strips at the bottom
-    // of the image can have less number of rows
-    uint32_t rows_in_strip = rows_per_strip;
-    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
-      {
-        // All strips have equal number of rows excpet strips at the bottom
-        // of the image can have less number of rows
-        if (strip_no == strip_count - 1)
-          rows_in_strip = image_data->height - rows_in_strip * strip_no;
-      }
-    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
-      {
-        // For separate planes, we should check the last strip of each channel
-        uint32_t strips_per_channel
-          = strip_count / image_data->samples_per_pixel;
-        for (uint32_t boundary_strip = strips_per_channel - 1;
-             boundary_strip <= strip_count;
-             boundary_strip += strips_per_channel)
-          if (strip_no == boundary_strip)
-            rows_in_strip = image_data->height
-                            - rows_in_strip * (strip_no % strips_per_channel);
-      }
-    else
-      error ("Planar Configuration not supported");
-    
-    return rows_in_strip;
-  }
-
-  template <typename T>
-  octave_value
-  read_strip (TIFF *tif, uint32_t strip_no, tiff_image_data * image_data)
-  {
-    // ASSUMES stripped image and strip_no is a valid zero-based strip
-    // index for the tif image
-    
-    uint32_t rows_in_strip;
-    if (! TIFFGetFieldDefaulted (tif, TIFFTAG_ROWSPERSTRIP, &rows_in_strip))
-      error ("Failed to obtain a value for RowsPerStrip");
-    
-    if (rows_in_strip > image_data->height)
-      rows_in_strip = image_data->height;
-    
-    uint32_t strip_count = TIFFNumberOfStrips (tif);
-    rows_in_strip = get_rows_in_strip (strip_no, strip_count,
-                                       rows_in_strip, image_data);
-    dim_vector strip_dims;
-    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
-      strip_dims = dim_vector (image_data->samples_per_pixel,
-                               image_data->width, rows_in_strip);
-    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
-      strip_dims = dim_vector (image_data->width, rows_in_strip, 1);
-    else
-      error ("Unsupported bit depth");
-    
-    T strip_data (strip_dims);
-    uint8_t *strip_fvec
-      = reinterpret_cast<uint8_t *> (strip_data.fortran_vec ());
-
-    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)
-      {
-        if (TIFFReadEncodedStrip (tif, strip_no, strip_fvec, -1) == -1)
-          error ("Failed to read strip data");
-      }
-    else if (image_data->bits_per_sample == 1)
-      {
-        if (image_data->samples_per_pixel != 1)
-          error ("Bi-Level images must have one channel only");
-        
-        // Create a buffer to hold the packed strip data
-        // Unique pointers are faster than vectors for constant size buffers
-        std::unique_ptr<uint8_t []> strip_ptr
-          = std::make_unique<uint8_t []> (TIFFStripSize (tif));
-        uint8_t *strip_buf = strip_ptr.get ();
-        if (TIFFReadEncodedStrip (tif, strip_no, strip_buf, -1) == -1)
-          error ("Failed to read strip data");
-        
-        // According to the format specification, the row should be byte
-        // aligned so the number of bytes is rounded up to the nearest byte
-        uint32_t padded_width = (image_data->width + 7) / 8;
-        // Packing the pixel data into bits
-        for (uint32_t row = 0; row < rows_in_strip; row++)
-          {
-            for (uint32_t col = 0; col < image_data->width; col++)
-            {
-              uint8_t shift = 7 - col % 8;
-              strip_fvec[col] = (strip_buf[col / 8] >> shift) & 0x1;
-            }
-            strip_fvec += image_data->width;
-            strip_buf += padded_width;
-          }
-      }
-    else
-      error ("Unsupported bit depth");
-    
-    Array<octave_idx_type> perm (dim_vector (3, 1));
-    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
-      {
-        perm(0) = 2;
-        perm(1) = 1;
-        perm(2) = 0;
-      }
-    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
-      {
-        perm(0) = 1;
-        perm(1) = 0;
-        perm(2) = 2;
-      }
-    
-    strip_data = strip_data.permute (perm);
-    return octave_value (strip_data);
-  }
-
-  template <typename T>
-  octave_value
-  read_tile (TIFF *tif, uint32_t tile_no, tiff_image_data *image_data)
-  {
-    // ASSUMES tiled image and tile_no is a valid zero-based tile
-    // index for the tif image
-
-    // TODO(maged): refactor into a function?
-    uint32_t tile_width, tile_height;
-    if (! TIFFGetField (tif, TIFFTAG_TILELENGTH, &tile_height))
-      error ("Filed to read tile length");
-    
-    if (! TIFFGetField (tif, TIFFTAG_TILEWIDTH, &tile_width))
-      error ("Filed to read tile length");
-
-    if (tile_height == 0 || tile_height % 16 != 0
-        || tile_width == 0 || tile_width % 16 != 0)
-      error ("Tile dimesion tags are invalid");
-
-    dim_vector tile_dims;
-    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
-      {
-        tile_dims = dim_vector (image_data->samples_per_pixel, tile_width,
-                                tile_height);
-      }
-    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
-      {
-        tile_dims = dim_vector (tile_width, tile_height, 1);
-      }
-    else
-      error ("Unsupported planar configuration");
-    
-    T tile_data (tile_dims);
-    uint8_t *tile_fvec
-      = reinterpret_cast<uint8_t *> (tile_data.fortran_vec ());
-
-    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)
-        {      
-          if (TIFFReadEncodedTile (tif, tile_no, tile_fvec, -1) == -1)
-            error ("Failed to read tile data");
-        }
-    else if (image_data->bits_per_sample == 1)
-      {
-        if (image_data->samples_per_pixel != 1)
-          error ("Bi-Level images must have one channel only");
-        
-        // Create a buffer to hold the packed tile data
-        // Unique pointers are faster than vectors for constant size buffers
-        std::unique_ptr<uint8_t []> tile_ptr
-          = std::make_unique<uint8_t []> (TIFFTileSize (tif));
-        uint8_t *tile_buf = tile_ptr.get ();
-        if (TIFFReadEncodedTile (tif, tile_no, tile_buf, -1) == -1)
-            error ("Failed to read tile data");
-        
-        // unpack tile bits into output matrix cells
-        for (uint32_t row = 0; row < tile_height; row++)
-          {
-            for (uint32_t col = 0; col < tile_width; col++)
-              {
-                uint8_t shift = 7 - col % 8;
-                tile_fvec[col] = (tile_buf [col / 8] >> shift) & 0x1;
-              }  
-            tile_fvec += tile_width;
-            tile_buf += tile_width / 8;
-          }
-      }
-    else
-      error ("Unsupported bit depth");
-    
-    Array<octave_idx_type> perm (dim_vector (3, 1));
-    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
-      {
-        perm(0) = 2;
-        perm(1) = 1;
-        perm(2) = 0;
-      }
-    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
-      {
-        perm(0) = 1;
-        perm(1) = 0;
-        perm(2) = 2;
-      }
-    
-    tile_data = tile_data.permute (perm);
-
-    // Get the actual tile dimensions
-    uint32_t tiles_across = (image_data->width + tile_width - 1)
-                            / tile_width;
-    uint32_t tiles_down = (image_data->height + tile_height - 1)
-                          / tile_height;
-    uint32_t corrected_width = tile_width;
-    uint32_t corrected_height = tile_height;
-    if (tile_no % tiles_across == tiles_across - 1)
-      corrected_width = image_data->width - tile_width * (tiles_across - 1);
-    if ((tile_no / tiles_across) % tiles_down == tiles_down - 1)
-      corrected_height = image_data->height - tile_height * (tiles_down - 1);
-    
-    dim_vector corrected_dims;
-    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
-      corrected_dims = dim_vector (corrected_height, corrected_width,
-                                   image_data->samples_per_pixel);
-    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
-      corrected_dims = dim_vector (corrected_height, corrected_width);
-    
-    tile_data.resize (corrected_dims);
-
-    return octave_value (tile_data);
-  }
-
-  template <typename T>
-  octave_value
-  read_strip_or_tile (TIFF *tif, uint32_t strip_tile_no,
-                      tiff_image_data *image_data)
-  {
-    if (image_data->is_tiled)
-      return read_tile<T> (tif, strip_tile_no, image_data);
-    else
-      return read_strip<T> (tif, strip_tile_no, image_data);
-  }
-
-  octave_value
-  handle_read_strip_or_tile (TIFF *tif, uint32_t strip_tile_no)
-  {
-    // Obtain all necessary tags
-    tiff_image_data image_data (tif);
-
-    uint16_t sample_format;
-    if (! TIFFGetFieldDefaulted(tif, TIFFTAG_SAMPLEFORMAT, &sample_format))
-      error ("Failed to obtain a value for sample format");
-
-    if (sample_format == 3)
-      {
-        if (image_data.bits_per_sample != 32
-            && image_data.bits_per_sample != 64)
-          error ("Floating point images are only supported for bit depths of 32 and 64");
-      }
-    else if (sample_format != 1 && sample_format != 4)
-      error ("Unsupported sample format");
-    
-    switch (image_data.bits_per_sample)
-      {
-      case 1:
-        return read_strip_or_tile<boolNDArray> (tif, strip_tile_no,
-                                                &image_data);
-        break;
-      case 8:
-        return read_strip_or_tile<uint8NDArray> (tif, strip_tile_no,
-                                                 &image_data);
-        break;
-      case 16:
-        return read_strip_or_tile<uint16NDArray> (tif, strip_tile_no,
-                                                  &image_data);
-        break;
-      case 32:
-        if (sample_format == 3)
-          return read_strip_or_tile<FloatNDArray> (tif, strip_tile_no,
-                                                   &image_data);
-        else
-          return read_strip_or_tile<uint32NDArray> (tif, strip_tile_no,
-                                                    &image_data);
-        break;
-      case 64:
-        if (sample_format == 3)
-          return read_strip_or_tile<NDArray> (tif, strip_tile_no,
-                                              &image_data);
-        else
-          return read_strip_or_tile<uint64NDArray> (tif, strip_tile_no,
-                                                    &image_data);
-        break;
-      default:
-        error ("Unsupported bit depth");
-      }
-  }
-
-  template <typename T>
-  octave_value
-  read_stripped_image (TIFF *tif, tiff_image_data *image_data)
-  {
-    typedef typename T::element_type P;
-
-    // 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
-    // size and the padding is removed at the end.
-    uint32_t padded_width = image_data->width;
-    uint8_t remove_padding = 0;
-    if ((image_data->bits_per_sample == 1 || image_data->bits_per_sample == 4)
-        && padded_width % 8 != 0)
-      {
-        padded_width += (8 - padded_width % 8) % 8;
-        remove_padding = 1;
-      }
-    
-    // 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));
-    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
-      img = T (dim_vector (padded_width, image_data->height,
-                           image_data->samples_per_pixel));
-    else
-      error ("Unsupported Planar Configuration");
-
-    P *img_fvec = img.fortran_vec ();
-
-    // Obtain the necessary data for handling the strips
-    uint32_t strip_count = TIFFNumberOfStrips (tif);
-    
-    // Can't rely on StripByteCounts because in compressed images
-    // the byte count reflect the actual number of bytes stored
-    // in the file not the size of the uncompressed strip
-    int64_t strip_size;
-    uint64_t written_size = 0;
-    uint64_t image_size = padded_width * image_data->height
-                          * image_data->samples_per_pixel
-                          * sizeof (P);
-    for (uint32_t strip = 0; strip < strip_count; strip++)
-      {
-        // Read the strip data into the matrix directly
-        strip_size = TIFFReadEncodedStrip (tif, strip, img_fvec, -1);
-
-        // Check if the strip read failed.
-        if (strip_size == -1)
-          error ("Failed to read strip data");
-        
-        // Check if the size being read exceeds the bounds of the matrix
-        // In case of a corrupt image with more data than needed
-        // This is probably redundant as LibTIFF checks sizes internally
-        if (written_size + strip_size > image_size)
-          error ("Strip data is larger than the image size");
-
-        if (image_data->bits_per_sample == 1)
-          {
-            if (image_data->samples_per_pixel != 1)
-              error ("Bi-Level images must have one channel only");
-            
-            // The strip size is multiplied by 8 to reflect tha actual
-            // number of bytes written to the matrix since each byte
-            // in the original strip contains 8 pixels of data
-            strip_size *= 8;
-
-            // Checking bounds again with the new size
-            if (written_size + strip_size > image_size)
-              error ("Strip 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 = strip_size - 1; pixel >= 0; pixel--)
-              {
-                // TODO(maged): is it necessary to check FillOrder?
-                uint8_t bit_number = 7 - pixel % 8;
-                uint8_t * img_u8 = reinterpret_cast<uint8_t *> (img_fvec);
-                img_fvec[pixel] = (img_u8[pixel / 8] >> bit_number) & 0x01;
-              }
-          }
-        else if (image_data->bits_per_sample == 4)
-          {
-            if (image_data->samples_per_pixel != 1)
-              error ("4-bit images are only supported for grayscale");
-            
-            // Strip size is multplied by as each byte contains 2 pixels
-            // and each pixel is later expanded into its own byte
-            strip_size *= 2;
-
-            // Checking bounds again with the ne size
-            if (written_size + strip_size > image_size)
-              error ("Strip 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 = strip_size - 1; pixel >= 0; pixel--)
-              {
-                uint8_t shift = pixel % 2 == 0? 4: 0;
-                uint8_t * img_u8 = reinterpret_cast<uint8_t *> (img_fvec);
-                img_fvec[pixel] = (img_u8[pixel / 2] >> shift) & 0x0F;
-              }
-          }
-        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");
-
-        // Advance the pointer by the amount of bytes read
-        img_fvec 
-          = reinterpret_cast<P *> (reinterpret_cast <uint8_t *> (img_fvec)
-                                   + strip_size);
-        written_size += strip_size;
-      }
-
-    // The matrices are permutated back to the shape expected by Octave
-    // which is height x width x channels
-    Array<octave_idx_type> perm (dim_vector (3, 1));
-    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
-      {
-        perm(0) = 2;
-        perm(1) = 1;
-        perm(2) = 0;
-      }
-    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
-      {
-        perm(0) = 1;
-        perm(1) = 0;
-        perm(2) = 2;
-      }
-    
-    img = img.permute (perm);
-
-    if (remove_padding)
-      img.resize (dim_vector (image_data->height, image_data->width));
-
-    return octave_value (img);
-  }
-
-  template <typename T>
-  octave_value
-  read_tiled_image (TIFF *tif, tiff_image_data *image_data)
-  {
-    typedef typename T::element_type P;
-
-    // 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;
-
-    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;    
-    uint64_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 (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;
-                uint8_t * img_u8 = reinterpret_cast<uint8_t *> (img_fvec);
-                img_fvec[pixel]= (img_u8[pixel / 8] >> bit_number) & 0x01;
-              }
-          }
-        else if (image_data->bits_per_sample == 4)
-          {
-            if (image_data->samples_per_pixel != 1)
-              error ("4-bit images are only supported for grayscale");
-            
-            // 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;
-                uint8_t * img_u8 = reinterpret_cast<uint8_t *> (img_fvec);
-                img_fvec[pixel] = (img_u8[pixel / 2] >> shift) & 0x0F;
-              }
-          }
-        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
-          = reinterpret_cast<P *> (reinterpret_cast<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);
-
-        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);
-  }
-
-  template <typename T>
-  octave_value
-  read_image (TIFF *tif, tiff_image_data *image_data)
-  {
-    if (image_data->is_tiled)
-      return read_tiled_image<T> (tif, image_data);
-    else
-      return read_stripped_image<T> (tif, image_data);
-  }
-
-  // Convert tag value to double
-  octave_value
-  interpret_scalar_tag_data (void *data, TIFFDataType tag_datatype)
-  {
-    double retval;
-
-    switch (tag_datatype)
-      {
-      case TIFF_BYTE:
-      case TIFF_UNDEFINED:
-        {
-          retval = static_cast<double> (*(reinterpret_cast<uint8_t *> (data)));
-          break;
-        }
-      case TIFF_SHORT:
-        {
-          retval = static_cast<double> (*(reinterpret_cast<uint16_t *> (data)));
-          break;
-        }
-      case TIFF_LONG:
-        {
-          retval = static_cast<double> (*(reinterpret_cast<uint32_t *> (data)));
-          break;
-        }
-      case TIFF_LONG8:
-        {
-          retval = static_cast<double> (*(reinterpret_cast<uint64_t *> (data)));
-          break;
-        }
-      case TIFF_RATIONAL:
-        {
-          error ("TIFF_RATIONAL should have at least 2 elements but got only 1");
-          break;
-        }
-      case TIFF_SBYTE:
-        {
-          retval = static_cast<double> (*(reinterpret_cast<int8_t *> (data)));
-          break;
-        }
-      case TIFF_SSHORT:
-        {
-          retval = static_cast<double> (*(reinterpret_cast<int16_t *> (data)));
-          break;
-        }
-      case TIFF_SLONG:
-        {
-          retval = static_cast<double> (*(reinterpret_cast<int32_t *> (data)));
-          break;
-        }
-      case TIFF_SLONG8:
-        {
-          retval = static_cast<double> (*(reinterpret_cast<int64_t *> (data)));
-          break;
-        }
-      case TIFF_FLOAT:
-        {
-          retval = *(reinterpret_cast<float *> (data));
-          break;
-        }
-      case TIFF_DOUBLE:
-        {
-          retval = *(reinterpret_cast<double *> (data));
-          break;
-        }
-      case TIFF_SRATIONAL:
-        {
-          error ("TIFF_SRATIONAL should have at least 2 elements but got only 1");
-          break;
-        }
-      case TIFF_IFD:
-      case TIFF_IFD8:
-        error ("Unimplemented IFFD data type");
-        break;
-      default:
-        error ("Unsupported tag data type");
-      }
-
-    return octave_value (retval);
-  }
-
-  // Convert memory buffer into a suitable octave value
-  // depending on tag_datatype
-  octave_value
-  interpret_tag_data (void *data, uint32_t count, TIFFDataType tag_datatype)
-  {
-    octave_value retval;
-    // Apparently matlab converts scalar numerical values into double
-    // but doesn't do the same for arrays
-    if (count == 1 && tag_datatype != TIFF_ASCII)
-      {
-        retval = interpret_scalar_tag_data (data, tag_datatype);
-      }
-    else
-      {
-        dim_vector arr_dims (1, count);
-
-        switch (tag_datatype)
-          {
-          case TIFF_BYTE:
-          case TIFF_UNDEFINED:
-            {
-              uint8NDArray arr (arr_dims);
-              for (uint32_t i = 0; i < count; i++)
-                {
-                  arr(i) = (reinterpret_cast<uint8_t *> (data))[i];
-                }
-              retval = octave_value (arr);
-              break;
-            }
-          case TIFF_ASCII:
-            {
-              retval = octave_value (*(reinterpret_cast<char **> (data)));
-              break;
-            }
-          case TIFF_SHORT:
-            {
-              uint16NDArray arr (arr_dims);
-              for (uint32_t i = 0; i < count; i++)
-                {
-                  arr(i) = (reinterpret_cast<uint16_t *> (data))[i];
-                }
-              retval = octave_value (arr);
-              break;
-            }
-          case TIFF_LONG:
-            {
-              uint32NDArray arr (arr_dims);
-              for (uint32_t i = 0; i < count; i++)
-                {
-                  arr(i) = (reinterpret_cast<uint32_t *> (data))[i];
-                }
-              retval = octave_value (arr);
-              break;
-            }
-          case TIFF_LONG8:
-            {
-              uint64NDArray arr (arr_dims);
-              for (uint32_t i = 0; i < count; i++)
-                {
-                  arr(i) = (reinterpret_cast<uint64_t *> (data))[i];
-                }
-              retval = octave_value (arr);
-              break;
-            }
-          case TIFF_RATIONAL:
-            {
-              NDArray arr (arr_dims);
-              for (uint32_t i = 0; i < count; i+=2)
-                {
-                  arr(i / 2) = static_cast<float> ((reinterpret_cast<uint32_t *> (data))[i])
-                               / static_cast<float> ((reinterpret_cast<uint32_t *> (data))[i+1]);
-                }
-              retval = octave_value (arr);
-              break;
-            }
-          case TIFF_SBYTE:
-            {
-              int8NDArray arr (arr_dims);
-              for (uint32_t i = 0; i < count; i++)
-                {
-                  arr(i) = (reinterpret_cast<int8_t *> (data))[i];
-                }
-              retval = octave_value (arr);
-              break;
-            }
-          case TIFF_SSHORT:
-            {
-              int16NDArray arr (arr_dims);
-              for (uint32_t i = 0; i < count; i++)
-                {
-                  arr(i) = (reinterpret_cast<int16_t *> (data))[i];
-                }
-              retval = octave_value (arr);
-              break;
-            }
-          case TIFF_SLONG:
-            {
-              int32NDArray arr (arr_dims);
-              for (uint32_t i = 0; i < count; i++)
-                {
-                  arr(i) = (reinterpret_cast<int32_t *> (data))[i];
-                }
-              retval = octave_value (arr);
-              break;
-            }
-          case TIFF_SLONG8:
-            {
-              int64NDArray arr (arr_dims);
-              for (uint32_t i = 0; i < count; i++)
-                {
-                  arr(i) = (reinterpret_cast<int64_t *> (data))[i];
-                }
-              retval = octave_value (arr);
-              break;
-            }
-          case TIFF_FLOAT:
-            {
-              NDArray arr (arr_dims);
-              for (uint32_t i = 0; i < count; i++)
-                {
-                  arr(i) = (reinterpret_cast<float *> (data))[i];
-                }
-              retval = octave_value (arr);
-              break;
-            }
-          case TIFF_DOUBLE:
-            {
-              NDArray arr (arr_dims);
-              for (uint32_t i = 0; i < count; i++)
-              {
-                  arr(i) = (reinterpret_cast<double *> (data))[i];
-              }
-              retval = octave_value (arr);
-              break;
-            }
-          case TIFF_SRATIONAL:
-            {
-              NDArray arr (arr_dims);
-              for (uint32_t i = 0; i < count; i+=2)
-                {
-                  arr(i / 2) = static_cast<float> ((reinterpret_cast<int32_t *> (data))[i])
-                               / static_cast<float> ((reinterpret_cast<int32_t *> (data))[i+1]);
-                }
-              retval = octave_value (arr);
-              break;
-            }
-          case TIFF_IFD:
-          case TIFF_IFD8:
-            // TODO(maged): implement IFD datatype?
-            error ("Unimplemented IFD data type");
-            break;
-          default:
-            error ("Unsupported tag data type");
-          }
-      }
-
-    return retval;
-  }
-
-  octave_value
-  get_scalar_field_data (TIFF *tif, const TIFFField *fip)
-  {
-    uint32_t tag_id = TIFFFieldTag (fip);
-
-    // TIFFFieldReadCount returns VARIABLE for some scalar tags
-    // (e.g. Compression) But TIFFFieldPassCount seems consistent
-    // Since scalar tags are the last to be handled, any tag that
-    // require a count to be passed is an unsupported tag.
-    if (TIFFFieldPassCount (fip))
-      error ("Unsupported tag");
-    
-    int type_size = TIFFDataWidth (TIFFFieldDataType (fip));
-    
-    void *data = _TIFFmalloc (type_size);
-    if (tag_id == TIFFTAG_PLANARCONFIG)
-      {
-        // Workaround for a bug in LibTIFF where it incorrectly returns
-        // zero as the default value for PlanaConfiguration
-        if (! TIFFGetField(tif, tag_id, data))
-          *reinterpret_cast<uint16_t *> (data) = 1;
-      }
-    else
-      validate_tiff_get_field (TIFFGetFieldDefaulted (tif, tag_id, data), data);
-    
-    octave_value tag_data_ov = interpret_tag_data (data, 1,
-                                                   TIFFFieldDataType (fip));
-    _TIFFfree (data);
-
-    return tag_data_ov;
-  }
-
-  octave_value
-  get_array_field_data (TIFF *tif, const TIFFField *fip, uint32_t array_size)
-  {
-    void *data;
-    validate_tiff_get_field (TIFFGetField (tif, TIFFFieldTag (fip), &data));
-      
-    return interpret_tag_data (data, array_size, TIFFFieldDataType (fip));
-  }
-
-  octave_value
-  get_field_data (TIFF *tif, const TIFFField *fip)
-  {
-    octave_value tag_data_ov;
-    uint32_t tag_id = TIFFFieldTag (fip);
-
-    // TODO(maged): find/create images to test the special tags
-    switch (tag_id)
-      {
-      case TIFFTAG_STRIPBYTECOUNTS:
-      case TIFFTAG_STRIPOFFSETS:
-        tag_data_ov = get_array_field_data (tif, fip, 
-                                            TIFFNumberOfStrips (tif));
-        break;
-      case TIFFTAG_TILEBYTECOUNTS:
-      case TIFFTAG_TILEOFFSETS:
-        tag_data_ov = get_array_field_data (tif, fip, TIFFNumberOfTiles (tif));
-        break;
-      case TIFFTAG_YCBCRCOEFFICIENTS:
-        tag_data_ov = get_array_field_data (tif, fip, 3);
-        break;
-      case TIFFTAG_REFERENCEBLACKWHITE:
-        tag_data_ov = get_array_field_data (tif, fip, 6);
-        break;
-      case TIFFTAG_GRAYRESPONSECURVE:
-        {
-          uint16_t bits_per_sample;
-          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
-                                            &bits_per_sample))
-            error ("Failed to obtain the bit depth");
-          
-          tag_data_ov = get_array_field_data (tif, fip, 1<<bits_per_sample);
-          break;
-        }
-      case TIFFTAG_COLORMAP:
-        {
-          uint16_t bits_per_sample;
-          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
-                                       &bits_per_sample))
-            error ("Failed to obtain the bit depth");
-          
-          // According to the format specification, this field should
-          // be 8 or 16 only.
-          if (bits_per_sample > 16)
-            error ("Too high bit depth for a palette image");
-
-          uint32_t count = 1 << bits_per_sample;
-          uint16_t *red, *green, *blue;
-          validate_tiff_get_field (TIFFGetField (tif, TIFFTAG_COLORMAP,
-                                                 &red, &green, &blue));
-          
-          // Retrieving the data of the three channels and concatenating
-          // them together
-          OCTAVE_LOCAL_BUFFER (NDArray, array_list, 3);
-          dim_vector col_dims(count, 1);
-          array_list[0] = NDArray (interpret_tag_data (red,
-                                                       count,
-                                                       TIFFFieldDataType(fip))
-                                                       .uint16_array_value ()
-                                                       .reshape (col_dims));
-          array_list[1] = NDArray (interpret_tag_data (green,
-                                                       count,
-                                                       TIFFFieldDataType(fip))
-                                                       .uint16_array_value ()
-                                                       .reshape (col_dims));
-          array_list[2] = NDArray (interpret_tag_data (blue,
-                                                       count,
-                                                       TIFFFieldDataType(fip))
-                                                       .uint16_array_value ()
-                                                       .reshape (col_dims));
-
-          NDArray mat_out = NDArray::cat(1, 3, array_list);
-          // Normalize the range to be between 0 and 1
-          mat_out /= UINT16_MAX;
-
-          tag_data_ov = octave_value (mat_out);
-          break;
-        }
-      case TIFFTAG_TRANSFERFUNCTION:
-        {
-          uint16_t samples_per_pixel;
-          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_SAMPLESPERPIXEL,
-                                       &samples_per_pixel))
-            error ("Failed to obtain the number of samples per pixel");
-
-          uint16_t bits_per_sample;
-          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
-                                       &bits_per_sample))
-            error ("Failed to obtain the number of samples per pixel");
-
-          uint32_t count = 1 << bits_per_sample;
-          uint16_t *ch1, *ch2, *ch3;
-          if (samples_per_pixel == 1)
-            {
-              validate_tiff_get_field (TIFFGetField (tif, TIFFTAG_COLORMAP, &ch1));
-              tag_data_ov = interpret_tag_data (ch1, count,
-                                                TIFFFieldDataType (fip));
-            }
-          else
-            {
-              validate_tiff_get_field (TIFFGetField (tif, TIFFTAG_COLORMAP,
-                                                     &ch1, &ch2, &ch3));
-              
-              OCTAVE_LOCAL_BUFFER (uint16NDArray, array_list, 3);
-              dim_vector col_dims(count, 1);
-              array_list[0] = interpret_tag_data (ch1,
-                                                  count,
-                                                  TIFFFieldDataType (fip))
-                                                  .uint16_array_value ()
-                                                  .reshape (col_dims);
-              array_list[1] = interpret_tag_data (ch2,
-                                                  count,
-                                                  TIFFFieldDataType (fip))
-                                                  .uint16_array_value ()
-                                                  .reshape (col_dims);
-              array_list[2] = interpret_tag_data (ch3,
-                                                  count,
-                                                  TIFFFieldDataType (fip))
-                                                  .uint16_array_value ()
-                                                  .reshape (col_dims);
-              
-              tag_data_ov = octave_value (uint16NDArray::cat (1, 3, array_list));
-            }
-          break;
-        }
-      case TIFFTAG_PAGENUMBER:
-      case TIFFTAG_HALFTONEHINTS:
-      case TIFFTAG_DOTRANGE:
-      case TIFFTAG_YCBCRSUBSAMPLING:
-        {
-          uint16_t tag_part1, tag_part2;
-          validate_tiff_get_field (TIFFGetField (tif, tag_id,
-                                                 &tag_part1, &tag_part2));
-          
-          NDArray mat_out (dim_vector (1, 2));
-          mat_out(0)
-            = interpret_tag_data (&tag_part1, 1,
-                                  TIFFFieldDataType (fip)).double_value ();
-          mat_out(1)
-            = interpret_tag_data (&tag_part2, 1,
-                                  TIFFFieldDataType (fip)).double_value ();
-          
-          tag_data_ov = octave_value (mat_out);
-          break;
-        }
-      case TIFFTAG_SUBIFD:
-        {
-          uint16_t count;
-          uint64_t *offsets;
-          validate_tiff_get_field (TIFFGetField (tif, tag_id, &count, &offsets));
-          tag_data_ov = interpret_tag_data (offsets, count,
-                                            TIFFFieldDataType (fip));
-          break;
-        }
-      case TIFFTAG_EXTRASAMPLES:
-        {
-          uint16_t count;
-          uint16_t *types;
-          validate_tiff_get_field (TIFFGetField (tif, tag_id, &count, &types));
-          tag_data_ov = interpret_tag_data (types, count,
-                                            TIFFFieldDataType (fip));
-          break;
-        }
-      // TODO(maged): These tags are more complex to implement
-      //  will be implemented and tested later.
-      case TIFFTAG_XMLPACKET:
-      case TIFFTAG_RICHTIFFIPTC:
-      case TIFFTAG_PHOTOSHOP:
-      case TIFFTAG_ICCPROFILE:
-        {
-          error ("Complex Tags not implemented");
-          break;
-        }
-      // These tags are not mentioned in the LibTIFF documentation
-      // but are handled correctly by the library
-      case TIFFTAG_ZIPQUALITY:
-      case TIFFTAG_SGILOGDATAFMT:
-      case TIFFTAG_GRAYRESPONSEUNIT:
-        {
-          tag_data_ov = get_scalar_field_data (tif, fip);
-          break;
-        }
-      default:
-        tag_data_ov = get_scalar_field_data (tif, fip);
-    }
-    
-    return tag_data_ov;
-  }
-
-  void
-  set_field_data (TIFF *tif, const TIFFField *fip, octave_value tag_ov)
-  {
-    // TODO(maged): prevent calling here for read-only images
-    
-    // TODO(maged): complete the implementation of this function
-    uint32_t tag_id = TIFFFieldTag (fip);
-    uint32_t tag_data = tag_ov.double_value ();
-
-    if (! TIFFSetField(tif, tag_id, tag_data))
-      error ("Failed to set tag value");
-  }
-  
-  template <typename T>
-  void
-  write_strip (TIFF *tif, uint32_t strip_no, T strip_data,
-               tiff_image_data *image_data)
-  { 
-    uint32_t rows_in_strip;
-    if (! TIFFGetFieldDefaulted (tif, TIFFTAG_ROWSPERSTRIP, &rows_in_strip))
-      error ("Failed to obtain the RowsPerStrip tag");
-    
-    // The tag has a default value of UINT32_MAX which means the entire
-    // image, but we need to cap it to the height for later calculations
-    if (rows_in_strip > image_data->height)
-      rows_in_strip = image_data->height;
-    
-    // LibTIFF uses zero-based indexing as opposed to Octave's 1-based
-    strip_no--;
-
-    uint32_t strip_count = TIFFNumberOfStrips (tif);
-    rows_in_strip = get_rows_in_strip (strip_no, strip_count,
-                                       rows_in_strip, image_data);
-
-    dim_vector strip_dimensions;
-    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
-      strip_dimensions = dim_vector (rows_in_strip, image_data->width,
-                                     image_data->samples_per_pixel);
-    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
-      strip_dimensions = dim_vector (rows_in_strip, image_data->width, 1);
-    else
-      error ("Planar configuration not supported");
-
-    if (strip_data.dim1 () > rows_in_strip)
-      warning ("The strip is composed of %u rows but the input has %ld rows.",
-               rows_in_strip,
-               strip_data.dim1 ());
-    
-    if (strip_data.dim2 () > image_data->width)
-      warning ("The image width is %u but the input has %ld columns.",
-               image_data->width,
-               strip_data.dim2 ());
-    
-    if (strip_data.ndims () > 2)
-      {
-        if (image_data->planar_configuration == PLANARCONFIG_CONTIG
-            && strip_data.dim3 () > image_data->samples_per_pixel)
-          warning ("The strip is composed of %u channels but the input has %ld channels.",
-                   image_data->samples_per_pixel,
-                   strip_data.dim3 ());
-        else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE
-                && strip_data.dim3 () > 1)
-          warning ("The strip is composed of %u channel but the input has %ld channels.",
-                   1, strip_data.dim3 ());
-      }
-
-    strip_data.resize (strip_dimensions);
-
-    // Permute the dimesions of the strip to match the expected memory
-    // arrangement of LibTIFF (channel x width x height)
-    Array<octave_idx_type> perm (dim_vector (3, 1));
-    perm(0) = 2;
-    perm(1) = 1;
-    perm(2) = 0;
-    strip_data = strip_data.permute (perm);
-
-    void *data_vec = strip_data.fortran_vec ();
-    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)
-      {
-        // Can't rely on LibTIFF's TIFFStripSize because boundary strips
-        // can be smaller in size
-        tsize_t strip_size = strip_data.numel ()
-                             * image_data->bits_per_sample / 8;
-        if (TIFFWriteEncodedStrip (tif, strip_no, data_vec, strip_size) == -1)
-          error ("Failed to write strip data to image");
-        
-      }
-    else if (image_data->bits_per_sample == 1)
-      {
-        if (image_data->samples_per_pixel != 1)
-          error ("Bi-Level images must have one channel only");
-        
-        // Create a buffer to hold the packed strip data
-        // Unique pointers are faster than vectors for constant size buffers
-        std::unique_ptr<uint8_t []> strip_ptr
-          = std::make_unique<uint8_t []> (TIFFStripSize (tif));
-        uint8_t *strip_buf = strip_ptr.get ();
-        uint8_t *data_u8 = reinterpret_cast<uint8_t *> (data_vec);
-        // According to the format specification, the row should be byte
-        // aligned so the number of bytes is rounded up to the nearest byte
-        uint32_t padded_width = (image_data->width + 7) / 8;
-        // Packing the pixel data into bits
-        for (uint32_t row = 0; row < rows_in_strip; row++)
-          {
-            for (uint32_t col = 0; col < image_data->width; col++)
-            {
-              uint8_t shift = 7 - col % 8;
-              strip_buf[row * padded_width + col / 8] |= data_u8[col] << shift;
-            }
-            data_u8 += image_data->width;
-          }
-        tsize_t strip_size = padded_width * rows_in_strip;
-        if (TIFFWriteEncodedStrip (tif, strip_no, strip_buf, strip_size) == -1)
-          error ("Failed to write strip data to image");
-      }
-    else
-      {
-        error ("Unsupported bit depth");
-      }
-  }
-
-  template <typename T>
-  void
-  write_tile (TIFF *tif, uint32_t tile_no, T tile_data,
-              tiff_image_data *image_data)
-  {
-    uint32_t tile_width, tile_height;
-    if (! TIFFGetField (tif, TIFFTAG_TILEWIDTH, &tile_width))
-      error ("Failed to get the tile width");
-    if (! TIFFGetField (tif, TIFFTAG_TILELENGTH, &tile_height))
-      error ("Failed to get the tile length");
-
-    if (tile_height == 0 || tile_height % 16 != 0
-        || tile_width == 0 || tile_width % 16 != 0)
-      error ("Tile dimesion tags are invalid");
-    
-    if (tile_no < 1 || tile_no > TIFFNumberOfTiles (tif))
-      error ("Tile number out of bounds");
-
-    if (tile_data.dim1 () > tile_height)
-      warning ("The tile is composed of %u rows but input has %ld rows",
-               tile_height, tile_data.dim1 ());
-    if (tile_data.dim2 () > tile_width)
-      warning ("The tile is composed of %u columns but input has %ld columns",
-               tile_width, tile_data.dim2 ());
-    if (tile_data.ndims () > 2)
-      {
-        if (image_data->planar_configuration == PLANARCONFIG_CONTIG
-            && tile_data.dim3 () > image_data->samples_per_pixel)
-          warning ("The tile is composed of %u channels but input has %ld channels",
-                  image_data->samples_per_pixel, tile_data.dim3 ());
-        else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE
-                 && tile_data.dim3 () > 1)
-          warning ("The tile is composed of %u channels but input has %ld channels",
-                   1, tile_data.dim3 ());
-      }
-    
-    dim_vector tile_dimensions;
-    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
-      tile_dimensions = dim_vector (tile_height, tile_width,
-                                    image_data->samples_per_pixel);
-    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
-      tile_dimensions = dim_vector (tile_height, tile_width, 1);
-    else
-      error ("Planar configuration not supported");
-    
-    tile_data.resize (tile_dimensions);
-    Array<octave_idx_type> perm (dim_vector (3, 1));
-    perm(0) = 2;
-    perm(1) = 1;
-    perm(2) = 0;
-    tile_data = tile_data.permute (perm);
-    
-    // Octave indexing is 1-based while LibTIFF is zero-based
-    tile_no--;
-    void *data_vec = tile_data.fortran_vec ();
-    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)
-      {
-        if (TIFFWriteEncodedTile (tif, tile_no, data_vec,
-                                  TIFFTileSize (tif)) == -1)
-          error ("Failed to write tile data to image");
-        
-      }
-    else if (image_data->bits_per_sample == 1)
-      {
-        if (image_data->samples_per_pixel != 1)
-          error ("Bi-Level images must have one channel only");
-        
-        // Create a buffer to hold the packed tile data
-        // Unique pointers are faster than vectors for constant size buffers
-        std::unique_ptr<uint8_t []> tile_ptr
-          = std::make_unique<uint8_t []> (TIFFTileSize (tif));
-        uint8_t *tile_buf = tile_ptr.get ();
-        uint8_t *data_u8 = reinterpret_cast<uint8_t *> (data_vec);
-        // Packing the pixel data into bits
-        for (uint32_t row = 0; row < tile_height; row++)
-          {
-            for (uint32_t col = 0; col < tile_width; col++)
-            {
-              uint8_t shift = 7 - col % 8;
-              tile_buf[row * tile_width/8 + col/8] |= data_u8[col] << shift;
-            }
-            data_u8 += tile_width;
-          }
-        if (TIFFWriteEncodedTile (tif, tile_no, tile_buf,
-                                  TIFFTileSize (tif)) == -1)
-          error ("Failed to write tile data to image");
-      }
-    else
-      {
-        error ("Unsupported bit depth");
-      }
-  }
-
-  template <typename T>
-  void
-  write_strip_or_tile (TIFF *tif, uint32_t strip_tile_no, T strip_data,
-                       tiff_image_data *image_data)
-  {
-    if (image_data->is_tiled)
-      write_tile<T> (tif, strip_tile_no, strip_data, image_data);
-    else
-      write_strip<T> (tif, strip_tile_no, strip_data, image_data);
-  }
-
-  void
-  handle_write_strip_or_tile (TIFF *tif, uint32_t strip_tile_no,
-                              octave_value data_ov,
-                              tiff_image_data *image_data)
-  {
-
-    // SampleFormat tag is not a required field and has a default value of 1
-    // So we need to use TIFFGetFieldDefaulted in case it is not present in
-    // the file
-    uint16_t sample_format;
-    if (! TIFFGetFieldDefaulted(tif, TIFFTAG_SAMPLEFORMAT, &sample_format))
-      error ("Failed to obtain a value for sample format");
-
-    // TODO(maged): add support for signed integer images
-    if (sample_format == 3)
-      {
-        if (image_data->bits_per_sample != 32
-            && image_data->bits_per_sample != 64)
-          error ("Floating point images are only supported for bit depths of 32 and 64");
-      }
-
-    // The standard specifies that a SampleFormat of 4 should be treated
-    // the same as 1 (unsigned integer)
-    else if (sample_format != 1 && sample_format != 4)
-      error ("Unsupported sample format");
-    
-    switch (image_data->bits_per_sample)
-      {
-      case 1:
-        // We need to check for both scalar and matrix types to handle single
-        // element strip
-        if (data_ov.is_bool_scalar () || data_ov.is_bool_matrix ())
-          write_strip_or_tile<boolNDArray> (tif, strip_tile_no,
-                                            data_ov.bool_array_value (),
-                                            image_data);
-        else
-          error ("Expected logical matrix for BiLevel image");
-        break;
-      case 8:
-        if (data_ov.is_uint8_type ())
-          write_strip_or_tile<uint8NDArray> (tif, strip_tile_no,
-                                             data_ov.uint8_array_value (),
-                                             image_data);
-        else
-          error ("Only uint8 data is allowed for uint images with bit depth of 8");
-        break;
-      case 16:
-        if (data_ov.is_uint16_type ())
-          write_strip_or_tile<uint16NDArray> (tif, strip_tile_no,
-                                              data_ov.uint16_array_value (),
-                                              image_data);
-        else
-          error ("Only uint16 data is allowed for uint images with bit depth of 16");
-        break;
-      case 32:
-        if (sample_format == 3)
-          if (data_ov.is_single_type () || data_ov.is_double_type ())
-            write_strip_or_tile<FloatNDArray> (tif, strip_tile_no,
-                                               data_ov.float_array_value (),
-                                               image_data);
-          else
-            error ("Only single and double data are allowed for floating-point images");
-        else
-          if (data_ov.is_uint32_type ())
-            write_strip_or_tile<uint32NDArray> (tif, strip_tile_no,
-                                                data_ov.uint32_array_value (),
-                                                image_data);
-          else
-            error ("Only uint32 data is allowed for uint images with bit depth of 32");
-        break;
-      case 64:
-        if (sample_format == 3)
-          if (data_ov.is_single_type () || data_ov.is_double_type ())
-            write_strip_or_tile<NDArray> (tif, strip_tile_no,
-                                          data_ov.array_value (),
-                                          image_data);
-          else
-            error ("Only single and double data are allowed for floating-point images");
-        else  
-          if (data_ov.is_uint64_type ())
-            write_strip_or_tile<uint64NDArray> (tif, strip_tile_no,
-                                                data_ov.uint64_array_value (),
-                                                image_data);
-          else
-            error ("Only uint64 data is allowed for uint images with bit depth of 64");
-        break;
-      default:
-        error ("Unsupported bit depth");
-      }
-  }
-
-  template <typename T>
-  void
-  write_stripped_image (TIFF *tif, T pixel_data, tiff_image_data *image_data)
-  {
-    // TODO(maged): remove this? ASSUMES pixel data dimensions are already validated
-
-    typedef typename T::element_type P;
-
-    // Permute pixel data to be aligned in memory to the way LibTIFF
-    // expects the data to be (i.e. channel x width x height for chunky
-    // and width x height x channel for separate planes)
-    Array<octave_idx_type> perm (dim_vector (3, 1));
-    if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
-      {
-        perm(0) = 1;
-        perm(1) = 0;
-        perm(2) = 2;
-      }
-    else
-      {
-        perm(0) = 2;
-        perm(1) = 1;
-        perm(2) = 0;
-      }
-    pixel_data = pixel_data.permute (perm);
-
-    uint32_t row_per_strip;
-    if (! TIFFGetFieldDefaulted (tif, TIFFTAG_ROWSPERSTRIP, &row_per_strip))
-      error ("Failed to obtain the RowPerStrip tag");
-    
-    // The default value is INT_MAX so we need to cap it to the image height
-    if (row_per_strip > image_data->height)
-      row_per_strip = image_data->height;
-
-    uint8_t *pixel_fvec = reinterpret_cast<uint8_t *> (pixel_data.fortran_vec ());
-    uint32_t strip_count = TIFFNumberOfStrips (tif);
-    tsize_t strip_size;
-    uint32_t rows_in_strip;
-    for (uint32_t strip = 0; strip < strip_count; strip++)
-      {
-        rows_in_strip = get_rows_in_strip (strip, strip_count,
-                                           row_per_strip, image_data);
-        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)
-          {
-            strip_size = rows_in_strip * image_data->width * sizeof (P);
-            if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
-              strip_size *= image_data->samples_per_pixel;
-            if (! TIFFWriteEncodedStrip (tif, strip, pixel_fvec, strip_size))
-              error ("Failed to rite strip data");
-            pixel_fvec += strip_size;
-          }
-        else if (image_data->bits_per_sample == 1)
-          {
-            if (image_data->samples_per_pixel != 1)
-              error ("Bi-Level images must have one channel only");
-            
-            // Create a buffer to hold the packed strip data
-            // Unique pointers are faster than vectors for constant size buffers
-            std::unique_ptr<uint8_t []> strip_ptr
-              = std::make_unique<uint8_t []> (TIFFStripSize (tif));
-            uint8_t *strip_buf = strip_ptr.get ();
-            // According to the format specification, the row should be byte
-            // aligned so the number of bytes is rounded up to the nearest byte
-            uint32_t padded_width = (image_data->width + 7) / 8;
-            // Packing the pixel data into bits
-            for (uint32_t row = 0; row < rows_in_strip; row++)
-              {
-                for (uint32_t col = 0; col < image_data->width; col++)
-                {
-                  uint8_t shift = 7 - col % 8;
-                  strip_buf[row * padded_width + col / 8]
-                    |= pixel_fvec[col] << shift;
-                }
-                pixel_fvec += image_data->width;
-              }
-            strip_size = padded_width * rows_in_strip;
-            if (TIFFWriteEncodedStrip (tif, strip, strip_buf, strip_size) == -1)
-              error ("Failed to write strip data to image");
-          }
-        else
-          error ("Unsupported bit depth");
-      }
-  }
-
-  template <typename T>
-  void
-  write_tiled_image (TIFF *tif, T pixel_data, tiff_image_data *image_data)
-  {
-    // TODO(maged): remove this? ASSUMES pixel data dimensions are already validated
-
-    uint32_t tile_width, tile_height;
-    if (! TIFFGetField (tif, TIFFTAG_TILEWIDTH, &tile_width))
-      error ("Failed to get the tile width");
-    if (! TIFFGetField (tif, TIFFTAG_TILELENGTH, &tile_height))
-      error ("Failed to get the tile length");
-
-    if (tile_height == 0 || tile_height % 16 != 0
-        || tile_width == 0 || tile_width % 16 != 0)
-      error ("Tile dimesion tags are invalid");
-    
-    uint32_t tiles_across = (image_data->width + tile_width - 1)
-                            / tile_width;
-    uint32_t tiles_down = (image_data->height + tile_height - 1)
-                          / tile_height;
-    
-    // Resize the image data to add tile padding
-    dim_vector padded_dims (tiles_down * tile_height,
-                            tiles_across * tile_width,
-                            image_data->samples_per_pixel);
-    pixel_data.resize (padded_dims);
-
-    // Reshape the image to separate tiles into their own dimension
-    // so we can permute them to the right order
-    dim_vector tiled_dims (tile_height, tiles_down, tile_width, tiles_across,
-                           image_data->samples_per_pixel);
-    pixel_data = pixel_data.reshape (tiled_dims);
-
-    // Permute the dimesnions to get the memory alignment to match LibTIFF
-    Array<octave_idx_type> perm (dim_vector (5, 1));
-    if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
-      {
-        // For separate planes, the data coming from libtiff will have all
-        // tiles of the first sample then all tiles of the second sample
-        // and so on. Tiles of each sample will be ordered from left to right
-        // and from top to bottom. And data inside each tile is organised as
-        // rows and each row contains columns.
-        // So the order for LibTIFF is:
-        //   samples x tiles_down x tiles_across x tile_height x tile_width
-        // But memory orientation of Octave arrays is reversed so we set it to
-        //   tile_width x tile_height x tiles_across x tiles_down x samples
-        perm(0) = 2;
-        perm(1) = 0;
-        perm(2) = 3;
-        perm(3) = 1;
-        perm(4) = 4;
-      }
-    else
-      {
-        // For chunky planes, the data coming from libtiff will be ordered
-        // from left to right and from top to bottom. And data inside each
-        // tile is organised as rows and each row contains columns and each
-        // column contains samples.
-        // So the order for LibTIFF is:
-        //   tiles_down x tiles_across x tile_height x tile_width x samples
-        // But memory orientation of Octave arrays is reversed so we set it to
-        //   samples x tile_width x tile_height x tiles_across x tiles_down
-        perm(0) = 4;
-        perm(1) = 2;
-        perm(2) = 0;
-        perm(3) = 3;
-        perm(4) = 1;
-      }
-    pixel_data = pixel_data.permute (perm);
-
-    uint8_t *pixel_fvec
-      = reinterpret_cast<uint8_t *> (pixel_data.fortran_vec ());
-    uint32_t tile_count = TIFFNumberOfTiles (tif);
-    tsize_t tile_size = TIFFTileSize (tif);
-
-    for (uint32_t tile = 0; tile <tile_count; tile++)
-      {
-        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)
-          {
-            if (! TIFFWriteEncodedTile (tif, tile, pixel_fvec, tile_size))
-              error ("Failed to write tile data");
-            pixel_fvec += tile_size;
-          }
-        else if (image_data->bits_per_sample == 1)
-          {
-            if (image_data->samples_per_pixel != 1)
-              error ("Bi-Level images must have one channel only");
-            
-            // Create a buffer to hold the packed tile data
-            // Unique pointers are faster than vectors for
-            // constant size buffers
-            std::unique_ptr<uint8_t []> tile_ptr
-              = std::make_unique<uint8_t []> (tile_size);
-            uint8_t *tile_buf = tile_ptr.get ();
-            // Packing the pixel data into bits
-            for (uint32_t row = 0; row < tile_height; row++)
-              {
-                for (uint32_t col = 0; col < tile_width; col++)
-                {
-                  uint8_t shift = 7 - col % 8;
-                  tile_buf[row * tile_width / 8 + col / 8]
-                    |= pixel_fvec[col] << shift;
-                }
-                pixel_fvec += tile_width;
-              }
-            if (TIFFWriteEncodedTile (tif, tile, tile_buf,
-                                      TIFFTileSize (tif)) == -1)
-              error ("Failed to write tile data to image");
-          }
-        else
-          error ("Unsupported bit depth");
-      }
-
-  }
-
-  template <typename T>
-  void
-  write_image (TIFF *tif, T pixel_data, tiff_image_data *image_data)
-  {
-    // TODO(maged): matlab sets the remaining to zero for less width and height
-    // and issues a warning for larger widths but not for larger height?
-    // and errors for less or more number of channels
-    if (image_data->height != pixel_data.dim1 ()
-        || image_data->width != pixel_data.dim2 ()
-        || pixel_data.ndims () > 3
-        || (image_data->samples_per_pixel > 1 && pixel_data.ndims () < 3)
-        || (pixel_data.ndims () > 2
-            && image_data->samples_per_pixel != pixel_data.dim3 ()))
-      error ("Dimensions of the input don't match image dimenions");
-    
-    if (image_data->is_tiled)
-      write_tiled_image<T> (tif, pixel_data, image_data);
-    else
-      write_stripped_image<T> (tif, pixel_data, image_data);
-
-  }
-
-
-#endif
-
-  DEFUN_DLD (__open_tiff__, args, ,
-             "Open a Tiff file and return its handle")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    if (nargin == 0 || nargin > 2)
-      {
-        // TODO(maged): return invalid object instead??
-        error ("No filename supplied\n");
-      }
-
-    std::string filename = args(0).string_value ();
-    std::string mode = "r";
-
-    if (nargin == 2)
-      mode = args(1).string_value ();
-
-    const std::vector<std::string> supported_modes {"r", "w", "w8", "a"};
-      
-    if (std::find (supported_modes.cbegin (), supported_modes.cend (), mode)
-          == supported_modes.cend ())
-      {
-        if (mode == "r+")
-          error ("Openning files in r+ mode is not yet supported");
-        else
-          error ("Invalid mode for openning Tiff file: %s", mode.c_str ());
-      }
-    
-    TIFF *tif = TIFFOpen (filename.c_str (), mode.c_str ());
-    
-    if (! tif)
-      error ("Failed to open Tiff file\n");
-
-    // TODO(maged): use inheritance of octave_base_value instead
-    octave_value tiff_ov = octave_value ((uint64_t)tif);
-    return octave_value_list (tiff_ov);
-#else
-    octave_unused_parameter (args);
-    err_disabled_feature ("Tiff", "Tiff");
-#endif
-  }
-
-
-  DEFUN_DLD (__close_tiff__, args, ,
-             "Close a tiff file")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    if (nargin == 0)
-      error ("No handle provided\n");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-    TIFFClose (tif);
-
-    return octave_value_list ();
-#else
-    err_disabled_feature ("close", "Tiff");
-#endif
-  }
-
-
-  DEFUN_DLD (__tiff_get_tag__, args, ,
-             "Get the value of a tag from a tiff image")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    if (nargin == 0)
-      error ("No handle provided\n");
-    
-    if (nargin < 2)
-      error ("No tag name provided\n");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-
-    uint32_t tag_id;
-    const TIFFField *fip;
-    
-    if (args(1).is_string ())
-      {
-        std::string tagName = args(1).string_value ();
-        fip = TIFFFieldWithName (tif, tagName.c_str ());
-        if (! fip)
-          error ("Tiff tag not found");
-        
-        tag_id = TIFFFieldTag (fip);
-      }
-    else
-      {
-        tag_id = args(1).int_value ();
-        fip = TIFFFieldWithTag (tif, tag_id);
-        
-        if (! fip)
-          error ("Tiff tag not found");
-      }
-
-    return octave_value_list (get_field_data (tif, fip));
-#else
-    err_disabled_feature ("getTag", "Tiff");
-#endif
-  }
-
-
-  DEFUN_DLD (__tiff_set_tag__, args, ,
-             "Set the value of a tag in a tiff image")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-    
-    if (nargin < 2)
-      error ("Too few arguments provided\n");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-    
-    if (args(1).isstruct ())
-      {
-        octave_scalar_map tags = args(1).xscalar_map_value ("__tiff_set_tag__: struct argument must be a scalar structure");
-        string_vector keys = tags.fieldnames ();
-        // Using iterators instead of this loop method seems to process
-        // the elements of the struct in a different order than they were
-        // assigned which makes the behavior here incompatible with matlab
-        // in case of errors.
-        for (octave_idx_type i = 0; i < keys.numel (); i++)
-          {
-            std::string key = keys[i];
-            const TIFFField *fip =  TIFFFieldWithName (tif, key.c_str ());
-            if (! fip)
-              error ("Tag %s not found", key.c_str ());
-            set_field_data (tif, fip, tags.contents (key));
-          }
-      }
-    else
-      {
-        if (nargin < 3)
-          error ("Too few arguments provided");
-
-        const TIFFField *fip;
-        // TODO(maged): matlab actually checks for its own strings not LibTIFF's
-        if (args(1).is_string ())
-          {
-            std::string tagName = args(1).string_value ();
-            fip = TIFFFieldWithName (tif, tagName.c_str ());
-            if (! fip)
-              error ("Tiff tag not found");
-          }
-        else
-          {
-            uint32_t tag_id = args(1).int_value ();
-            fip = TIFFFieldWithTag (tif, tag_id);
-            if (! fip)
-              error ("Tiff tag not found");
-          }
-        
-        set_field_data (tif, fip, args(2));
-      }
-
-    return octave_value_list ();
-#else
-    err_disabled_feature ("setTag", "Tiff");
-#endif
-  }
-
-  DEFUN_DLD (__tiff_read__, args, nargout,
-             "Read the image in the current IFD")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    if (nargin == 0)
-      error ("No handle provided\n");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-
-    // TODO(maged): nargout and ycbcr
-    octave_unused_parameter (nargout);
-
-    // Obtain all necessary tags
-    // The main purpose of this struct is to hold all the necessary tags that
-    // will be common among all the different cases of handling the the image
-    // data to avoid repeating the same calls to TIFFGetField in the different
-    // functions as each call is a possible point of failure
-    tiff_image_data image_data (tif);
-
-    uint16_t sample_format;
-    if (! TIFFGetFieldDefaulted(tif, TIFFTAG_SAMPLEFORMAT, &sample_format))
-      error ("Failed to obtain a value for sample format");
-
-    if (sample_format == 3)
-      {
-        if (image_data.bits_per_sample != 32 && image_data.bits_per_sample != 64)
-          error ("Floating point images are only supported for bit depths of 32 and 64");
-      }
-    else if (sample_format != 1 && sample_format != 4)
-      error ("Unsupported sample format");
-    
-    octave_value_list retval;
-    switch (image_data.bits_per_sample)
-      {
-      case 1:
-        retval(0) = read_image<boolNDArray> (tif, &image_data);
-        break;
-      case 4:
-      case 8:
-        retval(0) = read_image<uint8NDArray> (tif, &image_data);
-        break;
-      case 16:
-        retval(0) = read_image<uint16NDArray> (tif, &image_data);
-        break;
-      case 32:
-        if (sample_format == 3)
-          retval(0) = read_image<FloatNDArray> (tif, &image_data);
-        else
-          retval(0) = read_image<uint32NDArray> (tif, &image_data);
-        break;
-      case 64:
-        if (sample_format == 3)
-          retval(0) = read_image<NDArray> (tif, &image_data);
-        else
-          retval(0) = read_image<uint64NDArray> (tif, &image_data);
-        break;
-      default:
-        error ("Unsupported bit depth");
-      }
-    
-    return retval;
-#else
-    err_disabled_feature ("read", "Tiff");
-#endif
-  }
-
-  DEFUN_DLD (__tiff_read_encoded_strip__, args, nargout,
-             "Read a strip from the image in the current IFD")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    if (nargin != 2)
-      error ("rong number of arguments");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-
-    // TODO(maged): nargout and ycbcr
-    octave_unused_parameter (nargout);
-
-    if (TIFFIsTiled (tif))
-      error ("The image is tiled not stripped");
-    
-    uint32_t strip_no;
-    if (args(1).is_scalar_type ())
-      strip_no = args(1).uint32_scalar_value ();
-    else
-      error ("Expected scalar for strip number");
-    
-    if (strip_no < 1 || strip_no > TIFFNumberOfStrips (tif))
-      error ("Strip number out of bounds");
-    
-    // Convert from Octave's 1-based indexing to zero-based indexing
-    strip_no--;
-
-    return octave_value_list (handle_read_strip_or_tile (tif, strip_no));
-#else
-    err_disabled_feature ("readEncodedStrip", "Tiff");
-#endif
-  }
-
-  DEFUN_DLD (__tiff_read_encoded_tile__, args, nargout,
-             "Read a tile from the image in the current IFD")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    if (nargin != 2)
-      error ("rong number of arguments");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-
-    // TODO(maged): nargout and ycbcr
-    octave_unused_parameter (nargout);
-
-    if (! TIFFIsTiled (tif))
-      error ("The image is stripped not tiled");
-    
-    uint32_t tile_no;
-    if (args(1).is_scalar_type ())
-      tile_no = args(1).uint32_scalar_value ();
-    else
-      error ("Expected scalar for tile number");
-    
-    if (tile_no < 1 || tile_no > TIFFNumberOfTiles (tif))
-      error ("Tile number out of bounds");
-    
-    // Convert from Octave's 1-based indexing to zero-based indexing
-    tile_no--;
-
-    return octave_value_list (handle_read_strip_or_tile (tif, tile_no));
-#else
-    err_disabled_feature ("readEncodedStrip", "Tiff");
-#endif
-  }
-
-  DEFUN_DLD (__tiff_write__, args, ,
-             "Write the image data to the current IFD")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    if (nargin < 2)
-      error ("Wrong number of arguments\n");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-
-    // TODO(maged): check on windows
-    if (TIFFGetMode (tif) == O_RDONLY)
-      error ("Can't write data to a file opened in read-only mode");
-
-    // Obtain all necessary tags
-    tiff_image_data image_data (tif);
-
-    uint16_t sample_format;
-    if (! TIFFGetFieldDefaulted(tif, TIFFTAG_SAMPLEFORMAT, &sample_format))
-      error ("Failed to obtain a value for sample format");
-
-    if (sample_format == 3)
-      {
-        if (image_data.bits_per_sample != 32
-            && image_data.bits_per_sample != 64)
-          error ("Floating point images are only supported for bit depths of 32 and 64");
-      }
-    else if (sample_format != 1 && sample_format != 4)
-      error ("Unsupported sample format");
-    
-    switch (image_data.bits_per_sample)
-      {
-      case 1:
-        // We need to check for both scalar and matrix types to handle single
-        // pixel image
-        if (args (1).is_bool_scalar () || args (1).is_bool_matrix ())
-          write_image<boolNDArray> (tif, args (1).bool_array_value (),
-                                    &image_data);
-        else
-          error ("Expected logical matrix for BiLevel image");
-        break;
-      case 8:
-        if (args (1).is_uint8_type ())
-          write_image<uint8NDArray> (tif, args (1).uint8_array_value (),
-                                     &image_data);
-        else
-          error ("Only uint8 data is allowed for uint images with bit depth of 8");
-        break;
-      case 16:
-        if (args (1).is_uint16_type ())
-          write_image<uint16NDArray> (tif, args (1).uint16_array_value (),
-                                      &image_data);
-        else
-          error ("Only uint16 data is allowed for uint images with bit depth of 16");
-        break;
-      case 32:
-        if (sample_format == 3)
-          if (args (1).is_single_type () || args (1).is_double_type ())
-            write_image<FloatNDArray> (tif, args (1).float_array_value (),
-                                       &image_data);
-          else
-            error ("Only single and double data are allowed for floating-point images");
-        else
-          if (args (1).is_uint32_type ())
-            write_image<uint32NDArray> (tif, args (1).uint32_array_value (),
-                                        &image_data);
-          else
-            error ("Only uint32 data is allowed for uint images with bit depth of 32");
-        break;
-      case 64:
-        if (sample_format == 3)
-          if (args (1).is_single_type () || args (1).is_double_type ())
-            write_image<NDArray> (tif, args (1).array_value (), &image_data);
-          else
-            error ("Only single and double data are allowed for floating-point images");
-        else  
-          if (args (1).is_uint64_type ())
-            write_image<uint64NDArray> (tif, args (1).uint64_array_value (),
-                                        &image_data);
-          else
-            error ("Only uint64 data is allowed for uint images with bit depth of 64");
-        break;
-      default:
-        error ("Unsupported bit depth");
-      }
-    
-    return octave_value_list ();
-#else
-    err_disabled_feature ("write", "Tiff");
-#endif
-  }
-
-  DEFUN_DLD (__tiff_write_encoded_strip__, args, ,
-             "Write an encoded strip to the image")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    // TODO(maged): add support for YCbCr data
-    if (nargin < 3)
-      error ("Too few arguments provided\n");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-
-    // TODO(maged): check on windows
-    if (TIFFGetMode (tif) == O_RDONLY)
-      error ("Can't write data to a file opened in read-only mode");
-
-    // Obtain all necessary tags
-    tiff_image_data image_data (tif);
-
-    if (image_data.is_tiled)
-      error ("Can't write strips to a tiled image");
-
-    uint32_t strip_no = args(1).uint32_scalar_value ();
-    if (strip_no < 1 || strip_no > TIFFNumberOfStrips (tif))
-      error ("Strip number out of range");
-    
-    handle_write_strip_or_tile (tif, strip_no, args(2), &image_data);
-    
-    return octave_value_list ();
-#else
-    err_disabled_feature ("writeEncodedStrip", "Tiff");
-#endif
-  }
-
-  DEFUN_DLD (__tiff_write_encoded_tile__, args, ,
-             "Write an encoded tile to the image")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    // TODO(maged): add support for YCbCr data
-    if (nargin < 3)
-      error ("Too few arguments provided\n");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-
-    // TODO(maged): check on windows
-    if (TIFFGetMode (tif) == O_RDONLY)
-      error ("Can't write data to a file opened in read-only mode");
-
-    // Obtain all necessary tags
-    tiff_image_data image_data (tif);
-
-    if (! image_data.is_tiled)
-      error ("Can't write tiles to a stripped image");
-
-    uint32_t tile_no = args(1).uint32_scalar_value ();
-    if (tile_no < 1 || tile_no > TIFFNumberOfTiles (tif))
-      error ("Tile number out of range");
-
-    handle_write_strip_or_tile (tif, tile_no, args(2), &image_data);
-    
-    return octave_value_list ();
-#else
-    err_disabled_feature ("writeEncodedTile", "Tiff");
-#endif
-  }
-  
-  DEFUN_DLD (__tiff_is_tiled__, args, ,
-             "Get whether the image is tiled")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    if (nargin == 0)
-      error ("No handle provided\n");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-    bool is_tiled = static_cast<bool> (TIFFIsTiled (tif));
-    return octave_value_list (octave_value (is_tiled));
-#else
-    err_disabled_feature ("isTiled", "Tiff");
-#endif
-  }
-
-  DEFUN_DLD (__tiff_number_of_strips__, args, ,
-             "Get the number of strips in the image")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    if (nargin == 0)
-      error ("No handle provided\n");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-    
-    if (TIFFIsTiled (tif))
-      error ("The image is tiled not stripped");
-    
-    double strip_count = static_cast<double> (TIFFNumberOfStrips (tif));
-    return octave_value_list (octave_value (strip_count));
-#else
-    err_disabled_feature ("numberOfStrips", "Tiff");
-#endif
-  }
-
-  DEFUN_DLD (__tiff_number_of_tiles__, args, ,
-             "Get the number of tiles in the image")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    if (nargin == 0)
-      error ("No handle provided\n");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-    
-    if (! TIFFIsTiled (tif))
-      error ("The image is stripped not tiled");
-    
-    double tile_count = static_cast<double> (TIFFNumberOfTiles (tif));
-    return octave_value_list (octave_value (tile_count));
-#else
-    err_disabled_feature ("numberOfTiles", "Tiff");
-#endif
-  }
-
-  DEFUN_DLD (__tiff_compute_strip__, args, ,
-             "Get the strip index containing the given row")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    if (nargin < 2 || nargin > 3)
-      error ("Wrong number of arguments\n");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-    
-    if (TIFFIsTiled (tif))
-      error ("The image is tiled not stripped");
-    
-    tiff_image_data image_data (tif);
-
-    uint32_t row = args(1).uint32_scalar_value ();
-    if (row > image_data.height)
-      row = image_data.height;
-
-    // Convert from 1-based to zero-based indexing but avoid underflow
-    if (row > 0)
-      row--;
-    
-    uint16_t plane;
-    if (nargin > 2)
-      {
-        if (image_data.planar_configuration == PLANARCONFIG_CONTIG)
-          error ("Can't use plane argument for images with chunky PlanarConfiguration");
-        plane = args(2).uint16_scalar_value ();
-        if (plane > image_data.samples_per_pixel)
-          plane = image_data.samples_per_pixel;
-        if (plane > 0)
-          plane--;
-      }
-    else
-      {
-        if (image_data.planar_configuration == PLANARCONFIG_SEPARATE)
-          error ("The plane argument is required for images with separate PlanarConfiguration");
-        plane = 0;
-      }
-
-    double strip_number = TIFFComputeStrip (tif, row, plane) + 1;
-    if (strip_number > TIFFNumberOfStrips (tif))
-      strip_number = TIFFNumberOfStrips (tif);
-    return octave_value_list (octave_value (strip_number));
-#else
-    err_disabled_feature ("computeStrip", "Tiff");
-#endif
-  }
-
-  DEFUN_DLD (__tiff_compute_tile__, args, ,
-             "Get the tile index containing the given row and column")
-  {
-#if defined (HAVE_TIFF)
-    int nargin = args.length ();
-
-    if (nargin < 2 || nargin > 3)
-      error ("Wrong number of arguments\n");
-    
-    TIFF *tif = (TIFF *)(args(0).uint64_value ());
-    
-    if (! TIFFIsTiled (tif))
-      error ("The image is stripped not tiled");
-    
-    uint32NDArray coords = args(1).uint32_array_value ();
-    if (coords.dim2() < 2)
-      error ("Coordinates must be in the shape [row, col]");
-    uint32_t row = coords(0, 0);
-    uint32_t col = coords(0, 1);
-
-    tiff_image_data image_data (tif);
-
-    if (col > image_data.width)
-      col = image_data.width;
-    if (row > image_data.height)
-      row = image_data.height;
-
-    // Convert from 1-based to zero-based indexing but avoid underflow
-    if (row > 0)
-      row--;
-    if (col > 0)
-      col--;
-    
-    uint16_t plane;
-    if (nargin > 2)
-      {
-        if (image_data.planar_configuration == PLANARCONFIG_CONTIG)
-          error ("Can't use plane argument for images with chunky PlanarConfiguration");
-        plane = args(2).uint16_scalar_value ();
-        if (plane > image_data.samples_per_pixel)
-          plane = image_data.samples_per_pixel;
-        if (plane > 0)
-          plane--;
-      }
-    else
-      {
-        if (image_data.planar_configuration == PLANARCONFIG_SEPARATE)
-          error ("The plane argument is required for images with separate PlanarConfiguration");
-        plane = 0;
-      }
-
-    double tile_number = TIFFComputeTile (tif, col, row, 0, plane) + 1;
-    if (tile_number > TIFFNumberOfTiles (tif))
-      tile_number = TIFFNumberOfTiles (tif);
-    return octave_value_list (octave_value (tile_number));
-#else
-    err_disabled_feature ("computeTile", "Tiff");
-#endif
-  }
-
-  DEFUN_DLD (__tiff_version__, , ,
-             "Get the version stamp of LibTIFF")
-  {
-#if defined (HAVE_TIFF)
-    std::string version = std::string (TIFFGetVersion ());
-    return octave_value_list (octave_value (version));
-#else
-    err_disabled_feature ("getVersion", "Tiff");
-#endif
-  }
-
-}