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

#include <string>
#include <iostream>

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

#include "errwarn.h"

#if defined (HAVE_TIFF)
#  include <tiffio.h>
#endif

// TODO(maged): Fix warnings

namespace octve
{
#if defined (HAVE_TIFF)

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

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

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

      if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
                                   &bits_per_sample))
        error ("Failed to read the BitsPerSample tag");
      
      if (! TIFFGetFieldDefaulted (tif, TIFFTAG_PLANARCONFIG,
                                   &planar_configuration))
        error ("Failed to read the PlanarConfiguration tag");
      
      is_tiled = TIFFIsTiled(tif);
    }
  };

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

  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
        // TODO(maged): Are incorrect sized strips checked internally?
        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
        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;
                img_fvec[pixel]
                  = (((uint8_t *)img_fvec)[pixel / 8] >> bit_number) & 0x01;
              }
            break;
          }
        else if (image_data->bits_per_sample == 4)
          {
            if (image_data->samples_per_pixel != 1)
              error ("4-bit images are only supported for grayscale");
            
            // 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;
                img_fvec[pixel]
                  = (((uint8_t *)img_fvec)[pixel / 2] >> shift) & 0x0F;
              }
            break;
          }
        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 = (P*)((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;    
    uint32_t written_size = 0;
    for (uint32_t tile = 0; tile < tile_count; tile++)
      {
        tile_size = TIFFReadEncodedTile(tif, tile, img_fvec, -1);
        
        if (tile_size == -1)
          error ("Failed to read tile data");

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

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

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

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

        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 = (double)(*((uint8_t *)data));
          break;
        }
      case TIFF_SHORT:
        {
          retval = (double)(*((uint16_t *)data));
          break;
        }
      case TIFF_LONG:
        {
          retval = (double)(*((uint32_t *)data));
          break;
        }
      case TIFF_LONG8:
        {
          retval = (double)(*((uint64_t *)data));
          break;
        }
      case TIFF_RATIONAL:
        {
          error ("TIFF_RATIONAL should have at least 2 elements but got only 1");
          break;
        }
      case TIFF_SBYTE:
        {
          retval = (double)(*((int8_t *)data));
          break;
        }
      case TIFF_SSHORT:
        {
          retval = (double)(*((int16_t *)data));
          break;
        }
      case TIFF_SLONG:
        {
          retval = (double)(*((int32_t *)data));
          break;
        }
      case TIFF_SLONG8:
        {
          retval = (double)(*((int64_t *)data));
          break;
        }
      case TIFF_FLOAT:
        {
          retval = *((float *)data);
          break;
        }
      case TIFF_DOUBLE:
        {
          retval = *((double *)data);
          break;
        }
      case TIFF_SRATIONAL:
        {
          error ("TIFF_SRATIONAL should have at least 2 elements but got only 1");
          break;
        }
      case TIFF_IFD:
      case TIFF_IFD8:
        error ("Unimplemented IFFD data type");
        break;
      default:
        error ("Unsupported tag data type");
      }

    return octave_value (retval);
  }

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

        switch (tag_datatype)
          {
          case TIFF_BYTE:
          case TIFF_UNDEFINED:
            {
              uint8NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = ((uint8_t *)data)[i];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_ASCII:
            {
              retval = octave_value (*(char **)data);
              break;
            }
          case TIFF_SHORT:
            {
              uint16NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = ((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) = ((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) = ((uint64_t *)data)[i];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_RATIONAL:
            {
              NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i+=2)
                {
                  arr(i / 2) = (float)((uint32_t *)data)[i] 
                                / (float)((uint32_t *)data)[i+1];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_SBYTE:
            {
              int8NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i++)
                {
                  arr(i) = ((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) = ((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) = ((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) = ((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) = ((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) = ((double *)data)[i];
              }
              retval = octave_value (arr);
              break;
            }
          case TIFF_SRATIONAL:
            {
              NDArray arr (arr_dims);
              for (uint32_t i = 0; i < count; i+=2)
                {
                  arr(i / 2) = (float)((int32_t *)data)[i] 
                                / (float)((int32_t *)data)[i+1];
                }
              retval = octave_value (arr);
              break;
            }
          case TIFF_IFD:
          case TIFF_IFD8:
            // TODO(maged): implement IFD datatype?
            error ("Unimplemented IFFD 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));
    // TODO(maged): use shared pointer instead of malloc
    void *data = _TIFFmalloc (type_size);
    // TODO(maged): check if this should be GetFieldDefaulted instead
    validate_tiff_get_field (TIFFGetField (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");
                
          if (bits_per_sample > 24)
            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));
          
          // TODO(maged): use Array<T>::cat
          Matrix mat_out (count, 3);
          
          Matrix red_array (interpret_tag_data (red,
                                                count,
                                                TIFFFieldDataType (fip))
                                                .uint16_array_value ());
          Matrix green_array (interpret_tag_data (green,
                                                  count,
                                                  TIFFFieldDataType (fip))
                                                  .uint16_array_value ());
          Matrix blue_array (interpret_tag_data (blue,
                                                 count,
                                                 TIFFFieldDataType (fip))
                                                 .uint16_array_value ());
          
          double *out_ptr = mat_out.fortran_vec ();
          memcpy (out_ptr, red_array.fortran_vec (), sizeof(double)*count);
          out_ptr += count;
          memcpy (out_ptr, green_array.fortran_vec (), sizeof(double)*count);
          out_ptr += count;
          memcpy (out_ptr, blue_array.fortran_vec (), sizeof(double)*count);

          // Normalize the range to be between 0 and 1
          mat_out /= UINT16_MAX;

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

          uint16_t bits_per_sample;
          if (! TIFFGetFieldDefaulted (tif, TIFFTAG_BITSPERSAMPLE,
                                       &bits_per_sample))
            error ("Failed to obtain the number of samples per pixel");

          uint32_t count = 1 << bits_per_sample;
          uint16_t *ch1, *ch2, *ch3;
          if (samples_per_pixel == 1)
            {
              validate_tiff_get_field (TIFFGetField (tif, TIFFTAG_COLORMAP, &ch1));
              tag_data_ov = interpret_tag_data (ch1, count,
                                                TIFFFieldDataType (fip));
            }
          else
            {
              validate_tiff_get_field (TIFFGetField (tif, TIFFTAG_COLORMAP,
                                                     &ch1, &ch2, &ch3));
              
              uint16NDArray mat_out (dim_vector (count, 3));

              uint16NDArray ch1_array
                = interpret_tag_data (ch1,
                                      count,
                                      TIFFFieldDataType (fip)).uint16_array_value ();
              uint16NDArray ch2_array
                = interpret_tag_data (ch2,
                                      count,
                                      TIFFFieldDataType (fip)).uint16_array_value ();
              uint16NDArray ch3_array
                = interpret_tag_data (ch3,
                                      count,
                                      TIFFFieldDataType (fip)).uint16_array_value ();
              
              octave_uint16 *out_ptr = mat_out.fortran_vec ();
              memcpy (out_ptr, ch1_array.fortran_vec (), sizeof(uint16_t) * count);
              out_ptr += count;
              memcpy (out_ptr, ch2_array.fortran_vec (), sizeof(uint16_t) * count);
              out_ptr += count;
              memcpy (out_ptr, ch3_array.fortran_vec (), sizeof(uint16_t) * count);

              tag_data_ov = octave_value (mat_out);
            }
          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));
          
          Matrix mat_out (1, 2);
          mat_out(0)
            = interpret_tag_data (&tag_part1, 1,
                                  TIFFFieldDataType (fip)).double_value ();
          mat_out(1)
            = interpret_tag_data (&tag_part2, 1,
                                  TIFFFieldDataType (fip)).double_value ();
          
          tag_data_ov = octave_value (mat_out);
          break;
        }
      case TIFFTAG_SUBIFD:
        {
          uint16_t count;
          uint64_t *offsets;
          validate_tiff_get_field (TIFFGetField (tif, tag_id, &count, &offsets));
          tag_data_ov = interpret_tag_data (offsets, count,
                                            TIFFFieldDataType (fip));
          break;
        }
      case TIFFTAG_EXTRASAMPLES:
        {
          uint16_t count;
          uint16_t *types;
          validate_tiff_get_field (TIFFGetField (tif, tag_id, &count, &types));
          tag_data_ov = interpret_tag_data (types, count,
                                            TIFFFieldDataType (fip));
          break;
        }
      // TODO(maged): These tags are more complex to implement
      //  will be implemented and tested later.
      case TIFFTAG_XMLPACKET:
      case TIFFTAG_RICHTIFFIPTC:
      case TIFFTAG_PHOTOSHOP:
      case TIFFTAG_ICCPROFILE:
        {
          error ("Complex Tags not implemented");
          break;
        }
      // These tags are not mentioned in the LibTIFF documentation
      // but are handled correctly by the library
      case TIFFTAG_ZIPQUALITY:
      case TIFFTAG_SGILOGDATAFMT:
      case TIFFTAG_GRAYRESPONSEUNIT:
        {
          tag_data_ov = get_scalar_field_data (tif, fip);
          break;
        }
      default:
        tag_data_ov = get_scalar_field_data (tif, fip);
    }
    
    return tag_data_ov;
  }

  void
  set_field_data (TIFF *tif, const TIFFField *fip, octave_value tag_ov)
  {
    // TODO(maged): complete the implementation of this function
    uint32_t tag_id = TIFFFieldTag (fip);
    uint32_t tag_data = tag_ov.double_value ();

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

  template <typename T>
  void
  write_strip (TIFF *tif, uint32_t strip_no, octave_value strip_data_ov,
               tiff_image_data *image_data)
  {
    T strip_data = T (strip_data_ov.array_value ());
    
    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;
    
    uint32_t strip_count = TIFFNumberOfStrips (tif);
    uint32_t expected_numel;

    // Calculate the expected number of elements in the strip data array
    // All strips have equal number of rows excpet strips at the bottom
    // of the image can have less number of rows
    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)
          rows_in_strip = image_data->height - rows_in_strip * (strip_no - 1);
        expected_numel = rows_in_strip * image_data->width
                        * image_data->samples_per_pixel;
      }
    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;
             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 - 1)
                                                    % (strips_per_channel));
        expected_numel = rows_in_strip * image_data->width;
      }
    else
      error ("Planar configuration not supported");

    if (strip_data.numel () != expected_numel)
      error ("Size of strip data is different from the expected size of the strip");

    
    // Put the strip in a consistent shape for the subsequent permutation
    if (image_data->planar_configuration == PLANARCONFIG_CONTIG)
      {
        strip_data
          = strip_data.reshape (dim_vector (rows_in_strip, image_data->width,
                                            image_data->samples_per_pixel));
      }
    else if (image_data->planar_configuration == PLANARCONFIG_SEPARATE)
      {
        strip_data = strip_data.reshape (dim_vector (rows_in_strip, 
                                                     image_data->width, 1));
      }

    //TODO(maged): add suppot for 1-bit and 4-bit images

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

    // LibTIFF uses zero-based indexing as opposed to Octave's 1-based
    strip_no--;
    // Can't rely in LibTIFF's TIFFStripSize because boundary strips
    // can be smaller in size
    tsize_t strip_size = expected_numel * image_data->bits_per_sample / 8;
    if (TIFFWriteEncodedStrip (tif, strip_no,
                               strip_data.fortran_vec (), strip_size) == -1)
      error ("Failed to write strip data to image");
  }

#endif

  DEFUN_DLD (__open_tiff__, args, nargout,
             "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";

    // TODO(maged): check valid mode
    if (nargin == 2)
      mode = args (1).string_value ();

    std::vector<std::string> supported_modes {"r", "w", "w8", "a"};
      
    if (std::find(supported_modes.begin(), supported_modes.end(), mode)
        == supported_modes.end())
      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 ());
    
    // TODO(maged): Look into unwind action
    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
    err_disabled_feature ("Tiff", "Tiff");
#endif
  }


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

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

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


  DEFUN_DLD (__tiff_get_tag__, args, nargout,
             "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).type_name () == "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);
        // TODO(maged): Handle other types of errors (e.g. unsupported tags)
        if (! fip)
          error ("Tiff tag not found");
      }

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


  DEFUN_DLD (__tiff_set_tag__, args, nargout,
             "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 ());
    
    // TODO(maged): does matlab allow calling this function for images
    // opened for reading?
    if (args (1).type_name () == "struct")
      error ("setTag with struct is not yet supported");
    else
      {
        const TIFFField *fip;
        if (args (1).type_name () == "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");
          }
        
        if (nargin < 3)
          error ("Too few arguments provided");
        
        set_field_data (tif, fip, args (2));
      }

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

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

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

    // TODO(maged): nargout and ycbcr

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

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

    if (nargin < 3)
      error ("Too few arguments provided\n");
    
    TIFF *tif = (TIFF *)(args (0).uint64_value ());

    // Obtain all necessary tags
    tiff_image_data image_data (tif);

    if (image_data.is_tiled)
      error ("Can't rite strips to a tiled image");

    uint32_t strip_no = args (1).uint_value ();
    if (strip_no < 1 || strip_no > TIFFNumberOfStrips (tif))
      error ("Strip number out of range");
    
    // 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");

    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");
    
    std::string type_name = args(2).type_name ();
    switch (image_data.bits_per_sample)
      {
      case 1:
        if (type_name == "bool matrix")
          write_strip<boolNDArray> (tif, strip_no, args(2), &image_data);
        else
          error ("Expected logical matrix for BiLevel image");
        break;
      case 4:
      case 8:
        if (type_name == "uint8 matrix" || type_name == "int8 matrix")
          write_strip<uint8NDArray> (tif, strip_no, args(2), &image_data);
        else
          error ("Only uint8 and int8 data are allowed for images with bit depth of 8");
        break;
      case 16:
        // Type conversion from signed to unsigned is handled in the function
        // TODO(maged): what is the behavior if the input matrix has
        // negative numbers?
        if (type_name == "uint16 matrix" || type_name == "int16 matrix")
          write_strip<uint16NDArray> (tif, strip_no, args(2), &image_data);
        else
          error ("Only uint16 and int16 data are allowed for images with bit depth of 16");
        break;
      case 32:
        if (sample_format == 3)  
          if (type_name == "float matrix" || type_name == "matrix")
            write_strip<FloatNDArray> (tif, strip_no, args(2), &image_data);
          else
            error ("Only single and double data are allowed for floating-point images");
        else
          if (type_name == "uint32 matrix" || type_name == "int32 matrix")
            write_strip<uint32NDArray> (tif, strip_no, args(2), &image_data);
          else
            error ("Only uint32 and int32 data are allowed for images with bit depth of 32");
        break;
      case 64:
        if (sample_format == 3)          
          if (type_name == "float matrix" || type_name == "matrix")
            write_strip<NDArray> (tif, strip_no, args(2), &image_data);
          else
            error ("Only single and double data are allowed for floating-point images");
        else  
          if (type_name == "uint64 matrix" || type_name == "int64 matrix")
            write_strip<uint64NDArray> (tif, strip_no, args(2), &image_data);
          else
            error ("Only uint64 and int64 data are allowed for images with bit depth of 64");
        break;
      default:
        error ("Unsupported bit depth");
      }
    
    return octave_value_list ();
#else
    err_disabled_feature ("writeEncodedStrip", "Tiff");
#endif
  }
}