/*
 * WTIFC.MEX
 * 
 * This is a mex interface to the Sam Leffler's LIBTIFF library which
 * will allow images to be written in several of the many variants of
 * the TIFF image file format.   
 * 
 * The syntaxes are:
 *   
 *     wtifc(CMYK, [], filename, compression,description, resolution, overwrite);
 *     wtifc(A, [], filename, compression, description, resolution, overwrite, colorspace);
 *     wtifc(PACKED_RGB, [], filename, compression, description, resolution, overwrite);
 *     wtifc(GRAY, [], filename, compression, description, resolution, overwrite);
 *     wtifc(X, map, filename, compression, description, resolution, overwrite);
 *   
 * CMYK is a mxnx4 uint8 or uint16 array.
 *
 * A is an mxnx3 uint8 or uint16 array; colorspace can be 'rgb', 'icclab',
 * or 'cielab'.
 * 
 * PACKED_RGB is an mxn uint32 array [with each 4 bytes being RGBA, were A is 
 * currently garbage (maybe Alpha in the future)] containing a 24-bit image to
 * be stored in the TIFF file filename.  This form was implemented for HG to use.
 * 
 * GRAY is a mxn uint8 array containing a grayscale image to
 * be stored in the TIFF file filename.    
 * 
 * X is a uint8 array of indices into the colormap map.
 * map is an M-by-3 uint8 array, where M <= 256.
 *
 * The compression string allows a choice between the following
 * compression schemes:
 * "packbits"   Runlength encoding (default)
 * "ccitt"	CCITT Group 3 1-D Modified Huffman RLE alg.
 * "packbits"	Macintosh PackBits algorithm
 * "thunder"	ThunderScan 4-bit RLE algorithm
 * "next"	NeXT 2-bit RLE algorithm
 * "none"	No compression
 *
 * The resolution is the same for both X and Y. 
 * 
 * The overwrite flag indicates whether we should overwrite an existing
 * file (if overwrite is 1) or append to an existing file (if overwrite
 * is 0).
 *
 * 
 * KNOWN BUGS:
 * -----------
 *    Thunderscan compression isn't correctly implemented yet. 
 *    NeXT compression isn't correctly implemented yet.
 *    
 * ENHANCEMENTS UNDER CONSIDERATION:
 * ---------------------------------
 *    Add CCITT Fax3 and Fax4 writing support 
 *      (probably very simple since CCITTRLE is already implemented.)
 *    Allow tifwrite to write tiled tiff images.
 *
 * 
 * Sam  Leffler's LIBTIFF library, version 3.4 is available from:
 * ftp://ftp.sgi.com/graphics/tiff/tiff-v3.4-tar.gz
 *
 * Chris Griffin, June 1996
 * Copyright 1984-2001 The MathWorks, Inc. 
 * $Revision: 1.22 $  $Date: 2003/07/17 21:14:32 $
 */

#include "mex.h"
#include <stdio.h>
#include <string.h>
#include "tiffio.h"

/* Different types of images we can write: */
#define BINARY_IMG 0
#define GRAY_IMG   1
#define INDEX_IMG  2
#define RGB_IMG    3
#define RGBA_IMG   4
#define PACKED_RGB_IMG   5
#define CMYK_IMG   6
#define CIELAB_IMG 7
#define ICCLAB_IMG 8

#define STRLEN 1024

static void WriteBitsFromLogical(TIFF *tif, mxLogical *data, int w, int h, int rps);
static void WriteBytesFromUint8(TIFF *tif, uint8_T *data, int w, int h, int rps);
static void WriteBytesFromUint16(TIFF *tif, uint16_T *data, int w, int h, int rps);
static void WriteRGBFromUint8(TIFF *tif, uint8_T *data, int w, int h, int rps);
static void WriteRGBFromUint16(TIFF *tif, uint16_T *data, int w, int h, int rps);
static void WriteRGBFromUint32(TIFF *tif, uint32_T *data, int w, int h, int rps);
static void WriteCMYKFromUint8(TIFF *tif, uint8_T *data, int w, int h, int rps);
static void WriteCMYKFromUint16(TIFF *tif, uint16_T *data, int w, int h, int rps);
static void setColormap(TIFF *tif, const mxArray *cmap, int dataClass);
static void ErrHandler(const char *module, const char *fmt, va_list ap);
static void WarnHandler(const char *module, const char *fmt, va_list ap);

