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

  struct tiff_image_data
  {
  public:
    uint32_t width;
    uint32_t height;
    uint16_t samples_per_pixel;
    uint16_t bits_per_sample;
    uint16_t planar_configuration;
    uint16_t is_tiled;

    tiff_image_data (TIFF *tif)
    {
      if (! TIFFGetField (tif, TIFFTAG_IMAGEWIDTH, &width))
        error ("Failed to read image width");

      if (! TIFFGetField (tif, TIFFTAG_IMAGELENGTH, &height))
        error ("Failed to read image height");
      
      if (! TIFFGetFieldDefaulted (tif, TIFFTAG_SAMPLESPERPIXEL,
                                   &samples_per_pixel))
        error ("Failed to read the SamplesPerPixel tag");

      if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
                                   &bits_per_sample))
        error ("Failed to read the BitsPerSample tag");
      
      // TODO(maged): this doesn't really work for writing as LibTIFF will
      // refuse to write unless the tag is set
      if (! TIFFGetField (tif, TIFFTAG_PLANARCONFIG,
                          &planar_configuration))
        // LibTIFF has a bug where it incorrectly returns 0 as a default
        // value for PlanarConfiguration although the default value is 1
        planar_configuration = 1;
      
      is_tiled = TIFFIsTiled(tif);
    }
  };

  void
  check_readonly (TIFF *tif)
  {
    // This can't use O_RDONLY directly because its header file "fcntl.h"
    // isn't available on Windows. The wrapper, however, seems to return
    // the right thing even on Windows.
    if (TIFFGetMode (tif) == octave_o_rdonly_wrapper ())
      error ("Can't write data to a file opened in read-only mode");
  }

  // Error if status is not 1 (success status for TIFFGetField)
  void
  validate_tiff_get_field (bool status, void *p_to_free=NULL)
  {
      if (status != 1)
        {
          if (p_to_free != NULL)
            _TIFFfree (p_to_free);
          error ("Failed to read tag");
        }
  }

  uint32_t get_rows_in_strip (uint32_t strip_no, uint32_t strip_count,
                              uint32_t rows_per_strip,
                              tiff_image_data *image_data)
  {
    // Calculate the expected number of elements in the strip data array
    // All strips have equal number of rows except strips at the bottom
    // of the image can have less number of rows
    uint32_t rows_in_strip = rows_per_strip;
    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
      {
        // All strips have equal number of rows excpet strips at the bottom
        // of the image can have less number of rows
        if (strip_no == strip_count - 1)
          rows_in_strip = image_data->height - rows_in_strip * strip_no;
      }
    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      {
        // For separate planes, we should check the last strip of each channel
        uint32_t strips_per_channel
          = strip_count / image_data->samples_per_pixel;
        for (uint32_t boundary_strip = strips_per_channel - 1;
             boundary_strip <= strip_count;
             boundary_strip += strips_per_channel)
          if (strip_no == boundary_strip)
            rows_in_strip = image_data->height
                            - rows_in_strip * (strip_no % strips_per_channel);
      }
    else
      error ("Planar Configuration not supported");
    
    return rows_in_strip;
  }

  template <typename T>
  octave_value
  read_strip (TIFF *tif, uint32_t strip_no, tiff_image_data * image_data)
  {
    // ASSUMES stripped image and strip_no is a valid zero-based strip
    // index for the tif image
    
    uint32_t rows_in_strip;
    if (! TIFFGetFieldDefaulted (tif, TIFFTAG_ROWSPERSTRIP, &rows_in_strip))
      error ("Failed to obtain a value for RowsPerStrip");
    
    if (rows_in_strip > image_data->height)
      rows_in_strip = image_data->height;
    
    uint32_t strip_count = TIFFNumberOfStrips (tif);
    rows_in_strip = get_rows_in_strip (strip_no, strip_count,
                                       rows_in_strip, image_data);
    dim_vector strip_dims;
    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
      strip_dims = dim_vector (image_data->samples_per_pixel,
                               image_data->width, rows_in_strip);
    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      strip_dims = dim_vector (image_data->width, rows_in_strip, 1);
    else
      error ("Unsupported bit depth");
    
    T strip_data (strip_dims);
    uint8_t *strip_fvec
      = reinterpret_cast<uint8_t *> (strip_data.fortran_vec ());

    if (image_data->bits_per_sample == 8
        || image_data->bits_per_sample == 16
        || image_data->bits_per_sample == 32
        || image_data->bits_per_sample == 64)
      {
        if (TIFFReadEncodedStrip (tif, strip_no, strip_fvec, -1) == -1)
          error ("Failed to read strip data");
      }
    else if (image_data->bits_per_sample == 1)
      {
        if (image_data->samples_per_pixel != 1)
          error ("Bi-Level images must have one channel only");
        
        // Create a buffer to hold the packed strip data
        // Unique pointers are faster than vectors for constant size buffers
        std::unique_ptr<uint8_t []> strip_ptr
          = std::make_unique<uint8_t []> (TIFFStripSize (tif));
        uint8_t *strip_buf = strip_ptr.get ();
        if (TIFFReadEncodedStrip (tif, strip_no, strip_buf, -1) == -1)
          error ("Failed to read strip data");
        
        // According to the format specification, the row should be byte
        // aligned so the number of bytes is rounded up to the nearest byte
        uint32_t padded_width = (image_data->width + 7) / 8;
        // Packing the pixel data into bits
        for (uint32_t row = 0; row < rows_in_strip; row++)
          {
            for (uint32_t col = 0; col < image_data->width; col++)
            {
              uint8_t shift = 7 - col % 8;
              strip_fvec[col] = (strip_buf[col / 8] >> shift) & 0x1;
            }
            strip_fvec += image_data->width;
            strip_buf += padded_width;
          }
      }
    else
      error ("Unsupported bit depth");
    
    Array<octave_idx_type> perm (dim_vector (3, 1));
    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
      {
        perm(0) = 2;
        perm(1) = 1;
        perm(2) = 0;
      }
    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      {
        perm(0) = 1;
        perm(1) = 0;
        perm(2) = 2;
      }
    
    strip_data = strip_data.permute (perm);
    return octave_value (strip_data);
  }

  template <typename T>
  octave_value
  read_tile (TIFF *tif, uint32_t tile_no, tiff_image_data *image_data)
  {
    // ASSUMES tiled image and tile_no is a valid zero-based tile
    // index for the tif image

    // TODO(maged): refactor into a function?
    uint32_t tile_width, tile_height;
    if (! TIFFGetField (tif, TIFFTAG_TILELENGTH, &tile_height))
      error ("Filed to read tile length");
    
    if (! TIFFGetField (tif, TIFFTAG_TILEWIDTH, &tile_width))
      error ("Filed to read tile length");

    if (tile_height == 0 || tile_height % 16 != 0
        || tile_width == 0 || tile_width % 16 != 0)
      error ("Tile dimesion tags are invalid");

    dim_vector tile_dims;
    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
      {
        tile_dims = dim_vector (image_data->samples_per_pixel, tile_width,
                                tile_height);
      }
    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      {
        tile_dims = dim_vector (tile_width, tile_height, 1);
      }
    else
      error ("Unsupported planar configuration");
    
    T tile_data (tile_dims);
    uint8_t *tile_fvec
      = reinterpret_cast<uint8_t *> (tile_data.fortran_vec ());

    if (image_data->bits_per_sample == 8
        || image_data->bits_per_sample == 16
        || image_data->bits_per_sample == 32
        || image_data->bits_per_sample == 64)
        {      
          if (TIFFReadEncodedTile (tif, tile_no, tile_fvec, -1) == -1)
            error ("Failed to read tile data");
        }
    else if (image_data->bits_per_sample == 1)
      {
        if (image_data->samples_per_pixel != 1)
          error ("Bi-Level images must have one channel only");
        
        // Create a buffer to hold the packed tile data
        // Unique pointers are faster than vectors for constant size buffers
        std::unique_ptr<uint8_t []> tile_ptr
          = std::make_unique<uint8_t []> (TIFFTileSize (tif));
        uint8_t *tile_buf = tile_ptr.get ();
        if (TIFFReadEncodedTile (tif, tile_no, tile_buf, -1) == -1)
            error ("Failed to read tile data");
        
        // unpack tile bits into output matrix cells
        for (uint32_t row = 0; row < tile_height; row++)
          {
            for (uint32_t col = 0; col < tile_width; col++)
              {
                uint8_t shift = 7 - col % 8;
                tile_fvec[col] = (tile_buf [col / 8] >> shift) & 0x1;
              }  
            tile_fvec += tile_width;
            tile_buf += tile_width / 8;
          }
      }
    else
      error ("Unsupported bit depth");
    
    Array<octave_idx_type> perm (dim_vector (3, 1));
    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
      {
        perm(0) = 2;
        perm(1) = 1;
        perm(2) = 0;
      }
    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      {
        perm(0) = 1;
        perm(1) = 0;
        perm(2) = 2;
      }
    
    tile_data = tile_data.permute (perm);

    // Get the actual tile dimensions
    uint32_t tiles_across = (image_data->width + tile_width - 1)
                            / tile_width;
    uint32_t tiles_down = (image_data->height + tile_height - 1)
                          / tile_height;
    uint32_t corrected_width = tile_width;
    uint32_t corrected_height = tile_height;
    if (tile_no % tiles_across == tiles_across - 1)
      corrected_width = image_data->width - tile_width * (tiles_across - 1);
    if ((tile_no / tiles_across) % tiles_down == tiles_down - 1)
      corrected_height = image_data->height - tile_height * (tiles_down - 1);
    
    dim_vector corrected_dims;
    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
      corrected_dims = dim_vector (corrected_height, corrected_width,
                                   image_data->samples_per_pixel);
    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      corrected_dims = dim_vector (corrected_height, corrected_width);
    
    tile_data.resize (corrected_dims);

    return octave_value (tile_data);
  }

  template <typename T>
  octave_value
  read_strip_or_tile (TIFF *tif, uint32_t strip_tile_no,
                      tiff_image_data *image_data)
  {
    if (image_data->is_tiled)
      return read_tile<T> (tif, strip_tile_no, image_data);
    else
      return read_strip<T> (tif, strip_tile_no, image_data);
  }

  octave_value
  handle_read_strip_or_tile (TIFF *tif, uint32_t strip_tile_no)
  {
    // Obtain all necessary tags
    tiff_image_data image_data (tif);

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

    if (sample_format == 3)
      {
        if (image_data.bits_per_sample != 32
            && image_data.bits_per_sample != 64)
          error ("Floating point images are only supported for bit depths of 32 and 64");
      }
    else if (sample_format != 1 && sample_format != 4)
      error ("Unsupported sample format");
    
    switch (image_data.bits_per_sample)
      {
      case 1:
        return read_strip_or_tile<boolNDArray> (tif, strip_tile_no,
                                                &image_data);
        break;
      case 8:
        return read_strip_or_tile<uint8NDArray> (tif, strip_tile_no,
                                                 &image_data);
        break;
      case 16:
        return read_strip_or_tile<uint16NDArray> (tif, strip_tile_no,
                                                  &image_data);
        break;
      case 32:
        if (sample_format == 3)
          return read_strip_or_tile<FloatNDArray> (tif, strip_tile_no,
                                                   &image_data);
        else
          return read_strip_or_tile<uint32NDArray> (tif, strip_tile_no,
                                                    &image_data);
        break;
      case 64:
        if (sample_format == 3)
          return read_strip_or_tile<NDArray> (tif, strip_tile_no,
                                              &image_data);
        else
          return read_strip_or_tile<uint64NDArray> (tif, strip_tile_no,
                                                    &image_data);
        break;
      default:
        error ("Unsupported bit depth");
      }
  }

  template <typename T>
  octave_value
  read_stripped_image (TIFF *tif, tiff_image_data *image_data)
  {
    typedef typename T::element_type P;

    // For 1-bit and 4-bit images, each row must be byte aligned and padding
    // is added to the end of the row to reach the byte mark.
    // To facilitate reading the data, the matrix is defined with the padded
    // size and the padding is removed at the end.
    uint32_t padded_width = image_data->width;
    uint8_t remove_padding = 0;
    if ((image_data->bits_per_sample == 1 || image_data->bits_per_sample == 4)
        && padded_width % 8 != 0)
      {
        padded_width += (8 - padded_width % 8) % 8;
        remove_padding = 1;
      }
    
    // The matrix dimensions are defined in the order that corresponds to
    // the order of strip data read from LibTIFF.
    // At the end, the matrix is permutated to the order expected by Octave
    T img;
    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
      img = T (dim_vector (image_data->samples_per_pixel, padded_width,
                           image_data->height));
    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      img = T (dim_vector (padded_width, image_data->height,
                           image_data->samples_per_pixel));
    else
      error ("Unsupported Planar Configuration");

    P *img_fvec = img.fortran_vec ();

    // Obtain the necessary data for handling the strips
    uint32_t strip_count = TIFFNumberOfStrips (tif);
    
    // Can't rely on StripByteCounts because in compressed images
    // the byte count reflect the actual number of bytes stored
    // in the file not the size of the uncompressed strip
    int64_t strip_size;
    uint64_t written_size = 0;
    uint64_t image_size = padded_width * image_data->height
                          * image_data->samples_per_pixel
                          * sizeof (P);
    for (uint32_t strip = 0; strip < strip_count; strip++)
      {
        // Read the strip data into the matrix directly
        strip_size = TIFFReadEncodedStrip (tif, strip, img_fvec, -1);

        // Check if the strip read failed.
        if (strip_size == -1)
          error ("Failed to read strip data");
        
        // Check if the size being read exceeds the bounds of the matrix
        // In case of a corrupt image with more data than needed
        // This is probably redundant as LibTIFF checks sizes internally
        if (written_size + strip_size > image_size)
          error ("Strip data is larger than the image size");

        if (image_data->bits_per_sample == 1)
          {
            if (image_data->samples_per_pixel != 1)
              error ("Bi-Level images must have one channel only");
            
            // The strip size is multiplied by 8 to reflect tha actual
            // number of bytes written to the matrix since each byte
            // in the original strip contains 8 pixels of data
            strip_size *= 8;

            // Checking bounds again with the new size
            if (written_size + strip_size > image_size)
              error ("Strip data is larger than the image size");

            // Iterate over the memory region backwards to expand the bits
            // to their respective bytes without overwriting the read data
            for (int64_t pixel = strip_size - 1; pixel >= 0; pixel--)
              {
                // TODO(maged): is it necessary to check FillOrder?
                uint8_t bit_number = 7 - pixel % 8;
                uint8_t * img_u8 = reinterpret_cast<uint8_t *> (img_fvec);
                img_fvec[pixel] = (img_u8[pixel / 8] >> bit_number) & 0x01;
              }
          }
        else if (image_data->bits_per_sample == 4)
          {
            if (image_data->samples_per_pixel != 1)
              error ("4-bit images are only supported for grayscale");
            
            // Strip size is multplied by as each byte contains 2 pixels
            // and each pixel is later expanded into its own byte
            strip_size *= 2;

            // Checking bounds again with the ne size
            if (written_size + strip_size > image_size)
              error ("Strip data is larger than the image size");
            
            // Iterate over the memory region backwards to expand the nibbles
            // to their respective bytes without overwriting the read data
            for (int64_t pixel = strip_size - 1; pixel >= 0; pixel--)
              {
                uint8_t shift = pixel % 2 == 0? 4: 0;
                uint8_t * img_u8 = reinterpret_cast<uint8_t *> (img_fvec);
                img_fvec[pixel] = (img_u8[pixel / 2] >> shift) & 0x0F;
              }
          }
        else if (image_data->bits_per_sample != 8 &&
                 image_data->bits_per_sample != 16 &&
                 image_data->bits_per_sample != 32 &&
                 image_data->bits_per_sample != 64)
          error ("Unsupported bit depth");

        // Advance the pointer by the amount of bytes read
        img_fvec 
          = reinterpret_cast<P *> (reinterpret_cast <uint8_t *> (img_fvec)
                                   + strip_size);
        written_size += strip_size;
      }

    // The matrices are permutated back to the shape expected by Octave
    // which is height x width x channels
    Array<octave_idx_type> perm (dim_vector (3, 1));
    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
      {
        perm(0) = 2;
        perm(1) = 1;
        perm(2) = 0;
      }
    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      {
        perm(0) = 1;
        perm(1) = 0;
        perm(2) = 2;
      }
    
    img = img.permute (perm);

    if (remove_padding)
      img.resize (dim_vector (image_data->height, image_data->width));

    return octave_value (img);
  }

  template <typename T>
  octave_value
  read_tiled_image (TIFF *tif, tiff_image_data *image_data)
  {
    typedef typename T::element_type P;

    // Obtain the necessary data for handling the tiles
    uint32_t tile_width, tile_height;
    validate_tiff_get_field (TIFFGetField (tif, TIFFTAG_TILEWIDTH,
                                           &tile_width));
    validate_tiff_get_field (TIFFGetField (tif, TIFFTAG_TILELENGTH,
                                           &tile_height));
    uint32_t tile_count = TIFFNumberOfTiles (tif);
    uint32_t tiles_across = (image_data->width + tile_width - 1)
                            / tile_width;
    uint32_t tiles_down = (image_data->height + tile_height - 1)
                          / tile_height;

    T img;
    // The matrix dimensions are defined in the order that corresponds to
    // the order of the tile data read from LibTIFF.
    // At the end, the matrix is permutated, reshaped and resized to be in the
    // shape expected by Octave
    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
      img = T (dim_vector (image_data->samples_per_pixel, tile_width,
                           tile_height, tiles_across, tiles_down));
    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      img = T (dim_vector (tile_width, tile_height, tiles_across, tiles_down,
                           image_data->samples_per_pixel));
    else
      error ("Unsupported Planar Configuration");

    P *img_fvec = img.fortran_vec ();

    // image_size is calculated from the tile data and not from the
    // image parameters to account for the additional padding in the
    // boundary tiles
    uint64_t image_size = tile_width * tile_height * tile_count * sizeof(P);
    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
      image_size *= image_data->samples_per_pixel;

    // Can't rely on TileByteCounts because compressed images will report
    // the number of bytes present in the file as opposed to the actual
    // number of bytes of uncompressed data that is needed here
    int64_t tile_size;    
    uint64_t written_size = 0;
    for (uint32_t tile = 0; tile < tile_count; tile++)
      {
        tile_size = TIFFReadEncodedTile(tif, tile, img_fvec, -1);
        
        if (tile_size == -1)
          error ("Failed to read tile data");

        // Checking if the read bytes would exceed the size of the matrix
        if (tile_size + written_size > image_size)
          error ("Tile data is larger than image size");
        
        if (image_data->bits_per_sample == 1)
          {
            if (image_data->samples_per_pixel != 1)
              error ("Bi-Level images must have one channel only");
            
            // The tile size is multiplied by 8 to reflect tha actual
            // number of bytes written to the matrix since each byte
            // in the original tile contains 8 pixels of data
            tile_size *= 8;

            // Checking bounds again with the new size
            if (written_size + tile_size > image_size)
              error ("Tile data is larger than the image size");

            // Iterate over the memory region backwards to expand the bits
            // to their respective bytes without overwriting the read data
            for (int64_t pixel = tile_size - 1; pixel >= 0; pixel--)
              {
                // TODO(maged): is it necessary to check FillOrder?
                uint8_t bit_number = 7 - pixel % 8;
                uint8_t * img_u8 = reinterpret_cast<uint8_t *> (img_fvec);
                img_fvec[pixel]= (img_u8[pixel / 8] >> bit_number) & 0x01;
              }
          }
        else if (image_data->bits_per_sample == 4)
          {
            if (image_data->samples_per_pixel != 1)
              error ("4-bit images are only supported for grayscale");
            
            // tile size is multplied by as each byte contains 2 pixels
            // and each pixel is later expanded into its own byte
            tile_size *= 2;

            // Checking bounds again with the ne size
            if (written_size + tile_size > image_size)
              error ("Tile data is larger than the image size");
            
            // Iterate over the memory region backwards to expand the nibbles
            // to their respective bytes without overwriting the read data
            for (int64_t pixel = tile_size - 1; pixel >= 0; pixel--)
              {
                uint8_t shift = pixel % 2 == 0? 4: 0;
                uint8_t * img_u8 = reinterpret_cast<uint8_t *> (img_fvec);
                img_fvec[pixel] = (img_u8[pixel / 2] >> shift) & 0x0F;
              }
          }
        else if (image_data->bits_per_sample != 8 &&
                 image_data->bits_per_sample != 16 &&
                 image_data->bits_per_sample != 32 &&
                 image_data->bits_per_sample != 64)
          error ("Unsupported bit depth");
        
        img_fvec
          = reinterpret_cast<P *> (reinterpret_cast<uint8_t *> (img_fvec)
                                   + tile_size);
        written_size += tile_size;
      }
    
    // The data is now in the matrix but in a different order than expected
    // by Octave and with additional padding in boundary tiles.
    // To get it to the right order, the dimensions are permutated to
    // align tiles to their correct grid, then reshaped to remove the
    // extra dimensions (tiles_across, tiles_down), then resized to
    // remove any extra padding, and finally permutated to the correct
    // order that is: height x width x channels
    Array<octave_idx_type> perm1 (dim_vector (5, 1));
    Array<octave_idx_type> perm2 (dim_vector (3, 1));
    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
      {
        perm1(0) = 0;
        perm1(1) = 1;
        perm1(2) = 3;
        perm1(3) = 2;
        perm1(4) = 4;
        img = img.permute (perm1);

        img = img.reshape (dim_vector (image_data->samples_per_pixel,
                                       tile_width * tiles_across,
                                       tile_height * tiles_down));
        
        if (tile_width * tiles_across != image_data->width
            || tile_height * tiles_down != image_data->height)
          img.resize (dim_vector (image_data->samples_per_pixel,
                                  image_data->width, image_data->height));
        
        perm2(0) = 2;
        perm2(1) = 1;
        perm2(2) = 0;
        img = img.permute (perm2);
      }
    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      {
        perm1(0) = 0;
        perm1(1) = 2;
        perm1(2) = 1;
        perm1(3) = 3;
        perm1(4) = 4;
        img = img.permute (perm1);

        img = img.reshape (dim_vector (tile_width * tiles_across,
                                       tile_height * tiles_down,
                                       image_data->samples_per_pixel));
        
        if (tile_width * tiles_across != image_data->width
            || tile_height * tiles_down != image_data->height)
          img.resize (dim_vector (image_data->width, image_data->height,
                                  image_data->samples_per_pixel));
        
        perm2(0) = 1;
        perm2(1) = 0;
        perm2(2) = 2;
        img = img.permute (perm2);
      }

    return octave_value (img);
  }

  template <typename T>
  octave_value
  read_image (TIFF *tif, tiff_image_data *image_data)
  {
    if (image_data->is_tiled)
      return read_tiled_image<T> (tif, image_data);
    else
      return read_stripped_image<T> (tif, image_data);
  }

  // Convert tag value to double
  octave_value
  interpret_scalar_tag_data (void *data, TIFFDataType tag_datatype)
  {
    double retval;

    switch (tag_datatype)
      {
      case TIFF_BYTE:
      case TIFF_UNDEFINED:
        {
          retval = static_cast<double> (*(reinterpret_cast<uint8_t *> (data)));
          break;
        }
      case TIFF_SHORT:
        {
          retval = static_cast<double> (*(reinterpret_cast<uint16_t *> (data)));
          break;
        }
      case TIFF_LONG:
        {
          retval = static_cast<double> (*(reinterpret_cast<uint32_t *> (data)));
          break;
        }
      case TIFF_LONG8:
        {
          retval = static_cast<double> (*(reinterpret_cast<uint64_t *> (data)));
          break;
        }
      case TIFF_RATIONAL:
        {
          retval = *(reinterpret_cast<float *> (data));
          break;
        }
      case TIFF_SBYTE:
        {
          retval = static_cast<double> (*(reinterpret_cast<int8_t *> (data)));
          break;
        }
      case TIFF_SSHORT:
        {
          retval = static_cast<double> (*(reinterpret_cast<int16_t *> (data)));
          break;
        }
      case TIFF_SLONG:
        {
          retval = static_cast<double> (*(reinterpret_cast<int32_t *> (data)));
          break;
        }
      case TIFF_SLONG8:
        {
          retval = static_cast<double> (*(reinterpret_cast<int64_t *> (data)));
          break;
        }
      case TIFF_FLOAT:
        {
          retval = *(reinterpret_cast<float *> (data));
          break;
        }
      case TIFF_DOUBLE:
        {
          retval = *(reinterpret_cast<double *> (data));
          break;
        }
      case TIFF_SRATIONAL:
        {
          retval = *(reinterpret_cast<float *> (data));
          break;
        }
      case TIFF_IFD:
      case TIFF_IFD8:
        error ("Unimplemented IFFD data type");
        break;
      default:
        error ("Unsupported tag data type");
      }

    return octave_value (retval);
  }

  // Convert memory buffer into a suitable octave value
  // depending on tag_datatype
  octave_value
  interpret_tag_data (void *data, uint32_t count, TIFFDataType tag_datatype)
  {
    octave_value retval;
    // Apparently matlab converts scalar numerical values into double
    // but doesn't do the same for arrays
    if (count == 1 && tag_datatype != TIFF_ASCII)
      {
        retval = interpret_scalar_tag_data (data, tag_datatype);
      }
    else
      {
        dim_vector arr_dims (1, count);

        switch (tag_datatype)
          {
          case TIFF_BYTE:
          case TIFF_UNDEFINED:
            {
              uint8NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = (reinterpret_cast<uint8_t *> (data))[i];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_ASCII:
            {
              retval = octave_value (*(reinterpret_cast<char **> (data)));
              break;
            }
          case TIFF_SHORT:
            {
              uint16NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = (reinterpret_cast<uint16_t *> (data))[i];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_LONG:
            {
              uint32NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = (reinterpret_cast<uint32_t *> (data))[i];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_LONG8:
            {
              uint64NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = (reinterpret_cast<uint64_t *> (data))[i];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_RATIONAL:
            {
              NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = (reinterpret_cast<float *> (data))[i];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_SBYTE:
            {
              int8NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = (reinterpret_cast<int8_t *> (data))[i];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_SSHORT:
            {
              int16NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = (reinterpret_cast<int16_t *> (data))[i];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_SLONG:
            {
              int32NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = (reinterpret_cast<int32_t *> (data))[i];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_SLONG8:
            {
              int64NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = (reinterpret_cast<int64_t *> (data))[i];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_FLOAT:
            {
              NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = (reinterpret_cast<float *> (data))[i];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_DOUBLE:
            {
              NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
              {
                  arr(i) = (reinterpret_cast<double *> (data))[i];
              }
              retval = octave_value (arr);
              break;
            }
          case TIFF_SRATIONAL:
            {
              NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = (reinterpret_cast<float *> (data))[i];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_IFD:
          case TIFF_IFD8:
            // TODO(maged): implement IFD datatype?
            error ("Unimplemented IFD data type");
            break;
          default:
            error ("Unsupported tag data type");
          }
      }

    return retval;
  }

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

    // TIFFFieldReadCount returns VARIABLE for some scalar tags
    // (e.g. Compression) But TIFFFieldPassCount seems consistent
    // Since scalar tags are the last to be handled, any tag that
    // require a count to be passed is an unsupported tag.
    if (TIFFFieldPassCount (fip))
      error ("Unsupported tag");
    
    int type_size = TIFFDataWidth (TIFFFieldDataType (fip));
    
    void *data = _TIFFmalloc (type_size);
    if (tag_id == TIFFTAG_PLANARCONFIG)
      {
        // Workaround for a bug in LibTIFF where it incorrectly returns
        // zero as the default value for PlanaConfiguration
        if (! TIFFGetField(tif, tag_id, data))
          *reinterpret_cast<uint16_t *> (data) = 1;
      }
    else
      validate_tiff_get_field (TIFFGetFieldDefaulted (tif, tag_id, data), data);
    
    octave_value tag_data_ov = interpret_tag_data (data, 1,
                                                   TIFFFieldDataType (fip));
    _TIFFfree (data);

    return tag_data_ov;
  }

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

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

    // TODO(maged): find/create images to test the special tags
    switch (tag_id)
      {
      case TIFFTAG_STRIPBYTECOUNTS:
      case TIFFTAG_STRIPOFFSETS:
        tag_data_ov = get_array_field_data (tif, fip, 
                                            TIFFNumberOfStrips (tif));
        break;
      case TIFFTAG_TILEBYTECOUNTS:
      case TIFFTAG_TILEOFFSETS:
        tag_data_ov = get_array_field_data (tif, fip, TIFFNumberOfTiles (tif));
        break;
      case TIFFTAG_YCBCRCOEFFICIENTS:
        tag_data_ov = get_array_field_data (tif, fip, 3);
        break;
      case TIFFTAG_REFERENCEBLACKWHITE:
        tag_data_ov = get_array_field_data (tif, fip, 6);
        break;
      case TIFFTAG_GRAYRESPONSECURVE:
        {
          uint16_t bits_per_sample;
          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
                                            &bits_per_sample))
            error ("Failed to obtain the bit depth");
          
          tag_data_ov = get_array_field_data (tif, fip, 1<<bits_per_sample);
          break;
        }
      case TIFFTAG_COLORMAP:
        {
          uint16_t bits_per_sample;
          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
                                       &bits_per_sample))
            error ("Failed to obtain the bit depth");
          
          // According to the format specification, this field should
          // be 8 or 16 only.
          if (bits_per_sample > 16)
            error ("Too high bit depth for a palette image");

          uint32_t count = 1 << bits_per_sample;
          uint16_t *red, *green, *blue;
          validate_tiff_get_field (TIFFGetField (tif, TIFFTAG_COLORMAP,
                                                 &red, &green, &blue));
          
          // Retrieving the data of the three channels and concatenating
          // them together
          OCTAVE_LOCAL_BUFFER (NDArray, array_list, 3);
          dim_vector col_dims(count, 1);
          array_list[0] = NDArray (interpret_tag_data (red,
                                                       count,
                                                       TIFFFieldDataType(fip))
                                                       .uint16_array_value ()
                                                       .reshape (col_dims));
          array_list[1] = NDArray (interpret_tag_data (green,
                                                       count,
                                                       TIFFFieldDataType(fip))
                                                       .uint16_array_value ()
                                                       .reshape (col_dims));
          array_list[2] = NDArray (interpret_tag_data (blue,
                                                       count,
                                                       TIFFFieldDataType(fip))
                                                       .uint16_array_value ()
                                                       .reshape (col_dims));

          NDArray mat_out = NDArray::cat(1, 3, array_list);
          // Normalize the range to be between 0 and 1
          mat_out /= UINT16_MAX;

          tag_data_ov = octave_value (mat_out);
          break;
        }
      case TIFFTAG_TRANSFERFUNCTION:
        {
          uint16_t samples_per_pixel;
          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_SAMPLESPERPIXEL,
                                       &samples_per_pixel))
            error ("Failed to obtain the number of samples per pixel");

          uint16_t bits_per_sample;
          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
                                       &bits_per_sample))
            error ("Failed to obtain the bit depth");

          uint32_t count = 1 << bits_per_sample;
          uint16_t *ch1, *ch2, *ch3;
          validate_tiff_get_field (TIFFGetField (tif, tag_id,
                                                 &ch1, &ch2, &ch3));
          if (ch2 == NULL)
            tag_data_ov = interpret_tag_data (ch1, count,
                                              TIFFFieldDataType (fip));
          else
            { 
              OCTAVE_LOCAL_BUFFER (uint16NDArray, array_list, 3);
              dim_vector col_dims(count, 1);
              array_list[0] = interpret_tag_data (ch1,
                                                  count,
                                                  TIFFFieldDataType (fip))
                                                  .uint16_array_value ()
                                                  .reshape (col_dims);
              array_list[1] = interpret_tag_data (ch2,
                                                  count,
                                                  TIFFFieldDataType (fip))
                                                  .uint16_array_value ()
                                                  .reshape (col_dims);
              array_list[2] = interpret_tag_data (ch3,
                                                  count,
                                                  TIFFFieldDataType (fip))
                                                  .uint16_array_value ()
                                                  .reshape (col_dims);
              
              tag_data_ov = octave_value (uint16NDArray::cat (1, 3, array_list));
            }
          break;
        }
      case TIFFTAG_PAGENUMBER:
      case TIFFTAG_HALFTONEHINTS:
      case TIFFTAG_DOTRANGE:
      case TIFFTAG_YCBCRSUBSAMPLING:
        {
          uint16_t tag_part1, tag_part2;
          validate_tiff_get_field (TIFFGetField (tif, tag_id,
                                                 &tag_part1, &tag_part2));
          
          NDArray mat_out (dim_vector (1, 2));
          mat_out(0)
            = interpret_tag_data (&tag_part1, 1,
                                  TIFFFieldDataType (fip)).double_value ();
          mat_out(1)
            = interpret_tag_data (&tag_part2, 1,
                                  TIFFFieldDataType (fip)).double_value ();
          
          tag_data_ov = octave_value (mat_out);
          break;
        }
      case TIFFTAG_SUBIFD:
        {
          uint16_t count;
          uint64_t *offsets;
          validate_tiff_get_field (TIFFGetField (tif, tag_id, &count, &offsets));
          tag_data_ov = interpret_tag_data (offsets, count,
                                            TIFFFieldDataType (fip));
          break;
        }
      case TIFFTAG_EXTRASAMPLES:
        {
          uint16_t count;
          uint16_t *types;
          validate_tiff_get_field (TIFFGetField (tif, tag_id, &count, &types));
          tag_data_ov = interpret_tag_data (types, count,
                                            TIFFFieldDataType (fip));
          break;
        }
      // TODO(maged): These tags are more complex to implement
      //  will be implemented and tested later.
      case TIFFTAG_XMLPACKET:
      case TIFFTAG_RICHTIFFIPTC:
      case TIFFTAG_PHOTOSHOP:
      case TIFFTAG_ICCPROFILE:
        {
          error ("Complex Tags not implemented");
          break;
        }
      // These tags are not mentioned in the LibTIFF documentation
      // but are handled correctly by the library
      case TIFFTAG_ZIPQUALITY:
      case TIFFTAG_SGILOGDATAFMT:
      case TIFFTAG_GRAYRESPONSEUNIT:
        {
          tag_data_ov = get_scalar_field_data (tif, fip);
          break;
        }
      default:
        tag_data_ov = get_scalar_field_data (tif, fip);
    }
    
    return tag_data_ov;
  }

  void
  set_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");
    
    // According to matlab, the value must be a scalar double
    if (! tag_ov.is_scalar_type () || ! tag_ov.is_double_type ())
      error ("Tag value must be a scalar double");
    
    TIFFDataType tag_datatype = TIFFFieldDataType (fip);
    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:
        {
          // TODO(maged): Matlab doesnt error on setting this tag
          // but returns 0 for getTag
          error ("Unsupported 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 = 7 - col % 8;
              strip_buf[row * padded_width + col / 8] |= data_u8[col] << shift;
            }
            data_u8 += image_data->width;
          }
        tsize_t strip_size = padded_width * rows_in_strip;
        if (TIFFWriteEncodedStrip (tif, strip_no, strip_buf, strip_size) == -1)
          error ("Failed to write strip data to image");
      }
    else
      {
        error ("Unsupported bit depth");
      }
  }

  template <typename T>
  void
  write_tile (TIFF *tif, uint32_t tile_no, T tile_data,
              tiff_image_data *image_data)
  {
    uint32_t tile_width, tile_height;
    if (! TIFFGetField (tif, TIFFTAG_TILEWIDTH, &tile_width))
      error ("Failed to get the tile width");
    if (! TIFFGetField (tif, TIFFTAG_TILELENGTH, &tile_height))
      error ("Failed to get the tile length");

    if (tile_height == 0 || tile_height % 16 != 0
        || tile_width == 0 || tile_width % 16 != 0)
      error ("Tile dimesion tags are invalid");
    
    if (tile_no < 1 || tile_no > TIFFNumberOfTiles (tif))
      error ("Tile number out of bounds");

    if (tile_data.dim1 () > tile_height)
      warning ("The tile is composed of %u rows but input has %ld rows",
               tile_height, tile_data.dim1 ());
    if (tile_data.dim2 () > tile_width)
      warning ("The tile is composed of %u columns but input has %ld columns",
               tile_width, tile_data.dim2 ());
    if (tile_data.ndims () > 2)
      {
        if (image_data->planar_configuration == PLANARCONFIG_CONTIG
            && tile_data.dim3 () > image_data->samples_per_pixel)
          warning ("The tile is composed of %u channels but input has %ld channels",
                  image_data->samples_per_pixel, tile_data.dim3 ());
        else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE
                 && tile_data.dim3 () > 1)
          warning ("The tile is composed of %u channels but input has %ld channels",
                   1, tile_data.dim3 ());
      }
    
    dim_vector tile_dimensions;
    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
      tile_dimensions = dim_vector (tile_height, tile_width,
                                    image_data->samples_per_pixel);
    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      tile_dimensions = dim_vector (tile_height, tile_width, 1);
    else
      error ("Planar configuration not supported");
    
    tile_data.resize (tile_dimensions);
    Array<octave_idx_type> perm (dim_vector (3, 1));
    perm(0) = 2;
    perm(1) = 1;
    perm(2) = 0;
    tile_data = tile_data.permute (perm);
    
    // Octave indexing is 1-based while LibTIFF is zero-based
    tile_no--;
    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 = 7 - col % 8;
              tile_buf[row * tile_width/8 + col/8] |= data_u8[col] << shift;
            }
            data_u8 += tile_width;
          }
        if (TIFFWriteEncodedTile (tif, tile_no, tile_buf,
                                  TIFFTileSize (tif)) == -1)
          error ("Failed to write tile data to image");
      }
    else
      {
        error ("Unsupported bit depth");
      }
  }

  template <typename T>
  void
  write_strip_or_tile (TIFF *tif, uint32_t strip_tile_no, T strip_data,
                       tiff_image_data *image_data)
  {
    if (image_data->is_tiled)
      write_tile<T> (tif, strip_tile_no, strip_data, image_data);
    else
      write_strip<T> (tif, strip_tile_no, strip_data, image_data);
  }

  void
  handle_write_strip_or_tile (TIFF *tif, uint32_t strip_tile_no,
                              octave_value data_ov,
                              tiff_image_data *image_data)
  {

    // SampleFormat tag is not a required field and has a default value of 1
    // So we need to use TIFFGetFieldDefaulted in case it is not present in
    // the file
    uint16_t sample_format;
    if (! TIFFGetFieldDefaulted(tif, TIFFTAG_SAMPLEFORMAT, &sample_format))
      error ("Failed to obtain a value for sample format");

    // TODO(maged): add support for signed integer images
    if (sample_format == 3)
      {
        if (image_data->bits_per_sample != 32
            && image_data->bits_per_sample != 64)
          error ("Floating point images are only supported for bit depths of 32 and 64");
      }

    // The standard specifies that a SampleFormat of 4 should be treated
    // the same as 1 (unsigned integer)
    else if (sample_format != 1 && sample_format != 4)
      error ("Unsupported sample format");
    
    switch (image_data->bits_per_sample)
      {
      case 1:
        // We need to check for both scalar and matrix types to handle single
        // element strip
        if (data_ov.is_bool_scalar () || data_ov.is_bool_matrix ())
          write_strip_or_tile<boolNDArray> (tif, strip_tile_no,
                                            data_ov.bool_array_value (),
                                            image_data);
        else
          error ("Expected logical matrix for BiLevel image");
        break;
      case 8:
        if (data_ov.is_uint8_type ())
          write_strip_or_tile<uint8NDArray> (tif, strip_tile_no,
                                             data_ov.uint8_array_value (),
                                             image_data);
        else
          error ("Only uint8 data is allowed for uint images with bit depth of 8");
        break;
      case 16:
        if (data_ov.is_uint16_type ())
          write_strip_or_tile<uint16NDArray> (tif, strip_tile_no,
                                              data_ov.uint16_array_value (),
                                              image_data);
        else
          error ("Only uint16 data is allowed for uint images with bit depth of 16");
        break;
      case 32:
        if (sample_format == 3)
          if (data_ov.is_single_type () || data_ov.is_double_type ())
            write_strip_or_tile<FloatNDArray> (tif, strip_tile_no,
                                               data_ov.float_array_value (),
                                               image_data);
          else
            error ("Only single and double data are allowed for floating-point images");
        else
          if (data_ov.is_uint32_type ())
            write_strip_or_tile<uint32NDArray> (tif, strip_tile_no,
                                                data_ov.uint32_array_value (),
                                                image_data);
          else
            error ("Only uint32 data is allowed for uint images with bit depth of 32");
        break;
      case 64:
        if (sample_format == 3)
          if (data_ov.is_single_type () || data_ov.is_double_type ())
            write_strip_or_tile<NDArray> (tif, strip_tile_no,
                                          data_ov.array_value (),
                                          image_data);
          else
            error ("Only single and double data are allowed for floating-point images");
        else  
          if (data_ov.is_uint64_type ())
            write_strip_or_tile<uint64NDArray> (tif, strip_tile_no,
                                                data_ov.uint64_array_value (),
                                                image_data);
          else
            error ("Only uint64 data is allowed for uint images with bit depth of 64");
        break;
      default:
        error ("Unsupported bit depth");
      }
  }

  template <typename T>
  void
  write_stripped_image (TIFF *tif, T pixel_data, tiff_image_data *image_data)
  {
    // ASSUMES pixel data dimensions are already validated

    typedef typename T::element_type P;

    // Permute pixel data to be aligned in memory to the way LibTIFF
    // expects the data to be (i.e. channel x width x height for chunky
    // and width x height x channel for separate planes)
    Array<octave_idx_type> perm (dim_vector (3, 1));
    if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      {
        perm(0) = 1;
        perm(1) = 0;
        perm(2) = 2;
      }
    else
      {
        perm(0) = 2;
        perm(1) = 1;
        perm(2) = 0;
      }
    pixel_data = pixel_data.permute (perm);

    uint32_t row_per_strip;
    if (! TIFFGetFieldDefaulted (tif, TIFFTAG_ROWSPERSTRIP, &row_per_strip))
      error ("Failed to obtain the RowPerStrip tag");
    
    // The default value is INT_MAX so we need to cap it to the image height
    if (row_per_strip > image_data->height)
      row_per_strip = image_data->height;

    uint8_t *pixel_fvec
      = reinterpret_cast<uint8_t *> (pixel_data.fortran_vec ());
    uint32_t strip_count = TIFFNumberOfStrips (tif);
    tsize_t strip_size;
    uint32_t rows_in_strip;
    for (uint32_t strip = 0; strip < strip_count; strip++)
      {
        rows_in_strip = get_rows_in_strip (strip, strip_count,
                                           row_per_strip, image_data);
        if (image_data->bits_per_sample == 8
            || image_data->bits_per_sample == 16
            || image_data->bits_per_sample == 32
            || image_data->bits_per_sample == 64)
          {
            strip_size = rows_in_strip * image_data->width * sizeof (P);
            if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
              strip_size *= image_data->samples_per_pixel;
            if (! TIFFWriteEncodedStrip (tif, strip, pixel_fvec, strip_size))
              error ("Failed to rite strip data");
            pixel_fvec += strip_size;
          }
        else if (image_data->bits_per_sample == 1)
          {
            if (image_data->samples_per_pixel != 1)
              error ("Bi-Level images must have one channel only");
            
            // Create a buffer to hold the packed strip data
            // Unique pointers are faster than vectors for constant size buffers
            std::unique_ptr<uint8_t []> strip_ptr
              = std::make_unique<uint8_t []> (TIFFStripSize (tif));
            uint8_t *strip_buf = strip_ptr.get ();
            // According to the format specification, the row should be byte
            // aligned so the number of bytes is rounded up to the nearest byte
            uint32_t padded_width = (image_data->width + 7) / 8;
            // Packing the pixel data into bits
            for (uint32_t row = 0; row < rows_in_strip; row++)
              {
                for (uint32_t col = 0; col < image_data->width; col++)
                {
                  uint8_t shift = 7 - col % 8;
                  strip_buf[row * padded_width + col / 8]
                    |= pixel_fvec[col] << shift;
                }
                pixel_fvec += image_data->width;
              }
            strip_size = padded_width * rows_in_strip;
            if (TIFFWriteEncodedStrip (tif, strip, strip_buf, strip_size) == -1)
              error ("Failed to write strip data to image");
          }
        else
          error ("Unsupported bit depth");
      }
  }

  template <typename T>
  void
  write_tiled_image (TIFF *tif, T pixel_data, tiff_image_data *image_data)
  {
    // ASSUMES pixel data dimensions are already validated

    uint32_t tile_width, tile_height;
    if (! TIFFGetField (tif, TIFFTAG_TILEWIDTH, &tile_width))
      error ("Failed to get the tile width");
    if (! TIFFGetField (tif, TIFFTAG_TILELENGTH, &tile_height))
      error ("Failed to get the tile length");

    if (tile_height == 0 || tile_height % 16 != 0
        || tile_width == 0 || tile_width % 16 != 0)
      error ("Tile dimesion tags are invalid");
    
    uint32_t tiles_across = (image_data->width + tile_width - 1)
                            / tile_width;
    uint32_t tiles_down = (image_data->height + tile_height - 1)
                          / tile_height;
    
    // Resize the image data to add tile padding
    dim_vector padded_dims (tiles_down * tile_height,
                            tiles_across * tile_width,
                            image_data->samples_per_pixel);
    pixel_data.resize (padded_dims);

    // Reshape the image to separate tiles into their own dimension
    // so we can permute them to the right order
    dim_vector tiled_dims (tile_height, tiles_down, tile_width, tiles_across,
                           image_data->samples_per_pixel);
    pixel_data = pixel_data.reshape (tiled_dims);

    // Permute the dimesnions to get the memory alignment to match LibTIFF
    Array<octave_idx_type> perm (dim_vector (5, 1));
    if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      {
        // For separate planes, the data coming from libtiff will have all
        // tiles of the first sample then all tiles of the second sample
        // and so on. Tiles of each sample will be ordered from left to right
        // and from top to bottom. And data inside each tile is organised as
        // rows and each row contains columns.
        // So the order for LibTIFF is:
        //   samples x tiles_down x tiles_across x tile_height x tile_width
        // But memory orientation of Octave arrays is reversed so we set it to
        //   tile_width x tile_height x tiles_across x tiles_down x samples
        perm(0) = 2;
        perm(1) = 0;
        perm(2) = 3;
        perm(3) = 1;
        perm(4) = 4;
      }
    else
      {
        // For chunky planes, the data coming from libtiff will be ordered
        // from left to right and from top to bottom. And data inside each
        // tile is organised as rows and each row contains columns and each
        // column contains samples.
        // So the order for LibTIFF is:
        //   tiles_down x tiles_across x tile_height x tile_width x samples
        // But memory orientation of Octave arrays is reversed so we set it to
        //   samples x tile_width x tile_height x tiles_across x tiles_down
        perm(0) = 4;
        perm(1) = 2;
        perm(2) = 0;
        perm(3) = 3;
        perm(4) = 1;
      }
    pixel_data = pixel_data.permute (perm);

    uint8_t *pixel_fvec
      = reinterpret_cast<uint8_t *> (pixel_data.fortran_vec ());
    uint32_t tile_count = TIFFNumberOfTiles (tif);
    tsize_t tile_size = TIFFTileSize (tif);

    for (uint32_t tile = 0; tile <tile_count; tile++)
      {
        if (image_data->bits_per_sample == 8
            || image_data->bits_per_sample == 16
            || image_data->bits_per_sample == 32
            || image_data->bits_per_sample == 64)
          {
            if (! TIFFWriteEncodedTile (tif, tile, pixel_fvec, tile_size))
              error ("Failed to write tile data");
            pixel_fvec += tile_size;
          }
        else if (image_data->bits_per_sample == 1)
          {
            if (image_data->samples_per_pixel != 1)
              error ("Bi-Level images must have one channel only");
            
            // Create a buffer to hold the packed tile data
            // Unique pointers are faster than vectors for
            // constant size buffers
            std::unique_ptr<uint8_t []> tile_ptr
              = std::make_unique<uint8_t []> (tile_size);
            uint8_t *tile_buf = tile_ptr.get ();
            // Packing the pixel data into bits
            for (uint32_t row = 0; row < tile_height; row++)
              {
                for (uint32_t col = 0; col < tile_width; col++)
                {
                  uint8_t shift = 7 - col % 8;
                  tile_buf[row * tile_width / 8 + col / 8]
                    |= pixel_fvec[col] << shift;
                }
                pixel_fvec += tile_width;
              }
            if (TIFFWriteEncodedTile (tif, tile, tile_buf,
                                      TIFFTileSize (tif)) == -1)
              error ("Failed to write tile data to image");
          }
        else
          error ("Unsupported bit depth");
      }

  }

  template <typename T>
  void
  write_image (TIFF *tif, T pixel_data, tiff_image_data *image_data)
  {
    // 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);

  }


