#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include <string>
#include <iostream>
#include <fstream>

#include "defun.h"
#include "ov.h"
#include "ovl.h"
#include "error.h"

#include "errwarn.h"

#include "fcntl-wrappers.h"
#include "file-stat.h"
#include "file-ops.h"
#include "str-vec.h"
#include "oct-time.h"

#if defined (HAVE_TIFF)
#  include <tiffio.h>
bool handlers_set = true;
#endif

namespace octave
{
#if defined (HAVE_TIFF)
  class OCTINTERP_API octave_tiff_handle : public octave_base_value
  {
  public:
    octave_tiff_handle (TIFF *tif)
      : tiff_file (tif), closed (false)
    { }

    TIFF *get_file(void) { return tiff_file; }

    bool is_defined (void) const { return true; }
    bool is_constant (void) const { return true; }

    void close(void)
    {
      if (! closed)
        {
          closed = true;
          TIFFClose (tiff_file);
        }
    }

    bool is_closed (void) { return closed; }

    ~octave_tiff_handle (void)
    {
      if (! closed)
        close ();
    }

    static octave_tiff_handle *
    get_tiff_handle (const octave_value& ov)
    {
      octave_base_value *rep = ov.internal_rep ();
      octave_tiff_handle *handle = dynamic_cast<octave_tiff_handle *> (rep);
      if (! handle)
        error ("get_tiff_handle: dynamic_cast to octave_tiff_handle failed");
      
      return handle;
    }
  private:
    TIFF *tiff_file;
    bool closed;
  };

  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");
      
      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
        // This doesn't completely solve the issue as LibTIFF will refuse
        // to write data to images of the tag is not set.
        // see: https://www.asmail.be/msg0054918184.html
        // This is now solved by this commit:
        // https://gitlab.com/libtiff/libtiff/-/commit/ac6ddaf678fff597610cb217705b55508240f065
        // So this shouldn't be needed for future LibTIFF releases > 4.4.0
        planar_configuration = 1;
      