void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[]) 
{
    
    TIFF *tif;
    const mxArray *inputArray, *mapArray;
    const mxArray *fnameArray, *compArray, *desArray, *resArray;
    const mxArray *colorspaceArray = NULL;
    uint8_T *uint8Data;
    uint16_T *uint16Data;
    int ndims;                    /* Number of dimensions for input array */
    const int *size;
    char *filename;
    char *description;
    long stringlen;
    int w,h;                    /* width, height, and number of planes */
    int comp;                     /* Compression */
    int overwrite=1;            /* Will we overwrite or append? Default is overwrite*/
    char compStr[64];
    int imType, inClass;
    int rowsPerStrip;
    unsigned short resolutionUnit = 2;   /* Always */
    unsigned short samplesPerPixel = 1;  /* Default */
    float xres, yres;
    double *resPr;
    
    TIFFSetErrorHandler(ErrHandler);
    TIFFSetWarningHandler(WarnHandler);
    
    if (nrhs < 7 || nrhs > 8)
    {
        mexErrMsgTxt("Illegal number of input arguments.");
    }      
    
    inputArray = prhs[0];         /* First argument is the image data */
    mapArray = prhs[1];
    fnameArray = prhs[2];
    compArray = prhs[3];
    desArray = prhs[4];
    resArray = prhs[5];
    
    if (!mxIsDouble(resArray)) {
        mexErrMsgTxt("Resolution argument must be double.");
    }
    resPr = mxGetPr( resArray ); 
    switch (mxGetNumberOfElements(resArray))
    {
    case 1:
        xres = yres = (float) resPr[0];
        break;
        
    case 2:
        xres = (float) resPr[0];
        yres = (float) resPr[1];
        break;
        
    default:
        mexErrMsgTxt("Resolution must be a scalar or a 2-element vector.");
    }

    overwrite = mxGetScalar( prhs[6] ); 

    if (nrhs >= 8)
    {
        colorspaceArray = prhs[7];
    }

    inClass = mxGetClassID(inputArray); 
    ndims = mxGetNumberOfDimensions(inputArray);
    size = mxGetDimensions(inputArray);
    h = size[0];
    w = size[1];

    if (ndims == 3)
    {
        if ( (inClass != mxUINT8_CLASS) && (inClass != mxUINT16_CLASS) )
        {
            mexErrMsgTxt("Three-dimensional image array must be uint8 or uint16.");
        }

        if (size[2] == 4)
        {
            imType = CMYK_IMG;
        }
        else if (size[2] == 3)
        {
            if (colorspaceArray == NULL)
            {
                imType = RGB_IMG;
            }
            else
            {
                char str[STRLEN];
                if (mxGetString(colorspaceArray, str, STRLEN) != 0)
                {
                    mexErrMsgTxt("Invalid colorspace string.");
                }
                if (strcmp(str, "rgb") == 0)
                {
                    imType = RGB_IMG;
                }
                else if (strcmp(str, "cielab") == 0)
                {
                    imType = CIELAB_IMG;
                }
                else if (strcmp(str, "icclab") == 0)
                {
                    imType = ICCLAB_IMG;
                }
                else
                {
                    mexErrMsgTxt("Invalid colorspace string.");
                }
            }
        }
        else
        {
            mexErrMsgTxt("Invalid image array; third dimension must be 1, 3, or 4.");
        }
    }
    else if (ndims == 2)
    {
        if ( (inClass != mxLOGICAL_CLASS) && (inClass != mxUINT8_CLASS) &&
             (inClass != mxUINT16_CLASS)  && (inClass != mxUINT32_CLASS) )
        {
            /* uint32 case is undocumented, so we don't mention it in the error message. */
            mexErrMsgTxt("Two-dimensional image array must be logical, uint8, or uint16.");
        }

        if (inClass == mxUINT32_CLASS)
        {
            /* Packed, zbuffer-type truecolor data; used by HG */
            imType = PACKED_RGB_IMG;
        }
        else if (! mxIsEmpty(mapArray))
        {
            imType = INDEX_IMG;
        }
        else if (inClass == mxLOGICAL_CLASS)
        {
            imType = BINARY_IMG;
        }
        else
        {
            imType = GRAY_IMG;
        }
    }
    else
    {
        mexErrMsgTxt("Image array must be two- or three-dimensional.");
    }

    if (!mxIsChar(compArray))
    {
        mexErrMsgTxt("COMPRESSION must be a string");
    }
    mxGetString(compArray, compStr, 64);
    
    if(!strncmp(compStr, "packbits",64))
        comp = COMPRESSION_PACKBITS; 
    else if(!strncmp(compStr, "ccitt",64))
        comp = COMPRESSION_CCITTRLE;
    else if(!strncmp(compStr, "thunder",64))
        comp = COMPRESSION_THUNDERSCAN;
    else if(!strncmp(compStr, "next",64))
        comp = COMPRESSION_NEXT;
    else if(!strncmp(compStr, "fax3",64))
        comp = COMPRESSION_CCITTFAX3;
    else if(!strncmp(compStr, "fax4",64))
        comp = COMPRESSION_CCITTFAX4;
    else if(!strncmp(compStr, "none",64))
        comp = COMPRESSION_NONE;
    else
        mexErrMsgTxt("Invalid compression scheme requested.");
    
    if((comp == COMPRESSION_CCITTRLE || 
        comp == COMPRESSION_CCITTFAX3 ||
        comp == COMPRESSION_CCITTFAX4) && imType != BINARY_IMG)
        mexErrMsgTxt("CCITT compression is only valid for binary images.");


    /* Open Tiff file for writing */

    if (!mxIsChar(fnameArray))
    {
        mexErrMsgTxt("FILENAME must be a string");
    }
    stringlen = mxGetM(fnameArray) * mxGetN(fnameArray) * sizeof(mxChar) + 1;
    filename = (char *) mxCalloc(stringlen, sizeof(*filename));
    mxGetString(fnameArray, filename, stringlen);

    if(overwrite) {
        if ((tif = TIFFOpen(filename, "w")) == NULL) {
            mexErrMsgTxt("Couldn't open TIFF file for writing.");
        }
    }
    else {
        if ((tif = TIFFOpen(filename, "a")) == NULL) {
            mexErrMsgTxt("Couldn't open TIFF file for append.");
        }
    }
    mxFree((void *) filename);


    /* Is there a description? */
    if (!mxIsEmpty(desArray))
    {
        if (!mxIsChar(desArray))
        {
            mexErrMsgTxt("DESCRIPTION must be a string");
        }
        stringlen = mxGetM(desArray) * mxGetN(desArray) * sizeof(mxChar) + 1;
        description = (char *) mxCalloc(stringlen, sizeof(*description));
        mxGetString(desArray, description, stringlen);
        TIFFSetField(tif, TIFFTAG_IMAGEDESCRIPTION, description);
        mxFree((void *) description);
    }

    TIFFSetField(tif, TIFFTAG_IMAGEWIDTH, w);
    TIFFSetField(tif, TIFFTAG_IMAGELENGTH, h);
    TIFFSetField(tif, TIFFTAG_COMPRESSION, comp);
    TIFFSetField(tif, TIFFTAG_RESOLUTIONUNIT, resolutionUnit);    
    TIFFSetField(tif, TIFFTAG_XRESOLUTION, xres);
    TIFFSetField(tif, TIFFTAG_YRESOLUTION, yres);
    
    
    if (comp == COMPRESSION_CCITTFAX3)
        TIFFSetField(tif, TIFFTAG_GROUP3OPTIONS,
                     GROUP3OPT_2DENCODING+GROUP3OPT_FILLBITS);
    
    TIFFSetField(tif, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
    TIFFSetField(tif, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT);

    /* write the image data */
    if (imType == RGB_IMG || imType == PACKED_RGB_IMG || imType == CIELAB_IMG || imType == ICCLAB_IMG)
    {
        if(inClass == mxUINT16_CLASS) 
            TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE,   16);
        else   /* For mxUINT8_CLASS and mxUINT32_CLASS */
            TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE,   8);

        if (imType == CIELAB_IMG)
        {
            TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_CIELAB);
        }
        else if (imType == ICCLAB_IMG)
        {
            /* Reference: "Adobe Photoshop TIFF Technical Notes," March 22, 2002, www.adobe.com */
            TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, 9);
        }
        else
        {
            TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
        }

        samplesPerPixel = 3;
        TIFFSetField(tif, TIFFTAG_SAMPLESPERPIXEL, samplesPerPixel);

        /* Make each uncompressed strip fit in 8k if possible */
        if(w==0)
            rowsPerStrip = 1;
        else 
            rowsPerStrip = 8192 / (w * 3);
        
        /* Protect ourselves from empty strips */
        if (rowsPerStrip == 0)
            rowsPerStrip = 1;

        TIFFSetField(tif, TIFFTAG_ROWSPERSTRIP, rowsPerStrip);
        
        if ( imType == PACKED_RGB_IMG )
            WriteRGBFromUint32(tif, (uint32_T *) mxGetData(inputArray), w, h, rowsPerStrip);
        else if ( inClass == mxUINT8_CLASS )  /* It's RGB_IMG, ICCLAB_IMG, or CIELAB_IMG */
            WriteRGBFromUint8(tif, (uint8_T *) mxGetData(inputArray), w, h, rowsPerStrip);
        else if ( inClass == mxUINT16_CLASS ) /* It's RGB_IMG, ICCLAB_IMG, or CIELAB_IMG */
            WriteRGBFromUint16(tif, (uint16_T *) mxGetData(inputArray), w, h, rowsPerStrip);
        
    } 
    else if (imType == CMYK_IMG)
    {
        if(inClass == mxUINT16_CLASS) 
        {
            TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE,   16);
        }
        else
        {
            TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE,   8);
        }

        TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_SEPARATED);

        samplesPerPixel = 4;
        TIFFSetField(tif, TIFFTAG_SAMPLESPERPIXEL, samplesPerPixel);

        /* Make each uncompressed strip fit in 8k if possible */
        if(w==0)
            rowsPerStrip = 1;
        else 
            rowsPerStrip = 8192 / (w * 4);
        
        /* Protect ourselves from empty strips */
        if (rowsPerStrip == 0)
            rowsPerStrip = 1;

        TIFFSetField(tif, TIFFTAG_ROWSPERSTRIP, rowsPerStrip);

        if (inClass == mxUINT8_CLASS)
        {
            WriteCMYKFromUint8(tif, (uint8_T *) mxGetData(inputArray), w, h, rowsPerStrip);
        }
        else
        {
            WriteCMYKFromUint16(tif, (uint16_T *) mxGetData(inputArray), w, h, rowsPerStrip);
        }

    }
    else          /* Indexed, Gray or Binary image */
    {
        
        /* Make each uncompressed strip fit in 8k if possible */
        if(w==0)      /* Catch degenerate case - don't divide by 0 ! */
            rowsPerStrip = 1;
        else
            rowsPerStrip = 8192/w;
        
        if (rowsPerStrip == 0)
        {
            rowsPerStrip = 1;
        }
        TIFFSetField(tif, TIFFTAG_ROWSPERSTRIP, rowsPerStrip);
        TIFFSetField(tif, TIFFTAG_SAMPLESPERPIXEL, samplesPerPixel);        

        if(imType == INDEX_IMG)
        {
            if(inClass == mxUINT16_CLASS) 
                TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE,   16);
            else
                TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE,   8);

            TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_PALETTE);
            setColormap(tif, mapArray, inClass);
        }
        else if(imType == BINARY_IMG)
        {
            TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE,   1);
            TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISWHITE);
        }
        else /* Grayscale or Binary without CCITT compression */
        {
            if(inClass == mxUINT16_CLASS) 
                TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE,   16);
            else
                TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE,   8);

            TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK);
        }
        
        if ( inClass == mxLOGICAL_CLASS ) {
	    mxLogical *logicalData = mxGetLogicals(inputArray);
	    WriteBitsFromLogical(tif, logicalData, w, h, rowsPerStrip);
	}
        else if ( inClass == mxUINT8_CLASS ) {
            uint8Data  = (uint8_T *) mxGetData(inputArray);
            WriteBytesFromUint8(tif, uint8Data, w, h, rowsPerStrip);
        }
        else if ( inClass == mxUINT16_CLASS ) {
            uint16Data = (uint16_T *) mxGetData(inputArray);
            WriteBytesFromUint16(tif, uint16Data, w, h, rowsPerStrip);
        }
    }
  
    /*
     * Clean up
     */

    TIFFClose(tif);
    return;		
}