#endif

  DEFUN (__open_tiff__, args, ,
             "Open a Tiff file and return its handle")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin == 0 || nargin > 2)
      {
        // TODO(maged): return invalid object instead??
        error ("No filename supplied\n");
      }

    std::string filename = args(0).string_value ();
    std::string mode = "r";

    if (nargin == 2)
      mode = args(1).string_value ();

    const std::vector<std::string> supported_modes {"r", "w", "w8", "a"};
      
    if (std::find (supported_modes.cbegin (), supported_modes.cend (), mode)
          == supported_modes.cend ())
      {
        if (mode == "r+")
          error ("Openning files in r+ mode is not yet supported");
        else
          error ("Invalid mode for openning Tiff file: %s", mode.c_str ());
      }
    
    TIFF *tif = TIFFOpen (filename.c_str (), mode.c_str ());
    
    if (! tif)
      error ("Failed to open Tiff file\n");

    // TODO(maged): use inheritance of octave_base_value instead
    octave_value tiff_ov = octave_value ((uint64_t)tif);
    return octave_value_list (tiff_ov);
#else
    octave_unused_parameter (args);
    err_disabled_feature ("Tiff", "Tiff");
#endif
  }


  DEFUN (__close_tiff__, args, ,
             "Close a tiff file")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin == 0)
      error ("No handle provided\n");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());
    TIFFClose (tif);

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


  DEFUN (__tiff_get_tag__, args, ,
             "Get the value of a tag from a tiff image")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    if (nargin == 0)
      error ("No handle provided\n");
    
    if (nargin < 2)
      error ("No tag name provided\n");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());

    uint32_t tag_id;
    const TIFFField *fip;
    
    if (args(1).is_string ())
      {
        std::string tagName = args(1).string_value ();
        fip = TIFFFieldWithName (tif, tagName.c_str ());
        if (! fip)
          error ("Tiff tag not found");
        
        tag_id = TIFFFieldTag (fip);
      }
    else
      {
        tag_id = args(1).int_value ();
        fip = TIFFFieldWithTag (tif, tag_id);
        
        if (! fip)
          error ("Tiff tag not found");
      }

    return octave_value_list (get_field_data (tif, fip));
