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

#include <string>
#include <iostream>

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

#include "errwarn.h"

#include "fcntl-wrappers.h"

#if defined (HAVE_TIFF)
#  include <tiffio.h>
TIFFErrorHandler tiff_default_error_handler = NULL;
TIFFErrorHandler tiff_default_warning_handler = NULL;
#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
        planar_configuration = 1;
      
      is_tiled = TIFFIsTiled(tif);
    }
  };

  // A map of tag names supported by matlab, there are some differences
  // than LibTIFF's names (e.g. Photometric vs PhotometricInerpretation)
  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 = TIFFFieldWithName (tif, tag_name.c_str ());
    
    // If the tag name is not found, fallback to the names defined
    // in our name map.
    if (! fip && tag_name_map.find (tag_name) != tag_name_map.cend ())
      fip = TIFFFieldWithTag (tif, tag_name_map.at (tag_name));
    
    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);

    // TODO(maged): find/create images to test the special tags
    switch (tag_id)
      {
      case TIFFTAG_STRIPBYTECOUNTS:
      case TIFFTAG_STRIPOFFSETS:
        tag_data_ov = get_array_field_data (tif, fip, 
                                            TIFFNumberOfStrips (tif));
        break;
      case TIFFTAG_TILEBYTECOUNTS:
      case TIFFTAG_TILEOFFSETS:
        tag_data_ov
          = get_array_field_data (tif, fip, TIFFNumberOfTiles (tif));
        break;
      case TIFFTAG_YCBCRCOEFFICIENTS:
        tag_data_ov = get_array_field_data (tif, fip, 3, 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;
        }
      // TODO(maged): These tags are more complex to implement
      //  will be implemented and tested later.
      case TIFFTAG_XMLPACKET:
      case TIFFTAG_RICHTIFFIPTC:
      case TIFFTAG_PHOTOSHOP:
      case TIFFTAG_ICCPROFILE:
        {
          error ("Complex Tags not implemented");
          break;
        }
      // These tags are not mentioned in the LibTIFF documentation
      // but are handled correctly by the library
      case TIFFTAG_ZIPQUALITY:
      case TIFFTAG_SGILOGDATAFMT:
      case TIFFTAG_GRAYRESPONSEUNIT:
        {
          tag_data_ov = get_scalar_field_data (tif, fip);
          break;
        }
      default:
        tag_data_ov = get_scalar_field_data (tif, fip);
    }
    
    return tag_data_ov;
  }

  void
  set_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);
    
    // TODO(maged): find/create images to test the special tags
    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 ();
          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;
        }
      // TODO(maged): These tags are more complex to implement
      //  will be implemented and tested later.
      case TIFFTAG_XMLPACKET:
      case TIFFTAG_RICHTIFFIPTC:
      case TIFFTAG_PHOTOSHOP:
      case TIFFTAG_ICCPROFILE:
        {
          error ("Complex Tags not implemented");
          break;
        }
      // These tags are not mentioned in the LibTIFF documentation
      // but are handled correctly by the library
      case TIFFTAG_ZIPQUALITY:
      case TIFFTAG_SGILOGDATAFMT:
      case TIFFTAG_GRAYRESPONSEUNIT:
        {
          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))
              error ("Failed to rite strip data");
            pixel_fvec += strip_size;
          }
        else if (image_data->bits_per_sample == 1)
          {
            if (image_data->samples_per_pixel != 1)
              error ("Bi-Level images must have one channel only");
            
            // Create a buffer to hold the packed strip data
            // Unique pointers are faster than vectors for constant size buffers
            std::unique_ptr<uint8_t []> strip_ptr
              = std::make_unique<uint8_t []> (TIFFStripSize (tif));
            uint8_t *strip_buf = strip_ptr.get ();
            // According to the format specification, the row should be byte
            // aligned so the number of bytes is rounded up to the nearest byte
            uint32_t padded_width = (image_data->width + 7) / 8;
            // Packing the pixel data into bits
            for (uint32_t row = 0; row < rows_in_strip; row++)
              {
                for (uint32_t col = 0; col < image_data->width; col++)
                {
                  uint8_t shift = 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))
              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';
    
    TIFF *tif = TIFFOpen (filename.c_str (), mode.c_str ());
    
    if (! tif)
      error ("Failed to open Tiff file\n");
    
    // TODO(maged): check of there is more to be done for r+
    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));
    
    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);

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

    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);
    
    TIFF *tif = tiff_handle->get_file ();

    // TODO(maged): nargout and ycbcr
    octave_unused_parameter (nargout);

    // Obtain all necessary tags
    // The main purpose of this struct is to hold all the necessary tags that
    // will be common among all the different cases of handling the the image
    // data to avoid repeating the same calls to TIFFGetField in the different
    // functions as each call is a possible point of failure
    tiff_image_data image_data (tif);

    uint16_t sample_format;
    if (! TIFFGetFieldDefaulted(tif, TIFFTAG_SAMPLEFORMAT, &sample_format))
      error ("Failed to obtain a value for sample format");

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

    TIFF *tif = tiff_handle->get_file ();

    // TODO(maged): nargout and ycbcr
    octave_unused_parameter (nargout);

    if (TIFFIsTiled (tif))
      error ("The image is tiled not stripped");
    
    uint32_t strip_no;
    if (args(1).is_scalar_type ())
      strip_no = args(1).uint32_scalar_value ();
    else
      error ("Expected scalar for strip number");
    
    if (strip_no < 1 || strip_no > TIFFNumberOfStrips (tif))
      error ("Strip number out of bounds");
    
    // Convert from Octave's 1-based indexing to zero-based indexing
    strip_no--;

    return octave_value_list (handle_read_strip_or_tile (tif, strip_no));