/*
** For writing RGB images.
*/

static void 
WriteRGBFromUint8(TIFF *tif, uint8_T *data, int w, int h, int rps)
{
    
    uint8_T *buf;
    int i,j,p;
    int numStrips;
    int strip;
    int col, row;

    /* Allocate buffer for image write */
    buf = (uint8_T *) mxCalloc(w*rps*3, sizeof(uint8_T));

    numStrips = (h + rps - 1) / rps;

    for (strip = 0, row = 0; strip < numStrips; strip++)
    {
        /* Fill the strip buffer */
        for (p = 0; (p < rps) && (row < h); p++, row++)
        {
            for (col = 0; col < w; col++)
            {
                j = col + (p*w);    /* index into buf */
                i = (col*h) + row;  /* index into data */

                buf[j*3]   = data[i];         /* Red Pixel */
                buf[j*3+1] = data[i+(w*h)];   /* Green Pixel */
                buf[j*3+2] = data[i+(2*w*h)]; /* Blue Pixel */
            }
        }
    
        /* Write the strip buffer */
        TIFFWriteEncodedStrip(tif, strip, buf, w*p*3*sizeof(uint8));
    }

    mxFree((void *) buf);
}


static void 
WriteRGBFromUint16(TIFF *tif, uint16_T *data, int w, int h, int rps)
{
    
    uint16_T *buf;
    int i,j,p;
    int numStrips;
    int strip;
    int col, row;

    /* Allocate buffer for image write */
    buf = (uint16_T *) mxCalloc(w*rps*3, sizeof(uint16_T));

    numStrips = (h + rps - 1) / rps;

    for (strip = 0, row = 0; strip < numStrips; strip++)
    {
        /* Fill the strip buffer */
        for (p = 0; (p < rps) && (row < h); p++, row++)
        {
            for (col = 0; col < w; col++)
            {
                j = col + (p*w);    /* index into buf */
                i = (col*h) + row;  /* index into data */

                buf[j*3]   = data[i];         /* Red Pixel */
                buf[j*3+1] = data[i+(w*h)];   /* Green Pixel */
                buf[j*3+2] = data[i+(2*w*h)]; /* Blue Pixel */
            }
        }
    
        /* Write the strip buffer */
        TIFFWriteEncodedStrip(tif, strip, buf, w*p*3*sizeof(uint16_T));
    }

    mxFree((void *) buf);
}