#else
    err_disabled_feature ("getTag", "Tiff");
#endif
  }


  DEFUN (__tiff_set_tag__, args, ,
             "Set the value of a tag in a tiff image")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();
    
    if (nargin < 2)
      error ("Too few arguments provided\n");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());

    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 =  TIFFFieldWithName (tif, key.c_str ());
            if (! fip)
              error ("Tag %s not found", key.c_str ());
            set_field_data (tif, fip, tags.contents (key));
          }
      }
    else
      {
        if (nargin < 3)
          error ("Too few arguments provided");

        const TIFFField *fip;
        // TODO(maged): matlab actually checks for its own strings not LibTIFF's
        // Make it work in both ways
        if (args(1).is_string ())
          {
            std::string tagName = args(1).string_value ();
            fip = TIFFFieldWithName (tif, tagName.c_str ());
            if (! fip)
              error ("Tiff tag not found");
          }
        else
          {
            uint32_t tag_id = args(1).int_value ();
            fip = TIFFFieldWithTag (tif, tag_id);
            if (! fip)
              error ("Tiff tag not found");
          }
        
        set_field_data (tif, fip, args(2));
      }

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

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

    if (nargin == 0)
      error ("No handle provided\n");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());

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

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

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

    if (sample_format == 3)
      {
        if (image_data.bits_per_sample != 32 && image_data.bits_per_sample != 64)
          error ("Floating point images are only supported for bit depths of 32 and 64");
      }
    else if (sample_format != 1 && sample_format != 4)
      error ("Unsupported sample format");
    
    octave_value_list retval;
    switch (image_data.bits_per_sample)
      {
      case 1:
        retval(0) = read_image<boolNDArray> (tif, &image_data);
        break;
      case 4:
      case 8:
        retval(0) = read_image<uint8NDArray> (tif, &image_data);
        break;
      case 16:
        retval(0) = read_image<uint16NDArray> (tif, &image_data);
        break;
      case 32:
        if (sample_format == 3)
          retval(0) = read_image<FloatNDArray> (tif, &image_data);
        else
          retval(0) = read_image<uint32NDArray> (tif, &image_data);
        break;
      case 64:
        if (sample_format == 3)
          retval(0) = read_image<NDArray> (tif, &image_data);
        else
          retval(0) = read_image<uint64NDArray> (tif, &image_data);
        break;
      default:
        error ("Unsupported bit depth");
      }
    
    return retval;
