/* hauncomp.c - decompression of fits files.
 *
 * Andrew Ducore June 2003
 *
 * Usage: hauncomp <filename> <output directory>
 *
 * filename must not contain any extensions:
 * if the files "lspn0123.hdr", "lspn0123.imq" need to be decompressed,
 * filename should be "lspn0123".
 *
 * files will be saved in the output directory
 *
 * Compression Format: (documented Jun 2007, AMD)
 *
 * Data is interpreted as a series of 16 bit integers (16 bit pixels), and
 * since neighboring pixels tend to be close in value, this allows us to store
 * the difference between the old pixel and the new pixel as a 8 bit integer
 * rather than using the full 16 bits. When the difference is too large,
 * we have to store the enire 16 bits (plus a byte that tells us that those
 * 16 bits aren't two 'difference' values. For this, we use the following 
 * convention:
 * 
 *   If the difference is between -127 and +127, just store the difference
 *   plus 127, as an unsigned 8 bit integer, so a value of 0 means that the
 *   new pixel value is 127 less than the old value, a value of 127 means that
 *   the new pixel value is identical to the old value, and a value of 254
 *   means that the new value is 127 more than the old value.
 *
 *   If the difference is outside of the -127 to +127 range, we store the
 *   entire value, preceeded by a byte with the value 255 (0xFF).
 *
 *   The first value in a compressed file must, therefore, be 0xFF, since
 *   there is no previous value to modify. (In theory, this 0xFF could be
 *   omitted since it is implied by the start of the file, but it isn't.)
 *
 * To decompress,
 *
 *   Iterate through the file as follows:
 *   
 *     If the next byte B is 0xFF, discard it and read the two following
 *     bytes. Save these bytes in a variable V (so that we can use them to
 *     compute new values based on  'difference' bytes), and write these two
 *     bytes out to the file.
 *
 *     Otherwise, take the old value V and add to it 127 minus the byte B's
 *     value. Take this as the new value stored in V, and write the two bytes
 *     in V out to the file.
 *
 */

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

// function prototypes:

int  processHeader          (FILE* headerFile,          FILE* outputFile);
// Parses the header file, and constructs the appropriate header in the
// output file.
// input  [headerFile] handle to file containing header information
//        [outputFile] handle to file where all results are to be written
// output: number of bytes in the compressed image

void processCompressedImage (FILE* compressedImageFile, FILE* outputFile,
			     int   compressedDataLength);
// Takes a compressed image file and decompresses it into the output file.
// input: [compressedImageFile] handle to file containing compressed image
//        [outputFile] handle to file where uncompressed data should be written
//        [compressedDataLength] number of bytes to decompress
// output: none

int main (int argc, char** argv)
{
   // variables
   char* lspnfile;
   char* output_directory;
   char* filename;
   char* longfilename;
   int suffixPos;
   FILE* headerFile;
   FILE* compressedImageFile;
   FILE* outputFile;
   int compressedDataLength;

   // handle arguments (make sure there're 3 of them)
   if (argc != 3)
   {
      printf ("Usage: hauncomp <lspnfile> <output_directory>\n");
      printf ("   lspnfile should not have extensions or directories\n");
      printf ("example:\n");
      printf ("      hauncomp lspn0123 decompressed/\n");
      printf ("   would decompress lspn0123.hdr and lspn0123.imq in the current\n");
      printf ("   directory into decompressed/lspn0123.fits\n");
      exit (1);
   }
   lspnfile         = argv[1];
   output_directory = argv[2];

   // Create a holder for the filename
   suffixPos = strlen (lspnfile);
   filename = malloc (suffixPos + 5);
   strcpy (filename, lspnfile);
   // Add the .hdr suffix to the filename
   strcpy (filename + suffixPos, ".hdr");
   // open header file
   headerFile = fopen (filename, "rb");
   // Replace the .hdr suffix with .imq
   strcpy (filename + suffixPos, ".imq");
   // open compressed image file
   compressedImageFile = fopen (filename, "rb");
   // Create filename for output file
   longfilename = malloc (strlen(output_directory) + suffixPos + 6);
   longfilename[0] = 0;
   strcpy (longfilename, output_directory);
   strcpy (longfilename + strlen(output_directory), filename);
   strcpy (longfilename + strlen(output_directory) + suffixPos, ".fits");
   // open output file
   outputFile = fopen (longfilename, "wb");
   // free memory
   free (filename);
   free (longfilename);

   // process header; close header file.
   compressedDataLength = processHeader (headerFile, outputFile);
   fclose (headerFile);

   // process image
   processCompressedImage (compressedImageFile, outputFile,
			   compressedDataLength);

   // close remaining files.
   fclose (compressedImageFile);
   fclose (outputFile);

   // success.
   return 0;
}