static void 
WriteRGBFromUint32(TIFF *tif, uint32_T *data, int w, int h, int rps)
{
    
    uint8_T *buf, *data8;
    int j,p;
    int numStrips;
    int strip;
    int col, row;

    /* Index each byte of the 4 byte uint32 as packed [ Alpha(garbage) R G B ] */
    data8 = (uint8_T *)data;

    /* Allocate buffer for image write */
    buf = (uint8_T *) mxCalloc(w*rps*3, sizeof(uint8_T));

    numStrips = (h + rps - 1) / rps;

    for (strip = 0, row = 0; strip < numStrips; strip++)
    {
        /* Fill the strip buffer */
        for (p = 0; (p < rps) && (row < h); p++, row++)
        {
            for (col = 0; col < w; col++)
            {
                j = col + (p*w);      /* index into buf (why not just buf++?) */
                /* Shifts are hardcoded TCPI offsets of IMSAVE driver */
                buf[j*3]   = (*data >> 16) & (uint32_T)255; /* Red Pixel */
                buf[j*3+1] = (*data >>  8) & (uint32_T)255; /* Green Pixel */
                buf[j*3+2] =  *data        & (uint32_T)255; /* Blue Pixel */
                data++;
            }
        }
    
        /* Write the strip buffer */
        TIFFWriteEncodedStrip(tif, strip, buf, w*p*3);
    }

    mxFree((void *) buf);
}