#else
    err_disabled_feature ("read", "Tiff");
#endif
  }

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

    if (nargin != 2)
      error ("Wrong number of arguments");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());

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

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

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

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

    if (nargin != 2)
      error ("Wrong number of arguments");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());

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

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

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

  DEFUN (__tiff_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");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());

    tiff_image_data image_data (tif);

    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 ());
    // Matlab (R2022a) uses a top-left orientation ignoring the tag if present
    if (! TIFFReadRGBAImageOriented (tif, image_data.width, image_data.height,
                                     img_ptr, 1))
      error ("Failed to read image");
    
    Array<octave_idx_type> perm (dim_vector (3, 1));
    perm(0) = 2;
    perm(1) = 1;
    perm(2) = 0;
    img = img.permute (perm);

    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 (img.index (idx));
    idx(2) = idx_vector (3);
    uint8NDArray alpha = uint8NDArray (img.index (idx));
    
    octave_value_list retval (2);
    retval(0) = octave_value (rgb);
    retval(1) = octave_value (alpha);
    return retval;
#else
    err_disabled_feature ("readEncodedTile", "Tiff");
#endif
  }

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

    if (nargin < 2)
      error ("Wrong number of arguments\n");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());

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

    if (sample_format == 3)
      {
        if (image_data.bits_per_sample != 32
            && image_data.bits_per_sample != 64)
          error ("Floating point images are only supported for bit depths of 32 and 64");
      }
    else if (sample_format != 1 && sample_format != 4)
      error ("Unsupported sample format");
    
    switch (image_data.bits_per_sample)
      {
      case 1:
        // We need to check for both scalar and matrix types to handle single
        // pixel image
        if (args (1).is_bool_scalar () || args (1).is_bool_matrix ())
          write_image<boolNDArray> (tif, args (1).bool_array_value (),
                                    &image_data);
        else
          error ("Expected logical matrix for BiLevel image");
        break;
      case 8:
        if (args (1).is_uint8_type ())
          write_image<uint8NDArray> (tif, args (1).uint8_array_value (),
                                     &image_data);
        else
          error ("Only uint8 data is allowed for uint images with bit depth of 8");
        break;
      case 16:
        if (args (1).is_uint16_type ())
          write_image<uint16NDArray> (tif, args (1).uint16_array_value (),
                                      &image_data);
        else
          error ("Only uint16 data is allowed for uint images with bit depth of 16");
        break;
      case 32:
        if (sample_format == 3)
          if (args (1).is_single_type () || args (1).is_double_type ())
            write_image<FloatNDArray> (tif, args (1).float_array_value (),
                                       &image_data);
          else
            error ("Only single and double data are allowed for floating-point images");
        else
          if (args (1).is_uint32_type ())
            write_image<uint32NDArray> (tif, args (1).uint32_array_value (),
                                        &image_data);
          else
            error ("Only uint32 data is allowed for uint images with bit depth of 32");
        break;
      case 64:
        if (sample_format == 3)
          if (args (1).is_single_type () || args (1).is_double_type ())
            write_image<NDArray> (tif, args (1).array_value (), &image_data);
          else
            error ("Only single and double data are allowed for floating-point images");
        else  
          if (args (1).is_uint64_type ())
            write_image<uint64NDArray> (tif, args (1).uint64_array_value (),
                                        &image_data);
          else
            error ("Only uint64 data is allowed for uint images with bit depth of 64");
        break;
      default:
        error ("Unsupported bit depth");
      }
    
    return octave_value_list ();