      is_tiled = TIFFIsTiled(tif);
    }
  };

  void
  set_internal_handlers (void)
  {
    if (handlers_set)
      {
        TIFFSetErrorHandler (verror_with_id);
        TIFFSetWarningHandler (vwarning_with_id);
      }
    else
      {
        TIFFSetErrorHandler (NULL);
        TIFFSetWarningHandler (NULL);
      }
  }

  bool
  is_colormap (octave_value ov)
  {
    return ov.isnumeric () && ov.isreal () && ov.isfloat ()
           && ov.ndims () == 2 && ov.columns () == 3;
  }

  // A map of tag names supported by matlab, there are some differences
  // than LibTIFF's names (e.g. Photometric vs PhotometricInterpretation)
  static const std::map<std::string, ttag_t> tag_name_map = {
    {"SubFileType", 254},
    {"ImageWidth", 256},
    {"ImageLength", 257},
    {"BitsPerSample", 258},
    {"Compression", 259},
    {"Photometric", 262},
    {"Thresholding", 263},
    {"FillOrder", 266},
    {"DocumentName", 269},
    {"ImageDescription", 270},
    {"Make", 271},
    {"Model", 272},
    {"StripOffsets", 273},
    {"Orientation", 274},
    {"SamplesPerPixel", 277},
    {"RowsPerStrip", 278},
    {"StripByteCounts", 279},
    {"MinSampleValue", 280},
    {"MaxSampleValue", 281},
    {"XResolution", 282},
    {"YResolution", 283},
    {"PlanarConfiguration", 284},
    {"PageName", 285},
    {"XPosition", 286},
    {"YPosition", 287},
    {"GrayResponseUnit", 290},
    {"GrayResponseCurve", 291},
    {"Group3Options", 292},
    {"Group4Options", 293},
    {"ResolutionUnit", 296},
    {"PageNumber", 297},
    {"TransferFunction", 301},
    {"Software", 305},
    {"DateTime", 306},
    {"Artist", 315},
    {"HostComputer", 316},
    {"WhitePoint", 318},
    {"PrimaryChromaticities", 319},
    {"ColorMap", 320},
    {"HalfToneHints", 321},
    {"TileWidth", 322},
    {"TileLength", 323},
    {"TileOffsets", 324},
    {"TileByteCounts", 325},
    {"SubIFD", 330},
    {"InkSet", 332},
    {"InkNames", 333},
    {"NumberOfInks", 334},
    {"DotRange", 336},
    {"TargetPrinter", 337},
    {"ExtraSamples", 338},
    {"SampleFormat", 339},
    {"SMinSampleValue", 340},
    {"SMaxSampleValue", 341},
    {"YCbCrCoefficients", 529},
    {"YCbCrSubSampling", 530},
    {"YCbCrPositioning", 531},
    {"ReferenceBlackWhite", 532},
    {"XMP", 700},
    {"ImageDepth", 32997},
    {"Copyright", 33432},
    {"RichTIFFIPTC", 33723},
    {"Photoshop", 34377},
    {"ICCProfile", 34675},
    {"SToNits", 37439},
    {"JPEGQuality", 65537},
    {"JPEGColorMode", 65538},
    {"ZipQuality", 65557},
    {"SGILogDataFmt", 6556}
  };

  const TIFFField *get_fip_with_name (TIFF *tif, std::string& tag_name)
  {
    const TIFFField *fip = NULL;
    if (tag_name_map.find (tag_name) != tag_name_map.cend ())
      fip = TIFFFieldWithTag (tif, tag_name_map.at (tag_name));
    
    // If the tag name is not found, fallback to the names defined
    // in LibTIFF.
    if (fip == NULL)
      fip = TIFFFieldWithName (tif, tag_name.c_str ());
    
    return fip;
  }

  void
  check_readonly (TIFF *tif)
  {
    if (TIFFGetMode (tif) == octave_o_rdonly_wrapper ())
      error ("Can't write data to a file opened in read-only mode");
  }

  void
  check_closed (octave_tiff_handle *tiff_handle)
  {
    if (tiff_handle->is_closed ())
      error ("The image file was closed");
  }

  void
  validate_tiff_get_field (bool status)
  {
    // Check if the return status of TIFFGetField is not 1
    if (status != 1)
      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;
  }

  void
  get_tile_dimensions_validated (TIFF *tif, uint32_t& tile_width,
                                 uint32_t& 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");
  }

  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");
    
    // The default for RowsPerStrip is INT_MAX so it should be capped to the
    // image height
    if (rows_in_strip > image_data->height)
      rows_in_strip = image_data->height;
    
    uint32_t strip_count = TIFFNumberOfStrips (tif);
    // Get the actual number of rows in the strip
    rows_in_strip = get_rows_in_strip (strip_no, strip_count,
                                       rows_in_strip, image_data);
    
    // Set the strip dimensions depending on the planar configuration
    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 = col % 8;
                if (TIFFIsMSB2LSB (tif))
                  shift = 7 - shift;
              strip_fvec[col] = (strip_buf[col / 8] >> shift) & 0x1;
            }
            strip_fvec += image_data->width;
            strip_buf += padded_width;
          }
      }
    else
      error ("Unsupported bit depth");
    
    // Reorder the dimensions to the order expected by Octave
    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

    uint32_t tile_width, tile_height;
    get_tile_dimensions_validated (tif, tile_width, tile_height);

    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 = col % 8;
                if (TIFFIsMSB2LSB (tif))
                  shift = 7 - shift;
                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 for boundary tiles
    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
  read_unsigned_strip_or_tile (TIFF *tif, uint32_t strip_tile_no,
                               tiff_image_data *image_data)
  {
    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:
        return read_strip_or_tile<uint32NDArray> (tif, strip_tile_no,
                                                  image_data);
      case 64:
        return read_strip_or_tile<uint64NDArray> (tif, strip_tile_no,
                                                  image_data);
      default:
        error ("Unsupported bit depth");
      }
  }
  
  octave_value
  read_signed_strip_or_tile (TIFF *tif, uint32_t strip_tile_no,
                             tiff_image_data *image_data)
  {
    switch (image_data->bits_per_sample)
      {
      case 8:
        return read_strip_or_tile<int8NDArray> (tif, strip_tile_no,
                                                image_data);
        break;
      case 16:
        return read_strip_or_tile<int16NDArray> (tif, strip_tile_no,
                                                 image_data);
        break;
      case 32:
        return read_strip_or_tile<int32NDArray> (tif, strip_tile_no,
                                                 image_data);
      case 64:
        return read_strip_or_tile<int64NDArray> (tif, strip_tile_no,
                                                 image_data);
      default:
        error ("Unsupported bit depth for signed images");
      }
  }
  
  octave_value
  read_float_strip_or_tile (TIFF *tif, uint32_t strip_tile_no,
                            tiff_image_data *image_data)
  {
    switch (image_data->bits_per_sample)
      {
      case 32:
        return read_strip_or_tile<FloatNDArray> (tif, strip_tile_no,
                                                 image_data);
      case 64:
        return read_strip_or_tile<NDArray> (tif, strip_tile_no,
                                            image_data);
      default:
        error ("Unsupported bit depth for floating-point images");
      }
  }
  
  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");

    switch (sample_format)
      {
      case 1:
      case 4:
        return read_unsigned_strip_or_tile (tif, strip_tile_no, &image_data);
        break;
      case 2:
        return read_signed_strip_or_tile (tif, strip_tile_no, &image_data);
        break;
      case 3:
        return read_float_strip_or_tile (tif, strip_tile_no, &image_data);
        break;
      default:
        error ("Unsupported sample format");
      }
  }

  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--)
              {
                uint8_t bit_number = pixel % 8;
                if (TIFFIsMSB2LSB (tif))
                  bit_number = 7 - bit_number;
                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--)
              {
                uint8_t bit_number = pixel % 8;
                if (TIFFIsMSB2LSB (tif))
                  bit_number = 7 - bit_number;
                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);
  }

  octave_value
  read_unsigned_image (TIFF *tif, tiff_image_data *image_data)
  {
    switch (image_data->bits_per_sample)
      {
      case 1:
        return read_image<boolNDArray> (tif, image_data);
        break;
      case 4:
      case 8:
        return read_image<uint8NDArray> (tif, image_data);
        break;
      case 16:
        return read_image<uint16NDArray> (tif, image_data);
        break;
      case 32:
        return read_image<uint32NDArray> (tif, image_data);
      case 64:
        return read_image<uint64NDArray> (tif, image_data);
        break;
      default:
        error ("Unsupported bit depth");
      }
  }

  octave_value
  read_signed_image (TIFF *tif, tiff_image_data *image_data)
  {
    switch (image_data->bits_per_sample)
      {
      case 8:
        return read_image<int8NDArray> (tif, image_data);
        break;
      case 16:
        return read_image<int16NDArray> (tif, image_data);
        break;
      case 32:
        return read_image<int32NDArray> (tif, image_data);
      case 64:
        return read_image<int64NDArray> (tif, image_data);
        break;
      default:
        error ("Unsupported bit depth for signed images");
      }
  }

  octave_value
  read_float_image (TIFF *tif, tiff_image_data *image_data)
  {
    switch (image_data->bits_per_sample)
      {
      case 32:
        return read_image<FloatNDArray> (tif, image_data);
      case 64:
        return read_image<NDArray> (tif, image_data);
        break;
      default:
        error ("Unsupported bit depth for floating-point images");
      }
  }

  template <typename T>
  octave_value
  make_array (void *data, dim_vector& arr_dims)
  {
    typedef typename T::element_type P;

    T arr (arr_dims);
    for (uint32_t i = 0; i < arr_dims.numel (); i++)
      {
        arr(i) = (reinterpret_cast<P *> (data))[i];
      }
    return octave_value (arr);
  }

  octave_value
  interpret_tag_data (void *data, uint32_t count, TIFFDataType tag_datatype,
                      bool convert_to_double = false)
  {
    octave_value retval;dim_vector arr_dims (1, count);
    switch (tag_datatype)
      {
      case TIFF_BYTE:
      case TIFF_UNDEFINED:
        {
          retval = make_array<uint8NDArray> (data, arr_dims);
          break;
        }
      case TIFF_ASCII:
        {
          retval = octave_value (*(reinterpret_cast<char **> (data)));
          break;
        }
      case TIFF_SHORT:
        {
          retval = make_array<uint16NDArray> (data, arr_dims);
          break;
        }
      case TIFF_LONG:
      case TIFF_IFD:
        {
          retval = make_array<uint32NDArray> (data, arr_dims);
          break;
        }
      case TIFF_LONG8:
      case TIFF_IFD8:
        {
          retval = make_array<uint64NDArray> (data, arr_dims);
          break;
        }
      case TIFF_RATIONAL:
      case TIFF_SRATIONAL:
      case TIFF_FLOAT:
        {
          retval = make_array<FloatNDArray> (data, arr_dims);
          break;
        }
      case TIFF_SBYTE:
        {
          retval = make_array<int8NDArray> (data, arr_dims);
          break;
        }
      case TIFF_SSHORT:
        {
          retval = make_array<int16NDArray> (data, arr_dims);
          break;
        }
      case TIFF_SLONG:
        {
          retval = make_array<int32NDArray> (data, arr_dims);
          break;
        }
      case TIFF_SLONG8:
        {
          retval = make_array<int64NDArray> (data, arr_dims);
          break;
        }
      case TIFF_DOUBLE:
        {
          retval = make_array<NDArray> (data, arr_dims);
          break;
        }
      default:
        error ("Unsupported tag data type");
      }

    if (tag_datatype != TIFF_ASCII && convert_to_double)
      retval = retval.as_double ();
    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));
    
    std::unique_ptr<uint8_t []> data
      = std::make_unique<uint8_t []> (type_size);
    uint8_t *data_ptr = data.get ();
    if (tag_id == TIFFTAG_PLANARCONFIG)
      {
        // Workaround for a bug in LibTIFF where it incorrectly returns
        // zero as the default value for PlanaConfiguration
        // See: https://www.asmail.be/msg0054918184.html
        if (! TIFFGetField(tif, tag_id, data_ptr))
          *reinterpret_cast<uint16_t *> (data_ptr) = 1;
      }
    else
      validate_tiff_get_field (TIFFGetFieldDefaulted (tif, tag_id, data_ptr));
    
    // All scalar tags are returned as double
    octave_value tag_data_ov = interpret_tag_data (data_ptr, 1,
                                                   TIFFFieldDataType (fip),
                                                   true);

    return tag_data_ov;
  }

  octave_value
  get_array_field_data (TIFF *tif, const TIFFField *fip, uint32_t array_size,
                        bool convert_to_double = false)
  {
    void *data;
    validate_tiff_get_field (TIFFGetField (tif, TIFFFieldTag (fip), &data));
      
    return interpret_tag_data (data, array_size, TIFFFieldDataType (fip),
                               convert_to_double);
  }

  octave_value
  get_field_data (TIFF *tif, const TIFFField *fip)
  {
    octave_value tag_data_ov;
    uint32_t tag_id = TIFFFieldTag (fip);

    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, true);
        break;
      case TIFFTAG_REFERENCEBLACKWHITE:
        tag_data_ov = get_array_field_data (tif, fip, 6, true);
        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, true);
          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 bit depth");

          uint32_t count = 1 << bits_per_sample;
          uint16_t *ch1 = NULL, *ch2 = NULL, *ch3 = NULL;
          // TransferFunction can have 1 or 3 channels depending on the
          // SamplesPerPixel but the library still expects three pointers
          // to be passed and sets only the needed ones
          validate_tiff_get_field (TIFFGetField (tif, tag_id,
                                                 &ch1, &ch2, &ch3));
          if (ch2 == NULL)
            tag_data_ov = interpret_tag_data (ch1, count,
                                              TIFFFieldDataType (fip), true);
          else
            {
              // Concatenate the channels into an n by 3 array
              OCTAVE_LOCAL_BUFFER (NDArray, array_list, 3);
              dim_vector col_dims(count, 1);
              array_list[0] = interpret_tag_data (ch1,
                                                  count,
                                                  TIFFFieldDataType (fip))
                                                  .array_value ()
                                                  .reshape (col_dims);
              array_list[1] = interpret_tag_data (ch2,
                                                  count,
                                                  TIFFFieldDataType (fip))
                                                  .array_value ()
                                                  .reshape (col_dims);
              array_list[2] = interpret_tag_data (ch3,
                                                  count,
                                                  TIFFFieldDataType (fip))
                                                  .array_value ()
                                                  .reshape (col_dims);
              
              tag_data_ov
                = octave_value (NDArray::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;
        }
      case TIFFTAG_XMLPACKET:
      case TIFFTAG_RICHTIFFIPTC:
      case TIFFTAG_PHOTOSHOP:
      case TIFFTAG_ICCPROFILE:
        {
          uint32_t count;
          uint8_t *data;
          validate_tiff_get_field (TIFFGetField (tif, tag_id,
                                                 &count, &data));
          tag_data_ov = interpret_tag_data (data, count, TIFF_BYTE);
          // Matlab returns XMP tag as char array not byte array
          if (tag_id == TIFFTAG_XMLPACKET)
            tag_data_ov = octave_value (tag_data_ov.char_array_value ());
          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_scalar_field_data (TIFF *tif, const TIFFField *fip, octave_value tag_ov)
  {
    uint32_t tag_id = TIFFFieldTag (fip);
    
    // 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");
    
    TIFFDataType tag_datatype = TIFFFieldDataType (fip);

    // According to matlab, the value must be a scalar double
    // except for strings
    if (tag_datatype == TIFF_ASCII)
      {
        if (! tag_ov.is_string ())
          error ("Expected string for ascii tag");
      }
    else if (! tag_ov.is_scalar_type () || ! tag_ov.is_double_type ())
      error ("Tag value must be a scalar double");
    
    switch (tag_datatype)
      {
      case TIFF_BYTE:
      case TIFF_UNDEFINED:
        TIFFSetField (tif, tag_id, tag_ov.uint8_scalar_value ());
        break;
      case TIFF_ASCII:
        TIFFSetField (tif, tag_id, tag_ov.string_value ().c_str ());
        break;
      case TIFF_SHORT:
        TIFFSetField (tif, tag_id, tag_ov.uint16_scalar_value ());
        break;
      case TIFF_LONG:
        TIFFSetField (tif, tag_id, tag_ov.uint32_scalar_value ());
        break;
      case TIFF_LONG8:
        TIFFSetField (tif, tag_id, tag_ov.uint64_scalar_value ());
        break;
      case TIFF_RATIONAL:
        TIFFSetField (tif, tag_id, tag_ov.float_scalar_value ());
        break;
      case TIFF_SBYTE:
        TIFFSetField (tif, tag_id, tag_ov.int8_scalar_value ());
        break;
      case TIFF_SSHORT:
        TIFFSetField (tif, tag_id, tag_ov.int16_scalar_value ());
        break;
      case TIFF_SLONG:
        TIFFSetField (tif, tag_id, tag_ov.int32_scalar_value ());
        break;
      case TIFF_SLONG8:
        TIFFSetField (tif, tag_id, tag_ov.int64_scalar_value ());
        break;
      case TIFF_FLOAT:
        TIFFSetField (tif, tag_id, tag_ov.float_scalar_value ());
        break;
      case TIFF_DOUBLE:
        TIFFSetField (tif, tag_id, tag_ov.double_value ());
        break;
      case TIFF_SRATIONAL:
        TIFFSetField (tif, tag_id, tag_ov.float_scalar_value ());
        break;
      case TIFF_IFD:
      case TIFF_IFD8:
        error ("Unimplemented IFFD data type");
        break;
      default:
        error ("Unsupported tag data type");
      }
  }

  template <typename T>
  void
  set_array_field_helper (TIFF *tif, uint32_t tag_id, T data_arr)
  {
    const void *data_ptr = data_arr.data ();

    if (! TIFFSetField (tif, tag_id, data_ptr))
      error ("Failed to set tag");
  }

  void
  set_array_field_data (TIFF *tif, const TIFFField *fip,
                        octave_value tag_ov, uint32_t count)
  {
    uint32_t tag_id = TIFFFieldTag (fip);
    TIFFDataType tag_datatype = TIFFFieldDataType (fip);

    // Array type must be double in matlab
    if (! tag_ov.is_double_type ())
      error ("Tag data must be double");

    // Matlab checks number of elements not dimensions
    if (tag_ov.array_value ().numel () != count)
      error ("Expected an array with %u elements", count);

    switch (tag_datatype)
      {
      case TIFF_BYTE:
      case TIFF_UNDEFINED:
        set_array_field_helper<uint8NDArray> (tif, tag_id,
                                              tag_ov.uint8_array_value ());
        break;
      case TIFF_SHORT:
        set_array_field_helper<uint16NDArray> (tif, tag_id,
                                               tag_ov.uint16_array_value ());
        break;
      case TIFF_LONG:
        set_array_field_helper<uint32NDArray> (tif, tag_id,
                                               tag_ov.uint32_array_value ());
        break;
      case TIFF_LONG8:
        set_array_field_helper<uint64NDArray> (tif, tag_id,
                                               tag_ov.uint64_array_value ());
        break;
      case TIFF_RATIONAL:
        set_array_field_helper<FloatNDArray> (tif, tag_id,
                                               tag_ov.float_array_value ());
        break;
      case TIFF_SBYTE:
        set_array_field_helper<int8NDArray> (tif, tag_id,
                                             tag_ov.int8_array_value ());
        break;
      case TIFF_SSHORT:
        set_array_field_helper<int16NDArray> (tif, tag_id,
                                              tag_ov.int16_array_value ());
        break;
      case TIFF_SLONG:
        set_array_field_helper<int32NDArray> (tif, tag_id,
                                              tag_ov.int32_array_value ());
        break;
      case TIFF_SLONG8:
        set_array_field_helper<int64NDArray> (tif, tag_id,
                                              tag_ov.int64_array_value ());
        break;
      case TIFF_FLOAT:
        set_array_field_helper<FloatNDArray> (tif, tag_id,
                                              tag_ov.float_array_value ());
        break;
      case TIFF_DOUBLE:
        set_array_field_helper<NDArray> (tif, tag_id,
                                         tag_ov.array_value ());
        break;
      case TIFF_SRATIONAL:
        set_array_field_helper<FloatNDArray> (tif, tag_id,
                                              tag_ov.float_array_value ());
        break;
      case TIFF_IFD:
      case TIFF_IFD8:
        error ("Unimplemented IFFD data type");
        break;
      default:
        error ("Unsupported tag data type");
      }
  }

  void
  set_field_data (TIFF *tif, const TIFFField *fip, octave_value tag_ov)
  {
    uint32_t tag_id = TIFFFieldTag (fip);
    
    switch (tag_id)
      {
      case TIFFTAG_YCBCRCOEFFICIENTS:
        set_array_field_data (tif, fip, tag_ov, 3);
        break;
      case TIFFTAG_REFERENCEBLACKWHITE:
        set_array_field_data (tif, fip, tag_ov, 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");
          
          set_array_field_data (tif, fip, tag_ov, 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");
          
          if (! tag_ov.is_double_type ())
            error ("Tag data must be double");
          
          // 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;

          if (tag_ov.is_scalar_type ()
              || tag_ov.array_value ().dim1 () != count
              || tag_ov.array_value ().dim2 () != 3)
            error ("Expected matrix with %u rows and 3 columns", count);

          NDArray array_data = tag_ov.array_value ();
          array_data *= UINT16_MAX;
          uint16NDArray u16_array = uint16NDArray (array_data);
          const uint16_t *data_ptr
            = reinterpret_cast<const uint16_t *> (u16_array.data ());
          if (! TIFFSetField (tif, tag_id, data_ptr, data_ptr + count,
                              data_ptr + 2 * count))
            error ("Failed to set tag");
          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 bit depth");

          uint32_t count = 1 << bits_per_sample;
          
          // An intermediate variable is required for storing the return of
          // uint16_array_value so that the object remains in scope and the
          // pointer returned from data () is not a dangling pointer
          uint16NDArray data_arr = tag_ov.uint16_array_value ();

          if (data_arr.numel () != count)
            error ("Tag data should be and array of %u elements", count);

          const uint16_t *data_ptr
            = reinterpret_cast<const uint16_t *> (data_arr.data ());
          if (samples_per_pixel == 1)
            {
              if (! TIFFSetField (tif, tag_id, data_ptr))
                error ("Failed to set field");
            }
          else
            {
              if (! TIFFSetField (tif, tag_id, data_ptr, data_ptr + count,
                                  data_ptr + 2 * count))
                error ("Failed to set field");
            }
          break;
        }
      case TIFFTAG_PAGENUMBER:
      case TIFFTAG_HALFTONEHINTS:
      case TIFFTAG_DOTRANGE:
      case TIFFTAG_YCBCRSUBSAMPLING:
        {
          uint16NDArray array_data = tag_ov.uint16_array_value ();
          if (array_data.numel () != 2)
            error ("Tag data should be an array with 2 elements");
          
          if (! TIFFSetField (tif, tag_id, array_data (0), array_data (1)))
            error ("Failed to set field");
          break;
        }
      case TIFFTAG_SUBIFD:
        {
          // Setting the SubIFD invloves setting their count only
          // and the library will set the offsets correctly in
          // the next calls to writeDirectory, but an array of offsets
          // is still required so an array of zeroes is passed
          uint16_t subifd_count = tag_ov.uint16_scalar_value ();
          std::unique_ptr<uint64_t []> subifd_offsets
            = std::make_unique<uint64_t []> (subifd_count);
          uint64_t *offfsets_ptr = subifd_offsets.get ();
          if (! TIFFSetField (tif, tag_id, subifd_count, offfsets_ptr))
            error ("Failed to set tag");
          break;
        }
      case TIFFTAG_EXTRASAMPLES:
        {
          uint16_t samples_per_pixel;
          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_SAMPLESPERPIXEL,
                                       &samples_per_pixel))
            error ("Failed to obtain the number of samples per pixel");
          
          uint16NDArray data_array = tag_ov.uint16_array_value ();
          // FIXME: this only works for RGB images. Need to handle grayscale,
          // palette, cmyk and ycbcr
          if (data_array.numel () > samples_per_pixel - 3)
            error ("Failed to set field, too many values");
          
          if (! TIFFSetField (tif, tag_id,
                              static_cast<uint16_t> (data_array.numel ()),
                              data_array.data ()))
            error ("Failed to set field");
          break;
        }
      case TIFFTAG_XMLPACKET:
      case TIFFTAG_RICHTIFFIPTC:
      case TIFFTAG_PHOTOSHOP:
      case TIFFTAG_ICCPROFILE:
        {
          uint8NDArray data_array = tag_ov.uint8_array_value ();
          uint32_t count = data_array.numel ();
          if (! TIFFSetField (tif, tag_id, count, data_array.data ()))
            error ("Failed to set field");
          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:
        {
          set_scalar_field_data (tif, fip, tag_ov);
          break;
        }
      default:
        set_scalar_field_data (tif, fip, tag_ov);
    }
    
  }
  
  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);

    uint8_t *data_u8
      = 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)
      {
        // 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_u8, 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 ();
        // 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 = col % 8;
                if (TIFFIsMSB2LSB (tif))
                  shift = 7 - shift;
              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;
    get_tile_dimensions_validated (tif, tile_width, tile_height);
    
    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--;
    uint8_t *data_u8 = 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 (TIFFWriteEncodedTile (tif, tile_no, data_u8,
                                  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 ();
        // 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 = col % 8;
                if (TIFFIsMSB2LSB (tif))
                  shift = 7 - shift;
              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
  write_unsigned_strip_or_tile (TIFF *tif, uint32_t strip_tile_no,
                                octave_value data_ov,
                                tiff_image_data *image_data)
  {
    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 (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 (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");
      }
  }
  
  void
  write_signed_strip_or_tile (TIFF *tif, uint32_t strip_tile_no,
                              octave_value data_ov,
                              tiff_image_data *image_data)
  {
    switch (image_data->bits_per_sample)
      {
      case 8:
        if (data_ov.is_int8_type ())
          write_strip_or_tile<int8NDArray> (tif, strip_tile_no,
                                            data_ov.int8_array_value (),
                                            image_data);
        else
          error ("Only int8 data is allowed for int images with bit depth of 8");
        break;
      case 16:
        if (data_ov.is_int16_type ())
          write_strip_or_tile<int16NDArray> (tif, strip_tile_no,
                                             data_ov.int16_array_value (),
                                             image_data);
        else
          error ("Only int16 data is allowed for int images with bit depth of 16");
        break;
      case 32:
        if (data_ov.is_int32_type ())
          write_strip_or_tile<int32NDArray> (tif, strip_tile_no,
                                             data_ov.int32_array_value (),
                                             image_data);
        else
          error ("Only int32 data is allowed for int images with bit depth of 32");
        break;
      case 64:  
        if (data_ov.is_int64_type ())
          write_strip_or_tile<int64NDArray> (tif, strip_tile_no,
                                             data_ov.int64_array_value (),
                                             image_data);
        else
          error ("Only int64 data is allowed for int images with bit depth of 64");
        break;
      default:
        error ("Unsupported bit depth for signed images");
      }
  }
  
  void
  write_float_strip_or_tile (TIFF *tif, uint32_t strip_tile_no,
                             octave_value data_ov,
                             tiff_image_data *image_data)
  {
    switch (image_data->bits_per_sample)
      {
      case 32:
        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");
        break;
      case 64:  
        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");
        break;
      default:
        error ("Unsupported bit depth for floating-point images");
      }
  }
  
  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");

    switch (sample_format)
      {
      case 1:
      case 4:
        write_unsigned_strip_or_tile (tif, strip_tile_no, data_ov,
                                      image_data);
        break;
      case 2:
        write_signed_strip_or_tile (tif, strip_tile_no, data_ov, image_data);
        break;
      case 3:
        write_float_strip_or_tile (tif, strip_tile_no, data_ov, image_data);
        break;
      default:
        error ("Unsupported sample format");
      }
  }

  template <typename T>
  void
  write_stripped_image (TIFF *tif, T pixel_data, tiff_image_data *image_data)
  {
    // 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) == -1)
              error ("Failed to write 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 = col % 8;
                  if (TIFFIsMSB2LSB (tif))
                    shift = 7 - shift;
                  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)
  {
    // ASSUMES pixel data dimensions are already validated

    uint32_t tile_width, tile_height;
    get_tile_dimensions_validated (tif, tile_width, tile_height);
    
    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) == -1)
              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 = col % 8;
                  if (TIFFIsMSB2LSB (tif))
                    shift = 7 - shift;
                  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)
  {
    // This dimensions checking is intentially done in this non-homogeneous
    // way for matlab compatibility.
    // In Matlab R2022a, If the width or height of the input matrix is less
    // than needed, the data is silently padded with zeroes to fit.
    // If the width is larger than needed, a warning is issued and the excess
    // data is truncated. If the height is larger than needed, no warning is
    // issued and the image data is wrong. If the number of channels is less
    // or more than needed an error is produced.
    // We chose to deviate from matlab in the larger height case to avoid
    // errors resulting from silently corrupting image data, a warning is
    // produced instead.
    if ((image_data->samples_per_pixel > 1 && pixel_data.ndims () < 3)
        || (pixel_data.ndims () > 2
            && image_data->samples_per_pixel != pixel_data.dim3 ()))
      error ("Incorrect number of channels, expected %u",
             image_data->samples_per_pixel);
    
    if (pixel_data.dim1 () > image_data->height)
      warning ("Input has more rows than the image length, data will be truncated");
    if (pixel_data.dim2 () > image_data->width)
      warning ("Input has more columns than the image width, data will be truncated");
    
    pixel_data.resize (dim_vector (image_data->height, image_data->width,
                                   image_data->samples_per_pixel));
    
    if (image_data->is_tiled)
      write_tiled_image<T> (tif, pixel_data, image_data);
    else
      write_stripped_image<T> (tif, pixel_data, image_data);

  }

  void
  write_unsigned_image (TIFF *tif, octave_value image_ov,
                        tiff_image_data *image_data)
  {
    switch (image_data->bits_per_sample)
      {
      case 1:
        // We need to check for both scalar and matrix types to handle single
        // pixel image
        if (image_ov.is_bool_scalar () || image_ov.is_bool_matrix ())
          write_image<boolNDArray> (tif, image_ov.bool_array_value (),
                                    image_data);
        else
          error ("Expected logical matrix for BiLevel image");
        break;
      case 8:
        if (image_ov.is_uint8_type ())
          write_image<uint8NDArray> (tif, image_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 (image_ov.is_uint16_type ())
          write_image<uint16NDArray> (tif, image_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 (image_ov.is_uint32_type ())
          write_image<uint32NDArray> (tif, image_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 (image_ov.is_uint64_type ())
          write_image<uint64NDArray> (tif, image_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");
      }
  }

  void
  write_signed_image (TIFF *tif, octave_value image_ov,
                      tiff_image_data *image_data)
  {
    switch (image_data->bits_per_sample)
      {
      case 8:
        if (image_ov.is_int8_type ())
          write_image<int8NDArray> (tif, image_ov.int8_array_value (),
                                    image_data);
        else
          error ("Only int8 data is allowed for int images with bit depth of 8");
        break;
      case 16:
        if (image_ov.is_int16_type ())
          write_image<int16NDArray> (tif, image_ov.int16_array_value (),
                                     image_data);
        else
          error ("Only int16 data is allowed for int images with bit depth of 16");
        break;
      case 32:
        if (image_ov.is_int32_type ())
          write_image<int32NDArray> (tif, image_ov.int32_array_value (),
                                     image_data);
        else
          error ("Only int32 data is allowed for int images with bit depth of 32");
        break;
      case 64:
        if (image_ov.is_int64_type ())
          write_image<int64NDArray> (tif, image_ov.int64_array_value (),
                                     image_data);
        else
          error ("Only int64 data is allowed for int images with bit depth of 64");
        break;
      default:
        error ("Unsupported bit depth for signed images");
      }
  }

  void
  write_float_image (TIFF *tif, octave_value image_ov,
                     tiff_image_data *image_data)
  {
    switch (image_data->bits_per_sample)
      {
      case 32:
        if (image_ov.is_single_type () || image_ov.is_double_type ())
          write_image<FloatNDArray> (tif, image_ov.float_array_value (),
                                     image_data);
        else
          error ("Only single or double data are allowed for float images");
        break;
      case 64:
        if (image_ov.is_single_type () || image_ov.is_double_type ())
          write_image<NDArray> (tif, image_ov.array_value (),
                                image_data);
        else
          error ("Only single or double data are allowed for float images");
        break;
      default:
        error ("Unsupported bit depth for floating-point images");
      }
  }

  octave_value_list
  slice_rgba (uint8NDArray rgba_data)
  {
    Array<idx_vector> idx (dim_vector (3, 1));
    idx(0) = idx_vector (':');
    idx(1) = idx_vector (':');
    idx(2) = idx_vector (0, 3);
    uint8NDArray rgb = uint8NDArray (rgba_data.index (idx));
    idx(2) = idx_vector (3);
    uint8NDArray alpha = uint8NDArray (rgba_data.index (idx));
    
    octave_value_list retval (2);
    retval(0) = octave_value (rgb);
    retval(1) = octave_value (alpha);
    return retval;
  }
#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)
      {
        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", "r+"
    };
      
    if (std::find (supported_modes.cbegin (), supported_modes.cend (), mode)
          == supported_modes.cend ())
      error ("Invalid mode for openning Tiff file: %s", mode.c_str ());
    
    // LibTIFF doesn't support r+ mode but it appears that the difference
    // between "a" and "r+" in the context of LibTIFF is that in "a" mode
    // the directory is not set while in r+ the directory is set to the
    // first directory in the file
    bool is_rplus = false;
    if (mode == "r+")
      {
        is_rplus = true;
        mode = "a";
      }
    
    // LibTIFF does Strip chopping by default, which makes the organization of
    // data exposed to the user different from the actual file, so we need to
    // force disable this behavior by adding 'c' flag to the mode
    mode = mode + 'c';
    
    set_internal_handlers ();

    TIFF *tif = TIFFOpen (filename.c_str (), mode.c_str ());
    
    if (! tif)
      error ("Failed to open Tiff file\n");
    
    if (is_rplus && ! TIFFSetDirectory (tif, 0))
      error ("Failed to open Tiff file\n");

    octave_value tiff_ov = octave_value (new octave_tiff_handle (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");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    
    set_internal_handlers ();

    tiff_handle->close ();

    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");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    const TIFFField *fip;
    
    if (args(1).is_string ())
      {
        std::string tag_name = args(1).string_value ();
        fip = get_fip_with_name (tif, tag_name);
        if (! fip)
          error ("Tiff tag not found");
      }
    else
      {
        ttag_t 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");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();
    check_readonly (tif);
    
    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 =  get_fip_with_name (tif, key);
            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;
        if (args(1).is_string ())
          {
            std::string tag_name = args(1).string_value ();
            fip = get_fip_with_name (tif, tag_name);
            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");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);
    
    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    // FIXME: add support for reading in YCbCr mode
    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");

    octave_value_list retval;
    switch (sample_format)
      {
      case 1:
      case 4:
        retval (0) = read_unsigned_image (tif, &image_data);
        break;
      case 2:
        retval (0) = read_signed_image (tif, &image_data);
        break;
      case 3:
        retval (0) = read_float_image (tif, &image_data);
        break;
      default:
        error ("Unsupported sample format");
      }
    
    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 ("Wrong number of arguments");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    // FIXME: add support for reading in YCbCr mode
    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 ("Wrong number of arguments");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    // FIXME: add support for reading in YCbCr mode
    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_read_rgba_image__, args, ,
         "Read the image data in rgba mode")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin != 1)
      error ("Wrong number of arguments");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    tiff_image_data image_data (tif);

    uint16_t orientation;
    if (! TIFFGetFieldDefaulted (tif, TIFFTAG_ORIENTATION, &orientation))
      orientation = ORIENTATION_LEFTTOP;

    // Start with reversed dimensions to be aligned with LibTIFF and
    // permute to the correct order later
    dim_vector img_dims (4, image_data.width, image_data.height);
    uint8NDArray img (img_dims);
    uint32_t *img_ptr = reinterpret_cast <uint32_t *> (img.fortran_vec ());
    
    // The TIFFRGBAImageGet interface is used instead of TIFFReadRGBAImage
    // to have control on the requested orientation of the image
    // Matlab R2022a handles the orientation for each individual strip
    // or tile on its own, this results in a distorted image which is
    // not useful, and is probably a bug. So this deviated from matlab
    // by applying the orientation to the entire image to produce more
    // useful results.
    TIFFRGBAImage img_config;
    char emsg[1024];
    if (! TIFFRGBAImageOK (tif, emsg)
        || ! TIFFRGBAImageBegin (&img_config, tif, 0, emsg))
      error ("Failed to read image");
    
    // FIXME: rotated orientation don't work correctly (e.g. LeftTop)
    img_config.orientation = ORIENTATION_TOPLEFT;
    img_config.req_orientation = orientation;

    bool success = TIFFRGBAImageGet (&img_config, img_ptr, img_config.width,
                                     img_config.height);
                                     
    TIFFRGBAImageEnd (&img_config);
    if (!success)
      error ("Failed to read image");
    
    // Permute to the correct Octave dimension order
    Array<octave_idx_type> perm (dim_vector (3, 1));
    perm(0) = 2;
    perm(1) = 1;
    perm(2) = 0;
    img = img.permute (perm);

    // Slice the data into RGB and alpha
    return slice_rgba (img);
#else
    err_disabled_feature ("readRGBAImage", "Tiff");
#endif
  }

  DEFUN (__tiff_read_rgba_strip__, args, ,
         "Read the strip data containing the given row in rgba mode")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin != 2)
      error ("Wrong number of arguments");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    // Matlab (R2022a) requires row to be double
    if (! args(1).is_double_type ())
      error ("row should be of type double");
    
    uint32_t row = args(1).uint32_scalar_value ();
    
    tiff_image_data image_data (tif);
    if (image_data.is_tiled)
      error ("The image is tiled not stripped");
      
    if (row < 1 || row > image_data.height)
      error ("Row out of bounds of the image");
    
    // Convert from 1-based indexing to zero-based
    row--;

    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;
    
    // LibTIFF requires the row to be the first row in the strip
    // so this removes any offset to reach the first row
    row -= row % rows_in_strip;

    // The exact number of rows in the strip is needed for boundary strips
    uint32_t strip_no = TIFFComputeStrip (tif, row, 0);
    rows_in_strip = get_rows_in_strip (strip_no, TIFFNumberOfStrips (tif),
                                        rows_in_strip, &image_data);
    
    uint16_t orientation;
    if (! TIFFGetFieldDefaulted (tif, TIFFTAG_ORIENTATION, &orientation))
      orientation = ORIENTATION_LEFTTOP;

    // Start with reversed dimensions to be aligned with LibTIFF and
    // permute to the correct order later
    dim_vector strip_dims (4, image_data.width, rows_in_strip);
    uint8NDArray strip_data (strip_dims);
    uint32_t *strip_ptr
      = reinterpret_cast <uint32_t *> (strip_data.fortran_vec ());

    TIFFRGBAImage img_config;
    char emsg[1024];
    if (! TIFFRGBAImageOK (tif, emsg)
        || ! TIFFRGBAImageBegin (&img_config, tif, 0, emsg))
      error ("Failed to read strip");
    
    img_config.orientation = ORIENTATION_TOPLEFT;
    img_config.req_orientation = orientation;
    img_config.row_offset = row;
    img_config.col_offset = 0;

    bool success = TIFFRGBAImageGet (&img_config, strip_ptr, img_config.width,
                                     rows_in_strip);
                                     
    TIFFRGBAImageEnd (&img_config);
    if (!success)
      error ("Failed to read strip");
    
    // Permute to the correct order of dimensions for Octave
    Array<octave_idx_type> perm (dim_vector (3, 1));
    perm(0) = 2;
    perm(1) = 1;
    perm(2) = 0;
    strip_data = strip_data.permute (perm);

    // Slice the data into RGB and alpha
    return slice_rgba (strip_data);
#else
    err_disabled_feature ("readRGBAStrip", "Tiff");
#endif
  }

  DEFUN (__tiff_read_rgba_tile__, args, ,
         "Read the stile data containing the given row and column in rgba mode")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin != 3)
      error ("Wrong number of arguments");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();
    
    TIFF *tif = tiff_handle->get_file ();

    if (! args(1).is_double_type ())
      error ("row should be of type double");
    if (! args(2).is_double_type ())
      error ("col should be of type double");
    
    uint32_t row = args(1).uint32_scalar_value ();
    uint32_t col = args(2).uint32_scalar_value ();
    
    tiff_image_data image_data (tif);
    if (! image_data.is_tiled)
      error ("The image is stripped not tiled");
    
    if (row < 1 || row > image_data.height)
      error ("Row out of bounds of the image");
    if (col < 1 || col > image_data.width)
      error ("Column out of bounds of the image");
    
    // Convert from 1-based indexing to zero-based
    row--;
    col--;

    uint32_t tile_width, tile_height;
    if (! TIFFGetField (tif, TIFFTAG_TILELENGTH, &tile_height))
      error ("Failed to obtain a value for TileLength");
    if (! TIFFGetField (tif, TIFFTAG_TILEWIDTH, &tile_width))
      error ("Failed to obtain a value for TileWidth");
    
    // LibTIFF requires the row and columns to be the top-left corner of the
    // tile, so this removes any offset to reach the top-left row and column
    // of the tile
    row -= row % tile_height;
    col -= col % tile_width;
    
    uint16_t orientation;
    if (! TIFFGetFieldDefaulted (tif, TIFFTAG_ORIENTATION, &orientation))
      orientation = ORIENTATION_LEFTTOP;

    // Calculate the correct dimensions for boundary tiles
    uint32_t corrected_height = tile_height;
    uint32_t corrected_width = tile_width;
    if (row + tile_height > image_data.height)
      corrected_height = image_data.height - row;
    if (col + tile_width > image_data.width)
      corrected_width = image_data.width - col;

    // Start with reversed dimensions to be aligned with LibTIFF and
    // permute to the correct order later
    dim_vector tile_dims (4, corrected_width, corrected_height);
    uint8NDArray tile_data (tile_dims);
    uint32_t *tile_ptr
      = reinterpret_cast <uint32_t *> (tile_data.fortran_vec ());

    TIFFRGBAImage img_config;
    char emsg[1024];
    if (! TIFFRGBAImageOK (tif, emsg)
        || ! TIFFRGBAImageBegin (&img_config, tif, 0, emsg))
      error ("Failed to read tile");
    
    img_config.orientation = ORIENTATION_TOPLEFT;
    img_config.req_orientation = orientation;
    img_config.row_offset = row;
    img_config.col_offset = col;

    bool success = TIFFRGBAImageGet (&img_config, tile_ptr, corrected_width,
                                     corrected_height);
                                     
    TIFFRGBAImageEnd (&img_config);
    if (!success)
      error ("Failed to read tile");
    
    // Permute to the correct order of dimensions for Octave
    Array<octave_idx_type> perm (dim_vector (3, 1));
    perm(0) = 2;
    perm(1) = 1;
    perm(2) = 0;
    tile_data = tile_data.permute (perm);
    
    // Slice the data into RGB and alpha
    return slice_rgba (tile_data);
#else
    err_disabled_feature ("readRGBATile", "Tiff");
#endif
  }

  DEFUN (__tiff_write__, args, ,
             "Write the image data to the current IFD")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    // FIXME: add support for writing in YCbCr mode
    if (nargin < 2)
      error ("Wrong number of arguments\n");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    check_readonly (tif);

    // 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");
    
    switch (sample_format)
      {
      case 1:
      case 4:
        write_unsigned_image (tif, args(1), &image_data);
        break;
      case 2:
        write_signed_image (tif, args(1), &image_data);
        break;
      case 3:
        write_float_image (tif, args(1), &image_data);
        break;
      default:
        error ("Unsupported sample format");
      }
    
    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 ();

    // FIXME: add support for writing in YCbCr mode
    if (nargin < 3)
      error ("Too few arguments provided\n");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    check_readonly (tif);

    // 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 ();

    // FIXME: add support for writing in YCbCr mode
    if (nargin < 3)
      error ("Too few arguments provided\n");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    check_readonly (tif);

    // 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");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();
    
    bool is_tiled = static_cast<bool> (TIFFIsTiled (tif));
    return ovl (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");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();
    
    if (TIFFIsTiled (tif))
      error ("The image is tiled not stripped");
    
    double strip_count = static_cast<double> (TIFFNumberOfStrips (tif));
    return ovl (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");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();
    
    if (! TIFFIsTiled (tif))
      error ("The image is stripped not tiled");
    
    double tile_count = static_cast<double> (TIFFNumberOfTiles (tif));
    return ovl (tile_count);
#else
    err_disabled_feature ("numberOfTiles", "Tiff");
#endif
  }

  DEFUN (__tiff_number_of_directories__, 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");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();
    
    double dir_count = static_cast<double> (TIFFNumberOfDirectories (tif));
    return ovl (dir_count);
#else
    err_disabled_feature ("numberOfDirectories", "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");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();
    
    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 ovl (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");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();
    
    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 ovl (tile_number);
#else
    err_disabled_feature ("computeTile", "Tiff");
#endif
  }

  DEFUN (__tiff_current_directory__, args, ,
         "Get the index of the current directory")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin != 1)
      error ("Wrong number of arguments\n");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    uint16_t dir = TIFFCurrentDirectory (tif);
    dir++;
    
    return ovl (static_cast<double> (dir));
#else
    err_disabled_feature ("currentDirectory", "Tiff");
#endif
  }

  DEFUN (__tiff_last_directory__, args, ,
         "Get the whether the current directory is the last")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin != 1)
      error ("Wrong number of arguments\n");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    bool is_last = TIFFLastDirectory (tif);
    
    return ovl (is_last);
#else
    err_disabled_feature ("lastDirectory", "Tiff");
#endif
  }

  DEFUN (__tiff_next_directory__, args, ,
         "Set the next IFD as the current IFD")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin != 1)
      error ("Wrong number of arguments\n");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    bool is_last = TIFFLastDirectory (tif);
    if (is_last)
      error ("Current directory is the last directory");
    
    if (! TIFFReadDirectory (tif))
      error ("Failed to read the next directory");
    
    return octave_value_list ();