static void 
WriteCMYKFromUint8(TIFF *tif, uint8_T *data, int w, int h, int rps)
{
    
    uint8_T *buf;
    int i,j,p;
    int numStrips;
    int strip;
    int col, row;

    /* Allocate buffer for image write */
    buf = (uint8_T *) mxCalloc(w*rps*4, sizeof(uint8_T));

    numStrips = (h + rps - 1) / rps;

    for (strip = 0, row = 0; strip < numStrips; strip++)
    {
        /* Fill the strip buffer */
        for (p = 0; (p < rps) && (row < h); p++, row++)
        {
            for (col = 0; col < w; col++)
            {
                j = col + (p*w);    /* index into buf */
                i = (col*h) + row;  /* index into data */

                buf[j*4]   = data[i];         /* Cyan Pixel */
                buf[j*4+1] = data[i+(w*h)];   /* Magenta Pixel */
                buf[j*4+2] = data[i+(2*w*h)]; /* Yellow Pixel */
                buf[j*4+3] = data[i+(3*w*h)]; /* Black Pixel */
            }
        }
    
        /* Write the strip buffer */
        TIFFWriteEncodedStrip(tif, strip, buf, w*p*4*sizeof(uint8));
    }

    mxFree((void *) buf);
}


static void 
WriteCMYKFromUint16(TIFF *tif, uint16_T *data, int w, int h, int rps)
{
    
    uint16_T *buf;
    int i,j,p;
    int numStrips;
    int strip;
    int col, row;

    /* Allocate buffer for image write */
    buf = (uint16_T *) mxCalloc(w*rps*4, sizeof(uint16_T));

    numStrips = (h + rps - 1) / rps;

    for (strip = 0, row = 0; strip < numStrips; strip++)
    {
        /* Fill the strip buffer */
        for (p = 0; (p < rps) && (row < h); p++, row++)
        {
            for (col = 0; col < w; col++)
            {
                j = col + (p*w);    /* index into buf */
                i = (col*h) + row;  /* index into data */

                buf[j*4]   = data[i];         /* Cyan Pixel */
                buf[j*4+1] = data[i+(w*h)];   /* Magenta Pixel */
                buf[j*4+2] = data[i+(2*w*h)]; /* Yellow Pixel */
                buf[j*4+3] = data[i+(3*w*h)]; /* Black Pixel */
            }
        }
    
        /* Write the strip buffer */
        TIFFWriteEncodedStrip(tif, strip, buf, w*p*4*sizeof(uint16_T));
    }

    mxFree((void *) buf);
}