#else
    err_disabled_feature ("write", "Tiff");
#endif
  }

  DEFUN (__tiff_write_encoded_strip__, args, ,
             "Write an encoded strip to the image")
  {
#if defined (HAVE_TIFF)
    int nargin = args.length ();

    // TODO(maged): add support for YCbCr data
    if (nargin < 3)
      error ("Too few arguments provided\n");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());

    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");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());

    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");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());
    bool is_tiled = static_cast<bool> (TIFFIsTiled (tif));
    return octave_value_list (octave_value (is_tiled));
#else
    err_disabled_feature ("isTiled", "Tiff");
#endif
  }

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

    if (nargin == 0)
      error ("No handle provided\n");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());
    
    if (TIFFIsTiled (tif))
      error ("The image is tiled not stripped");
    
    double strip_count = static_cast<double> (TIFFNumberOfStrips (tif));
    return octave_value_list (octave_value (strip_count));
#else
    err_disabled_feature ("numberOfStrips", "Tiff");
#endif
  }

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

    if (nargin == 0)
      error ("No handle provided\n");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());
    
    if (! TIFFIsTiled (tif))
      error ("The image is stripped not tiled");
    
    double tile_count = static_cast<double> (TIFFNumberOfTiles (tif));
    return octave_value_list (octave_value (tile_count));
