/* [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved. */
#ifdef NEVERAGAIN
mexwritegif.c
plot(1:10);[xx,cc]=getframe(gcf);writegif('temp.gif',xx,cc);
#endif
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <mex.h>
#define STRLEN 128

#define BITS     12                     /* largest code size */
#define THELIMIT 4096                   /* NEVER generate this */
#define HSIZE    5003                   /* hash table size */
#define SHIFT    4                      /* shift for hashing */

struct gl {
   int      BitsPixel,		       /* number of bits per pixel */
   IniCodeSize,			       /* initial number of bits per code */
   CurCodeSize,			       /* current number of bits per code */
   CurMaxCode,			       /* maximum code, given CurCodeSize */
   ClearCode,			       /* reset code */
   EOFCode,			       /* end of file code */
   FreeCode;			       /* first unused entry */
   long      Nbyte;
   long     HashTab [HSIZE];        /* hash table */
   int      CodeTab [HSIZE];        /* code table */ 
   unsigned long cur_accum;
   int           cur_bits;
   int           a_count;
   char          accum[256];
   unsigned long *masks;
};


static void     output(int code, FILE *File, struct gl *glob);
static void     char_out(int c, FILE *File, struct gl *glob);
static void     char_flush(FILE *File, struct gl *glob);
static void     put_short(short word, FILE *File);

/***********************************************************************
 *                                                                     *
 * Name: GIFencode                                   Date:    02.10.92 *
 * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
 *                                                                     *
 * Function: Encode an image to GIF format                             *
 *                                                                     *
 * The Graphics Interchange Format(c) is the Copyright property of     *
 * CompuServe Incorporated. GIF(sm) is a Service Mark property of      *
 * CompuServe Incorporated.                                            *
 *                                                                     *
 * Input: Width      - image width  (must be >= 8)                     *
 *        Height     - image height (must be >= 8)                     *
 *        Ncol       - number of colors                                *
 *        R[]        - red components                                  *
 *        G[]        - green components                                *
 *        B[]        - blue components                                 *
 *        ScLine[]   - array for scan line (byte per pixel)            *
 *        get_scline - user routine to read scan line:                 *
 *                       get_scline(y, Width, ScLine)                  *
 *        pb         - user routine for "put_byte": pb(b)              *
 *                                                                     *
 * Return: size of GIF                                                 *
 *                                                                     *
 ***********************************************************************/