#else
    err_disabled_feature ("readEncodedStrip", "Tiff");
#endif
  }

  DEFUN (__tiff_read_encoded_tile__, args, nargout,
             "Read a tile from the image in the current IFD")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

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

    TIFF *tif = tiff_handle->get_file ();

    // TODO(maged): nargout and ycbcr
    octave_unused_parameter (nargout);

    if (! TIFFIsTiled (tif))
      error ("The image is stripped not tiled");
    
    uint32_t tile_no;
    if (args(1).is_scalar_type ())
      tile_no = args(1).uint32_scalar_value ();
    else
      error ("Expected scalar for tile number");
    
    if (tile_no < 1 || tile_no > TIFFNumberOfTiles (tif))
      error ("Tile number out of bounds");
    
    // Convert from Octave's 1-based indexing to zero-based indexing
    tile_no--;

    return octave_value_list (handle_read_strip_or_tile (tif, tile_no));
#else
    err_disabled_feature ("readEncodedTile", "Tiff");
#endif
  }

  DEFUN (__tiff_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);

    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
    TIFFRGBAImage img_config;
    char emsg[1024];
    if (! TIFFRGBAImageOK (tif, emsg)
        || ! TIFFRGBAImageBegin (&img_config, tif, 0, emsg))
      error ("Failed to read image");
    
    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);

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

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

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

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

    // TODO(maged): add support for YCbCr data
    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);

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

    // TODO(maged): add support for YCbCr data
    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);

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

    TIFF *tif = tiff_handle->get_file ();
    
    bool is_tiled = static_cast<bool> (TIFFIsTiled (tif));
    return octave_value_list (octave_value (is_tiled));
#else
    err_disabled_feature ("isTiled", "Tiff");
#endif
  }

  DEFUN (__tiff_number_of_strips__, args, ,
             "Get the number of strips in the image")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin == 0)
      error ("No handle provided\n");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    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 octave_value_list (octave_value (strip_count));
#else
    err_disabled_feature ("numberOfStrips", "Tiff");
#endif
  }

  DEFUN (__tiff_number_of_tiles__, args, ,
             "Get the number of tiles in the image")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin == 0)
      error ("No handle provided\n");
    
    octave_tiff_handle *tiff_handle
      = octave_tiff_handle::get_tiff_handle (args(0));
    check_closed (tiff_handle);

    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 octave_value_list (octave_value (tile_count));
#else
    err_disabled_feature ("numberOfTiles", "Tiff");
#endif
  }

  DEFUN (__tiff_compute_strip__, args, ,
             "Get the strip index containing the given row")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

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

    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 octave_value_list (octave_value (strip_number));
#else
    err_disabled_feature ("computeStrip", "Tiff");
#endif
  }

  DEFUN (__tiff_compute_tile__, args, ,
             "Get the tile index containing the given row and column")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

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

    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 octave_value_list (octave_value (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);

    TIFF *tif = tiff_handle->get_file ();

    uint16_t dir = TIFFCurrentDirectory (tif);
    dir++;
    
    return octave_value_list (octave_value (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);

    TIFF *tif = tiff_handle->get_file ();

    bool is_last = TIFFLastDirectory (tif);
    
    return octave_value_list (octave_value (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);

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

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

    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);
    
    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 octave_value_list (octave_value (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)
    // Get the default error handlers the first time this function is called
    if (tiff_default_error_handler == NULL)
      {
        tiff_default_error_handler = TIFFSetErrorHandler (NULL);
        tiff_default_warning_handler = TIFFSetWarningHandler (NULL);
      }

    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
    if (args(0).bool_value ())
      {
        TIFFSetErrorHandler (tiff_default_error_handler);
        TIFFSetWarningHandler (tiff_default_warning_handler);
      }
    else
      {
        TIFFSetErrorHandler (NULL);
        TIFFSetWarningHandler (NULL);
      }

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

  DEFUN (__tiff_make_tagid__, , ,
         "Create the TagID structure from the internal map")
  {
    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));
  }
}