#else
    err_disabled_feature ("numberOfTiles", "Tiff");
#endif
  }

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

    if (nargin < 2 || nargin > 3)
      error ("Wrong number of arguments\n");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());
    
    if (TIFFIsTiled (tif))
      error ("The image is tiled not stripped");
    
    tiff_image_data image_data (tif);

    uint32_t row = args(1).uint32_scalar_value ();
    if (row > image_data.height)
      row = image_data.height;

    // Convert from 1-based to zero-based indexing but avoid underflow
    if (row > 0)
      row--;
    
    uint16_t plane;
    if (nargin > 2)
      {
        if (image_data.planar_configuration == PLANARCONFIG_CONTIG)
          error ("Can't use plane argument for images with chunky PlanarConfiguration");
        plane = args(2).uint16_scalar_value ();
        if (plane > image_data.samples_per_pixel)
          plane = image_data.samples_per_pixel;
        if (plane > 0)
          plane--;
      }
    else
      {
        if (image_data.planar_configuration == PLANARCONFIG_SEPARATE)
          error ("The plane argument is required for images with separate PlanarConfiguration");
        plane = 0;
      }

    double strip_number = TIFFComputeStrip (tif, row, plane) + 1;
    if (strip_number > TIFFNumberOfStrips (tif))
      strip_number = TIFFNumberOfStrips (tif);
    return octave_value_list (octave_value (strip_number));