/* mapping of global variables into structure */
#define BitsPixel   glob->BitsPixel
#define IniCodeSize glob->IniCodeSize
#define CurCodeSize glob->CurCodeSize
#define CurMaxCode  glob->CurMaxCode
#define ClearCode   glob->ClearCode
#define EOFCode     glob->EOFCode
#define FreeCode    glob->FreeCode
#define Nbyte       glob->Nbyte
#define HashTab     glob->HashTab
#define CodeTab     glob->CodeTab
#define cur_accum   glob->cur_accum
#define cur_bits    glob->cur_bits
#define a_count     glob->a_count
#define accum       glob->accum
#define masks       glob->masks
void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{
   struct gl glroot;
   struct gl *glob=&glroot;
   FILE *File;
   double *ScLine;
   int Width,Height,Ncol;
   double *R,*G,*B;
   char info_str[STRLEN],fname[STRLEN];
   double *color;
   int colorn, colorm;
   int maxVal = 255;
   unsigned long mask[] = { 0x0000, 
      0x0001, 0x0003, 0x0007, 0x000F,
      0x001F, 0x003F, 0x007F, 0x00FF,
      0x01FF, 0x03FF, 0x07FF, 0x0FFF,
      0x1FFF, 0x3FFF, 0x7FFF, 0xFFFF };
  long          CodeK;
  int           ncol, i, x, disp, Code, K;

/* initialisations of global */
   masks     = mask;		       /* Assign masks to pointer */
   a_count   = 0;
   cur_accum = 0;
   cur_bits  = 0;

   if(nrhs < 3)
     mexErrMsgTxt("WriteGif(filename, data, colormap); 3 param required");
   else if (mxIsChar(prhs[0]))
     mxGetString(prhs[0],fname,STRLEN);
   else
     mexErrMsgTxt("WriteGif: First argument is a string");

/* open the file */
   if((File = fopen(fname,"wb")) == NULL)  {
      sprintf(info_str,"writeppm : File [%s] not Opened",fname);
      mexErrMsgTxt(info_str);
   }
   Width  = mxGetM(prhs[1]);
   Height = mxGetN(prhs[1]);
   ScLine  = mxGetPr(prhs[1]);	       /* data */
   color  = mxGetPr(prhs[2]);	       /* colormap */
   colorm = mxGetM(prhs[2]);
   colorn = mxGetN(prhs[2]);
   Ncol   = colorm;		       /* number of colors */
   printf("%d,%d %d,%d\n",Width, Height,colorm,colorn);
   R = color;G = R + colorm;B = G + colorm;

  /*   C H E C K   P A R A M E T E R S   */
  
  if (Width <= 0 || Width > 4096 || Height <= 0 || Height > 4096) {
     sprintf(info_str,
	     "WriteGif: incorrect image size: %d x %d", Width, Height);
     mexErrMsgTxt(info_str);
  }

   if (Ncol <= 0 || Ncol > 256) {
      sprintf(info_str,"WriteGif: wrong number of colors: %d", Ncol);
      mexErrMsgTxt(info_str);
   }

  /*   I N I T I A L I S A T I O N   */

   Nbyte  = 0;

  /*   F I N D   #   O F   B I T S   P E R    P I X E L   */

   BitsPixel = 1;  
   if (Ncol > 2)   BitsPixel = 2;  
   if (Ncol > 4)   BitsPixel = 3;  
   if (Ncol > 8)   BitsPixel = 4;  
   if (Ncol > 16)  BitsPixel = 5;  
   if (Ncol > 32)  BitsPixel = 6;  
   if (Ncol > 64)  BitsPixel = 7;  
   if (Ncol > 128) BitsPixel = 8;  

   ncol  = 1 << BitsPixel;
   IniCodeSize = BitsPixel;
   if (BitsPixel <= 1) IniCodeSize = 2;

  /*   W R I T E   H E A D E R  */

   fprintf(File, "GIF87a");	       /* Magic header */
   put_short((short)Width, File);                     /* screen size */
   put_short((short)Height, File);

   K  = 0x80;                            /* yes, there is a color map */
   K |= (8-1)<<4;                        /* OR in the color resolution */
   K |= (BitsPixel - 1);                 /* OR in the # of bits per pixel */

   fputc(K, File);
   fputc(0, File);		       /* background color */
   fputc(0, File);		       /* future expansion byte */

/* Write out the color map */
   for (i=0; i < Ncol; i++) {	       /* global colormap */
      fputc((int)(R[i]*maxVal), File);
      fputc((int)(G[i]*maxVal), File);
      fputc((int)(B[i]*maxVal), File);
   }
   for (; i<ncol; i++) {
      fputc(0, File);
      fputc(0, File);
      fputc(0, File);
   }

   fputc(',', File);		       /* image separator */
   put_short(0, File);                         /* left offset of image */
   put_short(0, File);                         /* top offset of image */
   put_short((short)Width, File);                     /* image size */
   put_short((short)Height, File); 
   fputc(0, File);		       /* no local colors, no interlace */
   fputc(IniCodeSize, File);	       /* initial code size */

  /*   L W Z   C O M P R E S S I O N   */

   CurCodeSize = ++IniCodeSize;
   CurMaxCode  = (1 << (IniCodeSize)) - 1;
   ClearCode   = (1 << (IniCodeSize - 1));
   EOFCode     = ClearCode + 1;
   FreeCode    = ClearCode + 2;
   output(ClearCode, File, glob);
   x     = 0;
   Code  = (int)(ScLine[x++]-0.5);
   while(x < (Width*Height)) {
      K     = (int)(ScLine[x++]-0.5);	       /* next symbol */
      CodeK = ((long) K << BITS) + Code;/* set full code */
      i     = (K << SHIFT) ^ Code;     /* xor hashing */
	 
      if (HashTab[i] == CodeK) {       /* full code found */
	 Code = CodeTab[i];
	 continue;
      }
      else if (HashTab[i] < 0 )	       /* empty slot */
      goto NOMATCH;
      
      disp  = HSIZE - i;	       /* secondary hash */
      if (i == 0) disp = 1;
	 
PROBE:
      if ((i -= disp) < 0)
      	i  += HSIZE;
	 
      if (HashTab[i] == CodeK) {       /* full code found */
	 Code = CodeTab[i];
	 continue;
      }
	 
      if (HashTab[i] > 0)	       /* try again */
      goto PROBE;
	 
NOMATCH:
      output(Code, File, glob);	       /* full code not found */
      Code = K;
	 
      if (FreeCode < THELIMIT) {
	 CodeTab[i] = FreeCode++;      /* code -> hashtable */
	 HashTab[i] = CodeK;
      } 
      else
      	output(ClearCode, File, glob);
   }
   /*   O U T P U T   T H E   R E S T  */

   output(Code, File, glob);
   output(EOFCode, File, glob);
   fputc(0, File);		       /* zero-length packet (EOF) */
   fputc(';', File);		       /* GIF file terminator */

   fclose(File);		       /* close the file */
}