#else
    err_disabled_feature ("nextDirectory", "Tiff");
#endif
  }

  DEFUN (__tiff_set_directory__, args, ,
         "Set the current IFD using the given index")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin != 2)
      error ("Wrong number of arguments\n");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    uint16_t dir = args(1).uint16_scalar_value ();
    if (dir < 1 || dir > TIFFNumberOfDirectories (tif))
      error ("Directory index out of bounds");
    
    dir--;

    if (! TIFFSetDirectory(tif, dir))
      error ("Failed to read directory");

    return octave_value_list ();
#else
    err_disabled_feature ("setDirectory", "Tiff");
#endif
  }

  DEFUN (__tiff_write_directory__, args, ,
         "Write the current IFD to file and create a new one")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin != 1)
      error ("Wrong number of arguments\n");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    TIFF *tif = tiff_handle->get_file ();
    
    if (! TIFFWriteDirectory(tif))
      error ("Failed to write directory");

    return octave_value_list ();
#else
    err_disabled_feature ("writeDirectory", "Tiff");
#endif
  }

  DEFUN (__tiff_rewrite_directory__, args, ,
         "Rewrite modifications to the current IFD")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin != 1)
      error ("Wrong number of arguments\n");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    set_internal_handlers ();

    TIFF *tif = tiff_handle->get_file ();

    if (! TIFFRewriteDirectory(tif))
      error ("Failed to rewrite directory");

    return octave_value_list ();