#else
    err_disabled_feature ("computeStrip", "Tiff");
#endif
  }

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

    if (nargin < 2 || nargin > 3)
      error ("Wrong number of arguments\n");
    
    TIFF *tif = (TIFF *)(args(0).uint64_value ());
    
    if (! TIFFIsTiled (tif))
      error ("The image is stripped not tiled");
    
    uint32NDArray coords = args(1).uint32_array_value ();
    if (coords.dim2() < 2)
      error ("Coordinates must be in the shape [row, col]");
    uint32_t row = coords(0, 0);
    uint32_t col = coords(0, 1);

    tiff_image_data image_data (tif);

    if (col > image_data.width)
      col = image_data.width;
    if (row > image_data.height)
      row = image_data.height;

    // Convert from 1-based to zero-based indexing but avoid underflow
    if (row > 0)
      row--;
    if (col > 0)
      col--;
    
    uint16_t plane;
    if (nargin > 2)
      {
        if (image_data.planar_configuration == PLANARCONFIG_CONTIG)
          error ("Can't use plane argument for images with chunky PlanarConfiguration");
        plane = args(2).uint16_scalar_value ();
        if (plane > image_data.samples_per_pixel)
          plane = image_data.samples_per_pixel;
        if (plane > 0)
          plane--;
      }
    else
      {
        if (image_data.planar_configuration == PLANARCONFIG_SEPARATE)
          error ("The plane argument is required for images with separate PlanarConfiguration");
        plane = 0;
      }

    double tile_number = TIFFComputeTile (tif, col, row, 0, plane) + 1;
    if (tile_number > TIFFNumberOfTiles (tif))
      tile_number = TIFFNumberOfTiles (tif);
    return octave_value_list (octave_value (tile_number));
#else
    err_disabled_feature ("computeTile", "Tiff");
#endif
  }

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

  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
  }

}