/* 
1 - Obtain bigtiff lib package from http://bigtiff.org/#API_CHANGES
    Source ZIP for libtiff_4.1 only.

2 - unzip libtiff-4.1_only.zip
    cd tiff-4.1/libtiff
    make
    cp libtiff.a to minc quarantine Linux-x86_64/lib/

3 - gcc -DMINC2 -I/home/claude/QUARANTINES/Linux-x86_64/include \
    -I/home/claude/QUARANTINES/TIFF/tiff-4.1/libtiff tiff2mnc.c \
    -L/home/claude/QUARANTINES/Linux-x86_64/lib -o tiff2mnc -ltiff \
    -ljpeg -lminc2 -lnetcdf -lhdf5 -lm -lz
*/


#include <stdio.h>
#include <stdlib.h>
#include <limits.h>

#include <tiffio.h>

#include <minc.h>

#define NDIMS 3
#define COMPRESS_LEVEL 4
#define TILEAVG
// #define INVERT

extern char *time_stamp(int, char **);

int
main(int argc, char **argv)
{
  TIFF *r_fp;
  int w_fd;
  uint32 x, y, yy, imageWidth, imageLength, tileWidth, tileLength, tilesize;
  uint32 outputImageWidth, outputImageLength, outputTileWidth, outputTileLength;
  uint32 vxmin, vxmax, vymin, vymax;
  uint16 bps;
  double max;
  double min;
  tstrip_t i;
  unsigned short **tiledata;
  unsigned short *mincdata;
  double vrange[2];
  int dimids[NDIMS];
  long count[NDIMS];
  long start[NDIMS];
  int j, ii, jj, flip;
  int var_id;
  char *history_str = time_stamp(argc, argv);
  double xstep, ystep;
  float xres, yres;
  uint16 resunit;

  if (argc < 3) {
    fprintf(stderr, "Not enough arguments\n");
  }

  r_fp = TIFFOpen(argv[1], "r");
  if (r_fp == NULL) {
    fprintf(stderr, "Can't open '%s' for reading\n", argv[1]);
    exit(-1);
  }

  if (TIFFGetField(r_fp, TIFFTAG_IMAGEWIDTH, &imageWidth)) {
    printf("(image width %d) ", imageWidth);
  }

  if (TIFFGetField(r_fp, TIFFTAG_IMAGELENGTH, &imageLength)) {
    printf("(image length %d) ", imageLength);
  }

  if (TIFFGetField(r_fp, TIFFTAG_TILEWIDTH, &tileWidth)) {
    printf("(tile width %d) ", tileWidth);
  }

  if (TIFFGetField(r_fp, TIFFTAG_TILELENGTH, &tileLength)) {
    printf("(tile length %d) ", tileLength);
  }

  if (TIFFGetField(r_fp, TIFFTAG_BITSPERSAMPLE, &bps)) {
    printf("(%d bps) ", bps);
  }

  if (TIFFGetField(r_fp, TIFFTAG_XRESOLUTION, &xres)) {
    // xres = 24000;
    printf("(xres %f) ", xres);
  }

  if (TIFFGetField(r_fp, TIFFTAG_YRESOLUTION, &yres)) {
    // yres = 24000;
    printf("(yres %f) ", yres);
  }

  if (TIFFGetField(r_fp, TIFFTAG_RESOLUTIONUNIT, &resunit)) {
    printf("(unit %d)", resunit);
  }

  printf("\n");

#if 0
  if (resunit == RESUNIT_INCH) {
    printf( "Units are in inches...\n" );
    xstep = 25.4 / xres;        /* convert to mm */
    ystep = 25.4 / yres;        /* convert to mm */
  } else if (resunit == RESUNIT_CENTIMETER) {
    printf( "Units are in cm...\n" );
    xstep = 10.0 / xres;
    ystep = 10.0 / yres;
  }
#else
  xstep = 0.001 * xres;  // units are read in microns, convert to mm
  ystep = 0.001 * yres;
#endif

  flip = ( imageWidth < imageLength );
  if( flip ) {
    outputImageWidth = imageLength;
    outputImageLength = imageWidth;
    outputTileWidth = tileLength;
    outputTileLength = tileWidth;
    double junk = xstep;
    xstep = ystep;
    ystep = -junk;
  } else {
    outputImageWidth = imageWidth;
    outputImageLength = imageLength;
    outputTileWidth = tileWidth;
    outputTileLength = tileLength;
  }
 
  // Open the minc file for writing

#ifdef MINC2
  struct mi2opts opts;
  opts.struct_version = MI2_OPTS_V1;
  opts.comp_type = MI2_COMP_ZLIB;
  opts.comp_param = COMPRESS_LEVEL;       // compress image
  opts.chunk_type = MI2_CHUNK_UNKNOWN;    // will use default chunking
  w_fd = micreatex(argv[2], NC_NOCLOBBER|MI2_CREATE_V2, &opts);
#else
  w_fd = micreate(argv[2], NC_NOCLOBBER);
#endif
  if (w_fd < 0) {
    fprintf(stderr, "Can't open '%s' for writing\n", argv[1]);
    exit(-1);
  }

#ifdef TILEAVG
  dimids[0] = ncdimdef(w_fd, MIzspace, 1);
  dimids[1] = ncdimdef(w_fd, MIyspace, 
                       (int)((outputImageLength+outputTileLength-1)/outputTileLength) );
  dimids[2] = ncdimdef(w_fd, MIxspace, 
                       (int)((outputImageWidth+outputTileWidth-1)/outputTileWidth) );
#else

#if 0
  // to create a subset image only
  vxmin = (int)( ( 0.0 * outputImageWidth ) / outputTileWidth ) * outputTileWidth;
  vxmax = (int)( ( 0.3571 * outputImageWidth ) / outputTileWidth ) * outputTileWidth - 1;
  vymin = (int)( ( 0.0 * outputImageLength ) / outputTileLength ) * outputTileLength;
  vymax = (int)( ( 1.0 * outputImageLength ) / outputTileLength ) * outputTileLength - 1;
#else
  // full image
  vxmin = 0;
  vxmax = outputImageWidth-1;
  vymin = 0;
  vymax = outputImageLength-1;
#endif

  dimids[0] = ncdimdef(w_fd, MIzspace, 1);
  dimids[1] = ncdimdef(w_fd, MIyspace, vymax - vymin + 1 );
  dimids[2] = ncdimdef(w_fd, MIxspace, vxmax - vxmin + 1 );
#endif

  if( bps == 16 ) {
    micreate_std_variable(w_fd, MIimage, NC_SHORT, 3, dimids);
  } else if( bps == 8 ) {
    micreate_std_variable(w_fd, MIimage, NC_BYTE, 3, dimids);
  }
  micreate_std_variable(w_fd, MIimagemin, NC_DOUBLE, 0, NULL);
  micreate_std_variable(w_fd, MIimagemax, NC_DOUBLE, 0, NULL);

  micreate_std_variable(w_fd, MIzspace, NC_INT, 0, NULL);
  micreate_std_variable(w_fd, MIyspace, NC_INT, 0, NULL);
  micreate_std_variable(w_fd, MIxspace, NC_INT, 0, NULL);

  miattputstr(w_fd, NC_GLOBAL, MIhistory, history_str);

  var_id = ncvarid(w_fd, MIimage);
  miattputstr(w_fd, var_id, MIsigntype, MI_UNSIGNED);
  miattputstr(w_fd, var_id, MIcomplete, MI_TRUE);

  if( bps == 16 ) {
    vrange[0] = 0;
    vrange[1] = USHRT_MAX;
  } else {
    vrange[0] = 0;
    vrange[1] = UCHAR_MAX;
  }

  miset_valid_range(w_fd, var_id, vrange);
  var_id = ncvarid(w_fd, MIxspace);
  miattputstr(w_fd, var_id, MIunits, "mm");
#ifdef TILEAVG
  miattputdbl(w_fd, var_id, MIstep, xstep*outputTileWidth);
#else
  miattputdbl(w_fd, var_id, MIstep, xstep);
#endif
  miattputdbl(w_fd, var_id, MIstart, -(xstep * outputImageWidth) / 2.0 );
  printf( "x from %f by %f for %d voxels\n", -(xstep * outputImageWidth) / 2.0,
          xstep, outputImageWidth );

  var_id = ncvarid(w_fd, MIyspace);
  miattputstr(w_fd, var_id, MIunits, "mm");
#ifdef TILEAVG
  miattputdbl(w_fd, var_id, MIstep, -ystep*outputTileLength);
#else
  miattputdbl(w_fd, var_id, MIstep, -ystep);
#endif
  miattputdbl(w_fd, var_id, MIstart, (ystep * outputImageLength) / 2.0 );
  printf( "y from %f by %f for %d voxels\n", (ystep * outputImageLength) / 2.0,
          -ystep, outputImageLength );

  var_id = ncvarid(w_fd, MIzspace);
  miattputstr(w_fd, var_id, MIunits, "mm");
  miattputdbl(w_fd, var_id, MIstep, 0.02);
  miattputdbl(w_fd, var_id, MIstart, 0.0);
  printf( "z from %f by %f for %d voxels\n", 0.0, 0.02, 1 );

  ncendef(w_fd);

  max = vrange[0];
  min = vrange[1];

  var_id = ncvarid(w_fd, MIimage);

  count[0] = 1;  // number of slices per tif
  count[1] = 1;  // number of rows per strip
#ifdef TILEAVG
  count[2] = (int)((outputImageWidth+outputTileWidth-1)/outputTileWidth);  // number of tiles per row
#else
  count[2] = vxmax - vxmin + 1;  // number of columns per row
#endif
  start[0] = 0;
  start[1] = 0;
  start[2] = 0;

  tilesize = TIFFTileSize(r_fp);
  printf("tilesize %d\n", tilesize);

  short numtiles = ( outputImageWidth + outputTileWidth - 1 ) / outputTileWidth;
  printf( "numtiles = %d\n", numtiles );
  tiledata = (unsigned short**)malloc( numtiles * sizeof( unsigned short * ) );
  for( j = 0; j < numtiles; j++ ) {
    tiledata[j] = (unsigned short *)malloc( tilesize );
    if( tiledata[j] == NULL ) {
      printf( "Cannot allocate memory for tile %d\n", j );
      exit(1);
    }
  }

#ifdef TILEAVG
  mincdata = (unsigned short *)malloc( count[2] * sizeof( unsigned short ) );
#else
  mincdata = (unsigned short * )malloc( outputImageWidth * sizeof( unsigned short ) );
#endif
  unsigned char * bytedata = (unsigned char *)mincdata;

  for( y = 0; y < outputImageLength; y += outputTileLength ) {
    printf( "Processing tiles at y = %d...\n", y );
    for( x = 0; x < outputImageWidth; x += outputTileWidth ) {
      // z=0 is the depth of the image
      if( flip ) {
        TIFFReadTile( r_fp, tiledata[x/outputTileWidth], y, x, 0, 0 );
      } else {
        TIFFReadTile( r_fp, tiledata[x/outputTileWidth], x, y, 0, 0 );
      }
    }
#ifdef TILEAVG
    // average each tile into a voxel. This assumes that the tile
    // is complete and ignores incomplete tiles near the edge of the
    // image. This is not too important since this image is only used
    // for quality control.
    for( x = 0; x < outputImageWidth; x += outputTileWidth ) {
      double avg = 0.0;
      j = x/outputTileWidth;

      // Note: internal ordering within a tile is NOT flipped
      if( bps == 8 ) {
        unsigned char * ptr = (unsigned char *)tiledata[j];
        for( ii = 0; ii < tilesize ; ii++ ) {
          avg += ptr[ii];
        }
      } else if( bps == 16 ) {
        for( ii = 0; ii < tilesize; ii++ ) {
          avg += tiledata[j][ii];
        }
      }
      avg /= (double)tilesize;
      if (avg > max) {
        max = avg;
      }
      if (avg < min) {
        min = avg;
      }
      if( bps == 8 ) {
#ifdef INVERT
        bytedata[j] = (unsigned char)( vrange[1] - avg );
#else
        bytedata[j] = (unsigned char)avg;
#endif
      } else {
#ifdef INVERT
        bytedata[j] = (unsigned short)( vrange[1] - avg );
#else
        mincdata[j] = (unsigned short)avg;
#endif
      }
    }
    if( bps == 8 ) {
      mivarput( w_fd, var_id, start, count, NC_BYTE, MI_UNSIGNED, bytedata );
    } else if( bps == 16 ) {
      mivarput( w_fd, var_id, start, count, NC_SHORT, MI_UNSIGNED, mincdata );
    }
    start[1] += count[1];
#else
    // read line-by-line from the tiles
    for( yy = 0; yy < outputTileLength && y + yy < outputImageLength; yy++ ) {
      if( y+yy >= vymin && y+yy <= vymax ) {
        for( x = 0; x < outputImageWidth; x++ ) {

          unsigned short val = 0;
          if (bps == 8) {
            unsigned char * ptr = (unsigned char *)tiledata[x/outputTileWidth];
            /* scale up from byte to short */

            // Note: internal ordering within a tile is NOT flipped
            if( flip ) {
              val = (unsigned short)ptr[(x%tileWidth)*tileWidth+yy];
            } else {
              val = (unsigned short)ptr[yy*tileWidth+(x%tileWidth)];
            }
          } else if( bps == 16 ) {
            if( flip ) {
              val = (unsigned short)tiledata[x/outputTileWidth][(x%tileWidth)*tileWidth+yy];
            } else {
              val = (unsigned short)tiledata[x/outputTileWidth][yy*tileWidth+(x%tileWidth)];
            }
          }
          if (val > max) {
            max = val;
          }
          if (val < min) {
            min = val;
          }
          if( bps == 8 ) {
#ifdef INVERT
            bytedata[x] = (unsigned char)( vrange[1] - val );
#else
            bytedata[x] = (unsigned char)val;
#endif
          } else {
#ifdef INVERT
            bytedata[x] = (unsigned short)( vrange[1] - val );
#else
            mincdata[x] = (unsigned short)val;
#endif
          }
        }
        if( bps == 8 ) {
          mivarput( w_fd, var_id, start, count, NC_BYTE, MI_UNSIGNED, &mincdata[vxmin] );
        } else if( bps == 16 ) {
          mivarput( w_fd, var_id, start, count, NC_SHORT, MI_UNSIGNED, &mincdata[vxmin] );
        }

        start[1] += count[1];
      }
    }
#endif
  }
  TIFFClose(r_fp);

  for( j = 0; j < numtiles; j++ ) {
    free( tiledata[j] );
  }
  free( tiledata );
  free( mincdata );

  start[0] = 0;
  mivarput1(w_fd, ncvarid(w_fd, MIimagemin), start, NC_DOUBLE, MI_SIGNED, 
            &vrange[0]);
//          &min);
  start[0] = 0;
  mivarput1(w_fd, ncvarid(w_fd, MIimagemax), start, NC_DOUBLE, MI_SIGNED, 
            &vrange[1]);
//          &max);
      
  miclose(w_fd);
}