int processHeader (FILE* headerFile,
		   FILE* outputFile)
{
   // variables
   char* inbuffer;
   char* outbuffer;
   char* line;
   int index;
   int outindex;
   int eoh;
   int started;
   int moreMod;
   int compressedDataLength;

   // allocate space for buffers
   inbuffer  = malloc (2881);
   outbuffer = malloc (2881);

   // set output buffer to all blanks.
   memset (outbuffer, ' ', 2880);
   // terminate with null.
   outbuffer[2880] = 0;

   // first, confirm valid primary header
   eoh = 0; // End Of Header flag.
   // Read the file one 2880 byte block at a time.
   while (!eoh && fread (inbuffer, 80, 36, headerFile))
   {
      // 36 80-character entries per 2880 byte read
      for (index = 0; index < 36; ++index)
      {
	 // point line to current entry
	 line = inbuffer + 80*index;
	 // check for "END", set eoh if found
	 if (!strncmp (line, "END", 3))
	 {
	    eoh = 1;
	 }
	 else
	 {
	    // make sure SIMPLE = T
	    if (!strncmp (line, "SIMPLE  ", 8) && line[29] != 'T')
	    {
	       printf ("SIMPLE keyword is not 'T'.\n");
	       exit (1);
	    }
	    // make sure BITPIX = 8
	    else if (!strncmp (line, "BITPIX  ", 8) && line[29] != '8' 
		     && line[30] != ' ')
	    {
	       printf ("BITPIX is not 8.\n");
	       exit (1);
	    }
	    // make sure NAXIS = 0
	    else if (!strncmp (line, "NAXIS   ", 8) && line[29] != '0')
	    {
	       printf ("NAXIS is not 0.\n");
	       exit (1);
	    }
	    // make sure EXTENDED = T
	    else if (!strncmp (line, "EXTEND  ", 8) && line[29] != 'T')
	    {
	       printf ("EXTEND is not 'T'.\n");
	       exit (1);
	    }
	 }
      }
   }

   // begin header
   strcpy (outbuffer, "SIMPLE  =                    T  / BASIC FITS FORMAT");
   // trim off null character
   outbuffer[strlen(outbuffer)] = ' ';

   // continue, intelligently, until end
   outindex = 1;
   eoh = 0;
   started = 0;
   moreMod = 1;
   // while more blocks in header file and not at end.
   while (fread (inbuffer, 80, 36, headerFile) && !eoh)
   {
      // 36 lines per 2880 block
      for (index = 0; index < 36; ++index)
      {
	 // point line to current line.
	 line = inbuffer + 80*index;
	 // if we've processed the compression lines, copy header lines
	 if (started)
	 {
	    if (moreMod) // more modifications, see below
	    {
	       // convert LAXIS* to NAXIS*
	       if (!strncmp (line, "LAXIS", 5))
		  line[0] = 'N';
	       // if just handled NAXIS2, no more modifications
	       if (line[5] == '2')
		  moreMod = 0;
	    }
	    // copy line to output buffer
	    memcpy (outbuffer + outindex*80, line, 80);
	    ++outindex;
	    // if output buffer full, push to disk and clear
	    if (outindex == 36)
	    {
	       fwrite (outbuffer, 80, 36, outputFile);
	       outindex = 0;
	       memset (outbuffer, ' ', 2880);
	    }
	 }
	 else // haven't processed compression lines
	 {
	    // NAXIS1 contains the number of bytes of compressed data
	    if (!strncmp (line, "NAXIS1  ", 8))
	    {
	       line[30] = 0;
	       line = line + 10;
	       compressedDataLength = strtol(line, (char**)NULL, 10);
	    }
	    // LBITPIX is the last compression line; start once we see this.
	    else if (!strncmp (line, "LBITPIX ", 8))
	    {
	       started = 1;
	       memcpy (line, "BITPIX  ", 8);
	       memcpy (outbuffer + outindex*80, line, 80);
	       ++outindex;
	    }
	 }
      }
   }
   // write final block, if not filled already
   if (outindex != 0) fwrite (outbuffer, 80, 36, outputFile);

   // cleanup
   free (inbuffer);
   free (outbuffer);
   return compressedDataLength;
}