#else
    err_disabled_feature ("rewriteDirectory", "Tiff");
#endif
  }

  DEFUN (__tiff_set_sub_directory__, args, ,
         "Set the given offset directory as the current IFD")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin != 2)
      error ("Wrong number of arguments\n");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);
    
    set_internal_handlers ();
    
    TIFF *tif = tiff_handle->get_file ();
    
    if (! args(1).is_double_type () && ! args(1).is_uint32_type ()
        && ! args(1).is_uint64_type ())
      error ("Expected offset of type double, uint32 or uint64");
    
    uint16_t count;
    uint64_t *offsets;
    if (! TIFFGetField (tif, TIFFTAG_SUBIFD, &count, &offsets))
      error ("The current directory doesn't have sub directories");
    
    uint64_t offset = args(1).uint64_scalar_value ();

    int found = 0;
    for (uint16_t i = 0; i < count; i++)
      if (offset == offsets[i])
        {
          found = 1;
          break;
        }
    
    if (! found)
      error ("No Sub directory with the given offset");
    
    if (! TIFFSetSubDirectory (tif, offset))
      error ("Failed to switch to the sub directory");

    return octave_value_list ();
#else
    err_disabled_feature ("setSubDirectory", "Tiff");
#endif
  }

  DEFUN (__tiff_version__, , ,
             "Get the version stamp of LibTIFF")
  {
#if defined (HAVE_TIFF)
    std::string version = std::string (TIFFGetVersion ());
    return ovl (version);
#else
    err_disabled_feature ("getVersion", "Tiff");
#endif
  }

  DEFUN (__tiff_set_errors_enabled__, args, ,
         "Enables or Disables error output from LibTIFF")
  {
#if defined (HAVE_TIFF)
    if (args.length () == 0)
      error ("No state argument provided");

    if (! args(0).is_bool_scalar ())
      error ("Expected logical value as argument");
    
    // Set the error and warning handlers according to the bool parameter
    handlers_set = args(0).bool_value ();
    set_internal_handlers ();

    return octave_value_list ();
#else
    err_disabled_feature ("setErrorEnabled", "Tiff");
#endif
  }

  DEFUN (__tiff_make_tagid__, , ,
         "Create the TagID structure from the internal map")
  {
#if defined (HAVE_TIFF)
    std::map<std::string, octave_value> tag_ov_map;
    for (auto it = tag_name_map.cbegin (); it != tag_name_map.cend (); it++)
      tag_ov_map[it->first] = octave_value (it->second);
    return octave_value_list(octave_scalar_map (tag_ov_map));
#else
    err_disabled_feature ("F__tiff_make_tagid__", "Tiff");
#endif
  }

  DEFUN (__tiff_imwrite__, args, ,
         "Handler for imwrite that uses Tiff interface")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin < 2)
      error ("imwrite: Wrong number of arguments");

    if (! (args(0).isnumeric () || args(0).is_bool_matrix ()
           || args(0).is_bool_scalar ()))
      error ("imwrite: Expected image matrix as first argument");

    // The argument checks here are very similar to the ones found in
    // __imwrite__.m and imwrite_filename.m, the code is duplicated here since
    // the original code is tailored to the generic library graphicsmagick
    // so it is not suitable to be used directly for this case
    uint16_t offset = 1;
    std::string filename;
    octave_value cmap = octave_value (NDArray ());
    if (args(1).is_string ())
      {
        filename = args(1).string_value ();
        offset++;
      }
    else if (nargin > 2 && is_colormap (args(1)) && args(2).is_string ())
      {
        filename = args(2).string_value ();
        cmap = args(1);
        offset += 2;
      }
    else
      error ("imwrite: no FILENAME specified");
    
    filename = sys::file_ops::tilde_expand (filename);
    
    if (nargin > offset && (nargin - offset) % 2 != 0
        && args(offset).is_string ())
      offset++;

    if ((nargin - offset) % 2 != 0)
      error ("imwrite: no pair for all arguments (odd number left)");

    // Setting the default value for the rest of the parameters
    bool append = false;
    uint16_t compression = COMPRESSION_NONE;
    uint32_t rows_per_strip = UINT32_MAX;
    float x_res = 72, y_res = 72;
    std::string description = "";

    // Extract the values for the paramters provided
    for (uint16_t arg_idx = offset; arg_idx < nargin; arg_idx+=2)
      {
        if (! args(arg_idx).is_string ())
          error ("imwrite: PARAM in PARAM/VALUE pair must be string");
        
        const char *param_cstr = args(arg_idx).string_value ().c_str ();
        if (strcasecmp (param_cstr, "colorspace") == 0)
          {
            // Matlab specifies three colorspaces: RGB, CIELAB and ICCLAB.
            // Of the three, only RGB is currently supported so this tag
            // can be ignored.
          }
        else if (strcasecmp (param_cstr, "compression") == 0)
          {
            if (! args(arg_idx + 1).is_string ())
              error ("imwrite: value for %s option must be a string",
                     param_cstr);
            std::string comp = args(arg_idx + 1).string_value ();
            if (comp == "packbits")
              compression = COMPRESSION_PACKBITS;
            else if (comp == "lzw")
              compression = COMPRESSION_LZW;
            else if (comp == "deflate")
              compression = COMPRESSION_DEFLATE;
            else if (comp == "jpeg")
              compression = COMPRESSION_JPEG;
            else if (comp == "ccitt")
              compression = COMPRESSION_CCITTRLE;
            else if (comp == "fax3")
              compression = COMPRESSION_CCITTFAX3;
            else if (comp == "fax4")
              compression = COMPRESSION_CCITTFAX4;
            else
              error ("imwrite: invalid compression '%s'", comp.c_str ());
          }
        else if (strcasecmp (param_cstr, "Description") == 0)
          {
            if (! args(arg_idx + 1).is_string ())
              error ("imwrite: value for %s option must be a string",
                     param_cstr);
            description = args(arg_idx + 1).string_value ();
          }
        else if (strcasecmp (param_cstr, "resolution") == 0)
          {
            if (! args(arg_idx + 1).isnumeric ())
              error ("imwrite: value for %s option must be numeric",
                     param_cstr);
            NDArray res = args(arg_idx + 1).array_value ();
            if (res.numel () == 1)
              {
                x_res = res(0);
                y_res = res(0);
              }
            else if (res.numel () == 2)
              {
                x_res = res(0);
                y_res = res(1);
              }
            else
              error ("imwrite: value for %s option must be either a scalar value or a two-element vector",
                     param_cstr);
          }
        else if (strcasecmp (param_cstr, "RowsPerStrip") == 0)
          {
            if (! args(arg_idx + 1).isnumeric ()
                || ! args(arg_idx + 1).is_scalar_type ())
              error ("imwrite: value for %s option must be a scalar value",
                     param_cstr);
            rows_per_strip = args(arg_idx + 1).uint32_scalar_value ();
          }
        else if (strcasecmp (param_cstr, "WriteMode") == 0)
          {
            if (! args(arg_idx + 1).is_string ())
              error ("imwrite: value for %s option must be a string",
                     param_cstr);
            std::string wmode = args(arg_idx + 1).string_value ();
            if (wmode == "overwrite")
              append = false;
            else if (wmode == "append")
              append = true;
            else
              error ("imwrite: value for %s option must be either 'overwrite' or 'append'",
                     param_cstr);
          }
        else
          error ("imread: invalid PARAMETER '%s'", param_cstr);
      }
    
    // Matlab specifies that JPEG compression must also specify RowsPerStrip
    // and must be a value divisible by 8
    if (compression == COMPRESSION_JPEG
        && (rows_per_strip == UINT32_MAX || rows_per_strip % 8 != 0))
      error ("imwrite: RowsPerStrip option must be specified for jpeg compression and must be divisible by 8");
    
    // These compression schemes are only available for binary images
    if (! (args(0).is_bool_matrix () || args(0).is_bool_scalar ())
        && (compression == COMPRESSION_CCITTRLE
            || compression == COMPRESSION_CCITTFAX3
            || compression == COMPRESSION_CCITTFAX4))
      error ("imwrite: the specified compression scheme is for binary images only");

    // Matlab rejects 4D data for TIFF images
    dim_vector img_dims = args(0).dims ();
    if (args(0).ndims () != 2 && args(0).ndims () != 3)
      error ("imwrite: expected 2-dimensional or 3-dimensional image matrix");

    if (args(0).ndims () > 2 && img_dims(2) > 1 && cmap.numel () > 0)
      error ("imwrite: palette images must be 1 channel only");

    TIFF *tif;
    if (append)
      tif = TIFFOpen (filename.c_str (), "a");
    else
      tif = TIFFOpen (filename.c_str (), "w");
    
    if (! tif)
      error ("Failed to open file %s", filename.c_str ());
    
    // A simple way to make sure the file will be closed when the function
    // returns or when an error occurs as the destructor will always be called
    octave_tiff_handle tiff_handle (tif);
    
    // Set all necessary tags
    if (! TIFFSetField (tif, TIFFTAG_IMAGELENGTH,
                        static_cast<uint32_t>(img_dims(0)))
        || ! TIFFSetField (tif, TIFFTAG_IMAGEWIDTH,
                           static_cast<uint32_t>(img_dims(1))))
      error ("Failed to set image dimensions");
  
    if (img_dims.ndims () > 2)
      if (! TIFFSetField (tif, TIFFTAG_SAMPLESPERPIXEL, img_dims(2)))
        error ("Failed to set the number of samples");
    
    if (! TIFFSetField (tif, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG))
      error ("Failed to set the planar configuration");
    
    if (! TIFFSetField (tif, TIFFTAG_COMPRESSION, compression))
      error ("Failed to set the compressoin tag");
    
    if (! TIFFSetField (tif, TIFFTAG_XRESOLUTION, x_res)
        || ! TIFFSetField (tif, TIFFTAG_YRESOLUTION, y_res))
      error ("Failed to set resolution tag");
    
    if (! TIFFSetField (tif, TIFFTAG_IMAGEDESCRIPTION, description.c_str ()))
      error ("Failed to set description tag");

    if (rows_per_strip == UINT32_MAX)
      rows_per_strip = TIFFDefaultStripSize (tif, 0);
    
    if (! TIFFSetField (tif, TIFFTAG_ROWSPERSTRIP, rows_per_strip))
      error ("Failed to set the RowsPerStrip tag");
    
    uint16_t bits_per_sample;
    if (args(0).is_bool_matrix () || args(0).is_bool_scalar ())
      bits_per_sample = 1;
    else if (args(0).is_uint8_type ())
      bits_per_sample = 8;
    else if (args(0).is_uint16_type ())
      bits_per_sample = 16;
    else if (args(0).is_uint32_type () || args(0).is_single_type ())
      bits_per_sample = 32;
    else if (args(0).is_double_type ())
      bits_per_sample = 64;
    else
      error ("imwrite: unsupported image data type");
    
    if (! TIFFSetField (tif, TIFFTAG_BITSPERSAMPLE, bits_per_sample))
      error ("Failed to set BitsPerSample tag");
    
    // Infer the photometric interpretation of the image from the color map
    // and the number of channels of the input
    uint16_t photometric = PHOTOMETRIC_MINISBLACK;
    if (cmap.numel () > 0)
      {
        photometric = PHOTOMETRIC_PALETTE;
        set_field_data (tif, TIFFFieldWithTag (tif, TIFFTAG_COLORMAP),
                        cmap);
      }
    else if (img_dims(2) == 3)
      photometric = PHOTOMETRIC_RGB;
    else if (img_dims(2) == 4)
      photometric = PHOTOMETRIC_SEPARATED;

    if (! TIFFSetField (tif, TIFFTAG_PHOTOMETRIC, photometric))
      error ("Failed to set the photometric tag");
    
    tiff_image_data image_data (tif);
    if (args(0).is_bool_scalar () || args(0).is_bool_matrix ()
        || args(0).is_uint8_type () || args(0).is_uint16_type ()
        || args(0).is_uint32_type ())
      write_unsigned_image (tif, args(0), &image_data);
    else
      {
        if (!TIFFSetField (tif, TIFFTAG_SAMPLEFORMAT, SAMPLEFORMAT_IEEEFP))
          error ("Failed to set SampleFormat tag");
        write_float_image (tif, args(0).as_double (), &image_data);
      }
    
    return octave_value_list ();