/*
** For writing Gray, Binary, or Indexed images.
*/

static void 
WriteBytesFromUint8(TIFF *tif, uint8_T *data, int w, int h, int rps)
{
    
    uint8_T *buf;
    int i,j,p;
    int numStrips;
    int strip;
    int col, row;

    /* Allocate buffer for image write */
    buf = (uint8_T *) mxCalloc(w*rps, sizeof(uint8_T));

    numStrips = (h + rps - 1) / rps;

    for (strip = 0, row = 0; strip < numStrips; strip++)
    {
        /* Fill the strip buffer */
        for (p = 0; (p < rps) && (row < h); p++, row++)
        {
            for (col = 0; col < w; col++)
            {
                j = col + (p*w);    /* index into buf */
                i = (col*h) + row;  /* index into data */

                buf[j] = data[i];
            }
        }
    
        /* Write the strip buffer */
        TIFFWriteEncodedStrip(tif, strip, buf, w*p*sizeof(uint8));
    }

    mxFree(buf);
}


static void 
WriteBytesFromUint16(TIFF *tif, uint16_T *data, int w, int h, int rps)
{
    
    uint16_T *buf;
    int i,j,p;
    int numStrips;
    int strip;
    int col, row;

    /* Allocate buffer for image write */
    buf = (uint16_T *) mxCalloc(w*rps, sizeof(uint16_T));

    numStrips = (h + rps - 1) / rps;

    for (strip = 0, row = 0; strip < numStrips; strip++)
    {
        /* Fill the strip buffer */
        for (p = 0; (p < rps) && (row < h); p++, row++)
        {
            for (col = 0; col < w; col++)
            {
                j = col + (p*w);    /* index into buf */
                i = (col*h) + row;  /* index into data */

                buf[j] = data[i];
            }
        }
    
        /* Write the strip buffer */
        TIFFWriteEncodedStrip(tif, strip, buf, w*p*sizeof(uint16_T));
    }

    mxFree(buf);
}

/*
 * WriteBitsFromLogical
 * Write from a binary image stored in mxLogicals.
 *
 * Write the image data with White-is-Zero (PHOTOMETRIX_MINISWHITE).  This
 * is non-intuitive, but it is the norm for binary TIFFS.  Many (supposed)
 * Tiff readers will not read the data correctly if it is written with
 * PHOTOMETRIC_MINISBLACK .
 */