void processCompressedImage (FILE* compressedImageFile,
			     FILE* outputFile,
			     int   compressedDataLength)
{
   // Variables
   unsigned char* inbuffer;
   unsigned char* outbuffer;
   int index;
   int outindex;
   unsigned int oldval = 0;
   unsigned char curr[2];
   int partialbyte;
   unsigned int bytesread;
   unsigned int outTotal;
   
   // Initialize buffers
   inbuffer  = malloc(2881);
   outbuffer = malloc(2881);

   outTotal = 0;

   bytesread = 0;
   partialbyte = 0;
   outindex = 0;
   // while we can read more of the file.
   while (fread (inbuffer, 80, 36, compressedImageFile))
   {
      index = 0;
      // if last block had an incomplete piece of data
      if (partialbyte)
      {
	 // just 0xFF?
	 if (partialbyte == 1)
	 {
	    // get the 2 bytes from the beginning of the buffer
	    curr[0] = inbuffer[0];
	    curr[1] = inbuffer[1];
	    index = 2;
	    bytesread+=2;
	 }
	 else // 0xFF + 1st byte
	 {
	    // grab second byte, that's all.
	    curr[1] = inbuffer[0];
	    ++bytesread;
	    index = 1;
	 }
	 partialbyte = 0;
	 // if outputing will split up the 2-byte value:
	 if (outindex == 2879)
	 {
	    // put 1st byte at end.
	    outbuffer[2879] = curr[0];
	    // output buffer
	    fwrite (outbuffer, 80, 36, outputFile);
	    // put 2nd byte at new beginning
	    outbuffer[0] = curr[1];
	    outindex = 1;
	    outTotal+=2;
	 }
	 else if (outindex == 2878)
	 {
	    // output will fill the buffer, fill and write.
	    outbuffer[2878] = curr[0];
	    outbuffer[2879] = curr[1];
	    fwrite (outbuffer, 80, 36, outputFile);
	    outindex = 0;
	    outTotal+=2;
	 }
	 else
	 {
	    // put values in output buffer.
	    outbuffer[outindex++] = curr[0];
	    outbuffer[outindex++] = curr[1];
	    outTotal+=2;
	 }
	 // compute new pixel value
	 oldval = curr[0] * 256 + curr[1];
      }
      else index = 0;
      // for the rest of the buffer:
      for (;index < 2880 && bytesread < compressedDataLength; ++index)
      {
	 if (inbuffer[index] != 0xFF) // just a difference byte
	 {
	    // compute new value.
	    oldval += (inbuffer[index] - 127);
	    curr[0] = oldval >> 8;
	    curr[1] = oldval % 256;
	    ++bytesread;
	 }
	 else // a new 2-byte value
	 {
	    if (index < 2878)
	    {
	       // space in input buffer, grab all and move on.
	       curr[0] = inbuffer[++index];
	       curr[1] = inbuffer[++index];
	       oldval = (int)(curr[0]) * 256 + curr[1];
	       bytesread+=3;
	    }
	    else if (index == 2878)
	    {
	       // only one byte available, save and move on to next 2880 block
	       curr[0] = inbuffer[++index];
	       partialbyte = 2; // 0xFF + 1st byte value
	       bytesread+=2;
	       break;
	    }
	    else
	    {
	       // just the 0xFF in this buffer, move on to next for
	       partialbyte = 1; // just 0xFF
	       ++bytesread;
	       break;
	    }
	 }
	 if (outindex < 2878)
	 {
	    // There is space in output buffer
	    outbuffer[outindex++] = curr[0];
	    outbuffer[outindex++] = curr[1];
	    outTotal+=2;
	 }
	 else if (outindex == 2878)
	 {
            // fill output buffer and write
	    outbuffer[2878] = curr[0];
	    outbuffer[2879] = curr[1];
	    fwrite (outbuffer, 80, 36, outputFile);
	    outindex = 0;
	    outTotal+=2;
	 }
	 else
	 {
	    // just room for one byte, put it in, write the output,
	    // and put the remaining byte at the beginning of the new buffer
	    outbuffer[2879] = curr[0];
	    fwrite (outbuffer, 80, 36, outputFile);
	    outbuffer[0] = curr[1];
	    outindex = 1;
	    outTotal+=2;
	 }
      }
      // if at end of data, quit
      if (bytesread >= compressedDataLength) break;
   }
   // if output buffer isn't empty:
   if (outindex != 0)
   {
      // zero-fill the rest of the buffer
      memset (outbuffer+outindex, 0, 2880-outindex);
      // write the buffer to disk
      fwrite (outbuffer, 80, 36, outputFile);
   }

   // free dynamic memory allocations
   free (inbuffer);
   free (outbuffer);
}