#else
    err_disabled_feature ("imwrite", "Tiff");
#endif
  }

  DEFUN (__tiff_imfinfo__, args, ,
         "Handler for imfinfo that uses Tiff interface")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin < 1)
      error ("imfinfo: missing filename argument");
    
    if (! args(0).is_string ())
      error ("imfinfo: filename must be a string");
    
    std::string filename = args(0).string_value ();

    const sys::file_stat fs (filename);
    if (! fs)
      error ("imfinfo: error reading '%s': %s", filename.c_str (),
             fs.error ().c_str ());
    
    TIFF *tif = TIFFOpen (filename.c_str (), "rc");
    if (! tif)
      error ("imfinfo: error reading '%s': LibTIFF failed to read file",
             filename.c_str ());
    
    // The destructor for this object will be called when the function returnd
    // or in case of an error so the file will always get closed at the end
    octave_tiff_handle tiff_handle (tif);

    // A lot of the logic here is copied from __magick_finfo__ due to the
    // great similarity between the two functions but this function is
    // format specific so a lot of the details are different
    uint16_t dir_count = TIFFNumberOfDirectories (tif);
    
    // Null terminated char* list to be used to create a string_vector
    static const char *fields[] = {
      "Filename",
      "FileModDate",
      "FileSize",
      "Format",
      "FormatVersion",
      "Width",
      "Height",
      "BitDepth",
      "ColorType",
      "FormatSignature",
      "ByteOrder",
      "NewSubFileType",
      "BitsPerSample",
      "Compression",
      "PhotometricInterpretation",
      "StripOffsets",
      "SamplesPerPixel",
      "RowsPerStrip",
      "StripByteCounts",
      "XResolution",
      "YResolution",
      "ResolutionUnit",
      "Colormap",
      "PlanarConfiguration",
      "TileWidth",
      "TileLength",
      "TileOffsets",
      "TileByteCounts",
      "Orientation",
      "FillOrder",
      "GrayResponseUnit",
      "MaxSampleValue",
      "MinSampleValue",
      "Thresholding",
      "Offset",
      "ImageDescription",
      nullptr
    };
    
    // A map to be used as a struct array to hold a scalar_map for each
    // directory
    octave_map info (dim_vector (dir_count, 1), string_vector (fields));

    // populate template_info with the info that is common between all
    // directories in the file
    octave_scalar_map template_info = (string_vector (fields));
    template_info.setfield ("Format", octave_value ("tif"));
    template_info.setfield ("FormatVersion", octave_value (""));
    const sys::localtime mtime (fs.mtime ());
    const std::string filetime = mtime.strftime ("%e-%b-%Y %H:%M:%S");
    template_info.setfield ("Filename",    octave_value (filename));
    template_info.setfield ("FileModDate", octave_value (filetime));
    template_info.setfield ("FileSize",    octave_value (fs.size ()));

    // Extract the image signature (first 4 bytes in file)
    std::ifstream tif_file;