static void 
WriteBitsFromLogical(TIFF *tif, mxLogical *data, int w, int h, int rps)
{
    
    uint8_T *buf;
    int i,j,k,p;
    int numStrips;
    int strip;
    uint8_T pixel, byte;
    int col, row;
    int stripBytes;

    /* Allocate buffer for image write */
    stripBytes = ((w+7)/8) * rps;
    buf = (uint8_T *) mxCalloc(stripBytes, sizeof(uint8_T));

    numStrips = (h + rps - 1) / rps;

    for (strip = 0, row = 0; strip < numStrips; strip++)
    {
        j = 0;    /* index into buf */
        k = 0;  /* bit within a byte */
        byte = 0;
        /* Fill the strip buffer */
        for (p = 0; (p < rps) && (row < h); p++, row++)
        {
            for (col = 0; col < w; col++)
            {
                i = (col*h) + row;          /* index into data */

                pixel = (data[i]==0) ? 1 : 0;  /* 0 is white */
                byte |= (pixel << (7-k));

                k++;
                if ((k == 8) || (col == (w-1)))
                {
                    /* either we've filled a byte or we've */
                    /* reached the end of a row.  Stuff the */
                    /* byte into the buffer. */
                    buf[j] = byte;
                    j++;
                    k = 0;
                    byte = 0;
                }
            }
        }
    
        TIFFWriteEncodedStrip(tif, strip, buf, j);
    }
    mxFree((void *) buf);
}



/*
 * setColormap
 * 
 * Put the colormap into the TIFF structure for the write operation.
 */

static void setColormap(TIFF *tif, const mxArray *cmap, int dataClass)
{
    uint16_T *red,*green,*blue;
    int i;
    const int *size;
    int rows, cols, ndims;
    mxClassID mapClass;
    uint8_T *u8_data;
    uint16_T *u16_data;

    red = (uint16_T *) mxCalloc(65536, sizeof(*red));
    green = (uint16_T *) mxCalloc(65536, sizeof(*green));
    blue = (uint16_T *) mxCalloc(65536, sizeof(*blue));
    
    for(i=0; i<256; i++)        /* Clear colormap */
    {
        red[i] = 0;
        green[i] = 0;
        blue[i] = 0;
    }
    
    mapClass = mxGetClassID(cmap);
    ndims = mxGetNumberOfDimensions(cmap);
    if(ndims != 2)
        mexErrMsgTxt("Invalid colormap, must be 2-D.");
    
    size = mxGetDimensions(cmap);
    rows = size[0];
    cols = size[1];

    if(cols!=3)
        mexErrMsgTxt("Invalid colormap, must be n X 3.");

    if( dataClass==mxUINT8_CLASS && rows>256)
        mexErrMsgTxt("Invalid colormap for 8-bit image, must be n X 3 (n<=256).");

    if( dataClass==mxUINT16_CLASS && rows>65536)
        mexErrMsgTxt("Invalid colormap for 16-bit image, must be n X 3 (n<=65536).");


    if (mapClass==mxUINT8_CLASS) {
        /* Multiply data by 65535/255 = 257 */
        u8_data = (uint8_T *) mxGetData(cmap);
        for(i=0; i<rows; i++)
        {
            red[i]   = (uint16_T) (257*u8_data[i]);
            green[i] = (uint16_T) (257*u8_data[i+rows]);
            blue[i]  = (uint16_T) (257*u8_data[i+(rows*2)]);
        }            
    }
    else if(mapClass==mxUINT16_CLASS) {
        u16_data = (uint16_T *) mxGetData(cmap);
        for(i=0; i<rows; i++)
        {
            red[i]   = (uint16_T) u16_data[i];
            green[i] = (uint16_T) u16_data[i+rows];
            blue[i]  = (uint16_T) u16_data[i+(rows*2)];
        }            
    } else
        mexErrMsgTxt("Class of colormap must be uint8 or uint16");

    TIFFSetField(tif, TIFFTAG_COLORMAP, red, green, blue);
}

/*******************************************/
static void ErrHandler(const char *module, 
                       const char *fmt, va_list ap)
{
  char buf[2048];
  char *cp = buf;

  if (module != NULL) {
    sprintf(cp, "%s: ", module);
    cp = (char *) strchr(cp, '\0');
  }

  vsprintf(cp, fmt, ap);
  strcat(cp, ".");

  mexErrMsgTxt(buf);
}


/*******************************************/
static void WarnHandler(const char *module, 
                       const char *fmt, va_list ap)
{
  char buf[2048];
  char *cp = buf;

  if (module != NULL) {
    sprintf(cp, "%s: ", module);
    cp = (char *) strchr(cp, '\0');
  }

  vsprintf(cp, fmt, ap);
  strcat(cp, ".");

  mexWarnMsgTxt(buf);
}