/***************************************************************
 *                                                             *
 * Name: output                                 Date: 02.10.92 *
 *                                                             *
 * Function: outpt GIF code                                    *
 *                                                             *
 * Input: code - GIF code                                      *
 *                                                             *
 ***************************************************************/
static void output(int code, FILE *File, struct gl *glob)
{
  /*   O U T P U T   C O D E   */

   cur_accum &= masks[cur_bits];
   if (cur_bits > 0)
     cur_accum |= ((long)code << cur_bits);
   else
     cur_accum = code;
   cur_bits += CurCodeSize;
   while( cur_bits >= 8 ) {
      char_out( (unsigned int) (cur_accum & 0xFF) , File, glob);
      cur_accum >>= 8;
      cur_bits -= 8;
   }

  /*   R E S E T   */

   if (code == ClearCode ) {
      memset((char *) HashTab, -1, sizeof(HashTab));
      FreeCode = ClearCode + 2;
      CurCodeSize = IniCodeSize;
      CurMaxCode  = (1 << (IniCodeSize)) - 1;
   }

  /*   I N C R E A S E   C O D E   S I Z E   */

   if (FreeCode > CurMaxCode ) {
      CurCodeSize++;
      if ( CurCodeSize == BITS )
      	CurMaxCode = THELIMIT;
      else
      	CurMaxCode = (1 << (CurCodeSize)) - 1;
   }

  /*   E N D   O F   F I L E :  write the rest of the buffer  */

   if( code == EOFCode ) {
      while( cur_bits > 0 ) {
	 char_out( (unsigned int)(cur_accum & 0xff) ,File, glob);
	 cur_accum >>= 8;
	 cur_bits -= 8;
      }
     char_flush(File, glob);
  }
}

static void char_out(int c, FILE *File, struct gl *glob)
{
   accum[a_count++] = c;
   if (a_count >= 254) 
     char_flush(File, glob);
}

static void char_flush(FILE *File, struct gl *glob)
{
   if (a_count)  {
      fputc(a_count, File);
      fwrite(accum,1,a_count,File);
      a_count = 0;
   }
}

static void put_short(short word, FILE *File)
{
  fputc(word & 0xFF, File);
  fputc((word>>8) & 0xFF, File);
}