#if defined (OCTAVE_USE_WINDOWS_API)
    std::wstring wname = sys::u8_to_wstring (filename);
    tif_file.open (wname.c_str (), std::ios_base::binary);
#else
    tif_file.open (filename.c_str (), std::ios_base::binary);
#endif
    uint8NDArray signature (dim_vector (1, 4));;
    tif_file.read (reinterpret_cast<char *> (signature.fortran_vec ()), 4);
    tif_file.close ();
    template_info.setfield ("FormatSignature", octave_value (signature));

    std::string byte_order
      = TIFFIsBigEndian (tif)? "big-endian": "little-endian";
    template_info.setfield ("ByteOrder", octave_value (byte_order));

    // Extract directory specific information
    for (uint16_t dir = 0; dir < dir_count; dir++)
      {
        octave_scalar_map dir_info (template_info);

        // Switch to the directory
        if (! TIFFSetDirectory (tif, dir))
          error ("imfinfo: Failed to access frame %d\n", dir);
        
        tiff_image_data image_data (tif);
        dir_info.setfield ("Width", octave_value (image_data.width));
        
        dir_info.setfield ("Height", octave_value (image_data.height));
        
        uint16_t bit_depth = image_data.samples_per_pixel
                             * image_data.bits_per_sample;
        dir_info.setfield ("BitDepth", octave_value (bit_depth));
        
        std::string planar = "unrecognized";
        if (image_data.planar_configuration == 1)
          planar = "Chunky";
        else if (image_data.planar_configuration == 2)
          planar = "Separate";
        dir_info.setfield ("PlanarConfiguration", octave_value (planar));

        // Extract photometric information as well as color map if exists
        std::string color_str, photometric_str;
        uint16_t photometric = PHOTOMETRIC_MINISBLACK;
        octave_value cmap = octave_value (Matrix ());
        const TIFFField *fip;
        TIFFGetField (tif, TIFFTAG_PHOTOMETRIC, &photometric);
        switch (photometric)
          {
          case PHOTOMETRIC_MINISBLACK:
            color_str = "grayscale";
            photometric_str = "BlasckIsZero";
            break;
          case PHOTOMETRIC_MINISWHITE:
            color_str = "grayscale";
            photometric_str = "WhiteIsZero";
            break;
          case PHOTOMETRIC_RGB:
            color_str = "truecolor";
            photometric_str = "RGB";
            break;
          case PHOTOMETRIC_PALETTE:
            color_str = "indexed";
            photometric_str = "RGB Palette";
            fip = TIFFFieldWithTag (tif, TIFFTAG_COLORMAP);
            cmap = get_field_data (tif, fip);
            break;
          case PHOTOMETRIC_SEPARATED:
            color_str = "CMYK";
            photometric_str = "CMYK";
            break;
          default:
            color_str = "undefined";
            photometric_str = "undefined";
          }
        dir_info.setfield ("ColorType", octave_value(color_str));
        dir_info.setfield ("PhotometricInterpretation",
                           octave_value(photometric_str));
        dir_info.setfield ("Colormap", octave_value (cmap));

        fip = TIFFFieldWithTag (tif, TIFFTAG_SUBFILETYPE);
        dir_info.setfield ("NewSubFileType", get_field_data (tif, fip));

        fip = TIFFFieldWithTag (tif, TIFFTAG_BITSPERSAMPLE);
        dir_info.setfield ("BitsPerSample", get_field_data (tif, fip));

        fip = TIFFFieldWithTag (tif, TIFFTAG_SAMPLESPERPIXEL);
        dir_info.setfield ("SamplesPerPixel", get_field_data (tif, fip));

        // Use LibTIFF's compression codec to extract compression scheme name
        uint16_t compression;
        TIFFGetFieldDefaulted (tif, TIFFTAG_COMPRESSION, &compression);
        std::string comp_str = "unrecognized";
        const TIFFCodec *codec = TIFFFindCODEC (compression);
        if (codec)
          comp_str = codec->name;
        dir_info.setfield ("Compression", octave_value (comp_str));

        // Set strip-specific and tile-specific fields accordingly
        bool tiled = TIFFIsTiled (tif);
        fip = TIFFFieldWithTag (tif, TIFFTAG_TILELENGTH);
        dir_info.setfield ("TileLength", tiled? get_field_data (tif, fip)
                                         : octave_value (Matrix ()));
        fip = TIFFFieldWithTag (tif, TIFFTAG_TILEWIDTH);
        dir_info.setfield ("TileWidth", tiled? get_field_data (tif, fip)
                                        : octave_value (Matrix ()));
        fip = TIFFFieldWithTag (tif, TIFFTAG_TILEOFFSETS);
        dir_info.setfield ("TileOffsets", tiled? get_field_data (tif, fip)
                                          : octave_value (Matrix ()));
        fip = TIFFFieldWithTag (tif, TIFFTAG_TILEBYTECOUNTS);
        dir_info.setfield ("TileByteCounts", tiled? get_field_data (tif, fip)
                                             : octave_value (Matrix ()));
        fip = TIFFFieldWithTag (tif, TIFFTAG_ROWSPERSTRIP);
        dir_info.setfield ("RowsPerStrip", tiled? octave_value (Matrix ())
                                           : get_field_data (tif, fip));
        fip = TIFFFieldWithTag (tif, TIFFTAG_STRIPOFFSETS);
        dir_info.setfield ("StripOffsets", tiled? octave_value (Matrix ())
                                           : get_field_data (tif, fip));
        fip = TIFFFieldWithTag (tif, TIFFTAG_STRIPBYTECOUNTS);
        dir_info.setfield ("StripByteCounts", tiled? octave_value (Matrix ())
                                              : get_field_data (tif, fip));
        
        uint16_t res;
        if (TIFFGetField (tif, TIFFTAG_XRESOLUTION, &res))
          dir_info.setfield ("XResolution", octave_value (res));
        else
          dir_info.setfield ("XResolution", octave_value (Matrix ()));

        if (TIFFGetField (tif, TIFFTAG_YRESOLUTION, &res))
          dir_info.setfield ("YResolution", octave_value (res));
        else
          dir_info.setfield ("YResolution", octave_value (Matrix ()));
        
        TIFFGetFieldDefaulted (tif, TIFFTAG_RESOLUTIONUNIT, &res);
        std::string res_unit = "Inch";
        if (res == 1)
          res_unit = "None";
        else if (res == 3)
          res_unit = "Centimeter";
        dir_info.setfield ("ResolutionUnit", octave_value(res_unit));

        fip = TIFFFieldWithTag (tif, TIFFTAG_ORIENTATION);
        dir_info.setfield ("Orientation", get_field_data (tif, fip));

        fip = TIFFFieldWithTag (tif, TIFFTAG_FILLORDER);
        dir_info.setfield ("FillOrder", get_field_data (tif, fip));
        
        // The current version of LibTIFF (4.4.0) doesn't set the deafult
        // value for GrayResponseUnit corectly, so we can't use
        // TIFFGetFieldDefaulted, instead we set the default value ourselves
        double gray_response_unit = 0.01;
        uint16_t gray_unit_val;
        if (TIFFGetField (tif, TIFFTAG_GRAYRESPONSEUNIT, &gray_unit_val))
          {
            switch (gray_unit_val)
              {
              case GRAYRESPONSEUNIT_10S:
                gray_response_unit = 0.1;
                break;
              case GRAYRESPONSEUNIT_100S:
                gray_response_unit = 0.01;
                break;
              case GRAYRESPONSEUNIT_1000S:
                gray_response_unit = 0.001;
                break;
              case GRAYRESPONSEUNIT_10000S:
                gray_response_unit = 0.0001;
                break;
              case GRAYRESPONSEUNIT_100000S:
                gray_response_unit = 0.00001;
                break;
              }
          }
        dir_info.setfield ("GrayResponseUnit",
                           octave_value (gray_response_unit));

        // The current version of LibTIFF (4.4.0) doesn't set the deafult
        // value for MinSampleValue and MaxSampleValue, so we can't use
        // TIFFGetFieldDefaulted, instead we set the default value ourselves
        uint16_t min_sample_value = 0;
        uint16_t max_sample_value = (1<<image_data.bits_per_sample) - 1;
        TIFFGetField (tif, TIFFTAG_MINSAMPLEVALUE, &min_sample_value);
        TIFFGetField (tif, TIFFTAG_MAXSAMPLEVALUE, &max_sample_value);
        dim_vector vector_dims = dim_vector (1, image_data.samples_per_pixel);
        NDArray min_sample_values (vector_dims, min_sample_value);
        NDArray max_sample_values (vector_dims, max_sample_value);
        dir_info.setfield ("MinSampleValue",
                           octave_value (min_sample_values));
        dir_info.setfield ("MaxSampleValue",
                           octave_value (max_sample_values));

        fip = TIFFFieldWithTag (tif, TIFFTAG_THRESHHOLDING);
        dir_info.setfield ("Thresholding", get_field_data (tif, fip));

        dir_info.setfield ("Offset",
                           octave_value (TIFFCurrentDirOffset (tif)));
        
        char *desc = NULL;
        if (TIFFGetField (tif, TIFFTAG_IMAGEDESCRIPTION, &desc))
          dir_info.setfield ("ImageDescription", octave_value (desc));
        else
          dir_info.setfield ("ImageDescription", octave_value (""));

        // Insert the directory information into the map
        info.fast_elem_insert (dir, dir_info);
      }

    return ovl (info);
#else
    err_disabled_feature ("imfinfo", "Tiff");
#endif
  }
}