#include "example.h"

/*
   This gives a very crude C example that contours a "real world"
   dataset using the CONREC algorthm. This is by no means a general
   solution and uses many specific features of the particular data
   being contoured.
*/

/* The dimension of the grid */
int NX = 0;
int NY = 0;

/* 
   The data will be scaled up by this amount to give some
   resolution for the contour lines.
*/
#define SCALE 5

/* Two dimensional array of data */
double **data;

/* Arrays for the x axis and y axis coordinates */
double *xaxis,*yaxis;

/* Array for the contour levels, 5 of them */
#define NCONTOUR 5
double contours[NCONTOUR];

/* Image on which the contours will be drawn, see bitmaplib.c */
BITMAP4 *image;

/* Debugging - count the number of line segments drawn */
int vectorsdrawn = 0;

int main(int argc,char **argv)
{
   int i,j,ii,jj;
   int n=0,lines = 0;
   double sum;
   double x,y,z;
   double xmin=1e32,xmax=-1e32;
   double ymin=1e32,ymax=-1e32;
   double zmin=1e32,zmax=-1e32;
   COLOUR colour;
   BITMAP4 col,grey = {128,128,128,0};
   FILE *fptr;

   /* 
      Open the data file and read the header.
      This is specific to the test dataset being used here,
      the header consists of the x and y dimensions of the data grid
      (which need not be fill)
   */
   if ((fptr = fopen("example.data","r")) == NULL) {
      fprintf(stderr,"Failed to open data file\n");
      exit(-1);
   }
   if (fscanf(fptr,"%d %d",&NX,&NY) != 2) {
      fprintf(stderr,"Expected to be able to read width and height of grid\n");
      exit(-1);
   }
   if (NX < 4 || NY < 4 || NX > 1000 || NY > 1000) {
      fprintf(stderr,"Got an unexpected grid size (%d x %d)\n",NX,NY);
      exit(-1);
   }

   /* 
      Malloc space for the two dimensional data 
      Initialise it to 0
   */
   if ((data = malloc(SCALE*NX*sizeof(double *))) == NULL) {
      fprintf(stderr,"Failed to malloc space for the data\n");
      exit(-1);
   }
   for (i=0;i<SCALE*NX;i++) {
      if ((data[i] = malloc(SCALE*NY*sizeof(double))) == NULL) {
         fprintf(stderr,"Failed to malloc space for the data\n");
         exit(-1);
      }
   }
   for (i=0;i<SCALE*NX;i++)
      for (j=0;j<SCALE*NY;j++)
         data[i][j] = 0;

   /*
      Read the data.
      This is VERY specific to the details of this example data.
   */
   while (fscanf(fptr,"%lf %lf %lf",&x,&y,&z) == 3) {

      /* Spread the data over a block, for future smoothing */
      for (ii=0;ii<SCALE;ii++)
         for (jj=0;jj<SCALE;jj++)
            data[SCALE*(int)x+ii][SCALE*(int)y+jj] = z;

      /* Check the bounds - debugging purposes */
      xmin = MIN(xmin,x);
      xmax = MAX(xmax,x);
      ymin = MIN(ymin,y);
      ymax = MAX(ymax,y);
      zmin = MIN(zmin,z);
      zmax = MAX(zmax,z);
      lines++;
   }
   fclose(fptr);
   fprintf(stderr,"Read %d lines\n",lines);
   fprintf(stderr,"X Range: %lf -> %lf\n",xmin,xmax);
   fprintf(stderr,"Y Range: %lf -> %lf\n",ymin,ymax);
   fprintf(stderr,"Z Range: %lf -> %lf\n",zmin,zmax);

   /* 
      Smooth the data
      This is just a simple 4x4 rectangular filter, it isn't
      actually necessary but makes the result sexier. 
   */
   for (i=0;i<SCALE*NX;i++) {
      for (j=0;j<NY*SCALE;j++) {
         n = 0;
         sum = 0;
         for (ii=-4;ii<=4;ii++) {
            for (jj=-4;jj<=4;jj++) {
               if (i + ii < 0 || i + ii >= SCALE*NX)
                  continue;
               if (j + jj < 0 || j + jj >= SCALE*NY)
                  continue;
               sum += data[i+ii][j+jj];
               n++;
            }
         }
         if (n <= 0) {
            fprintf(stderr,"No cells averaged, this shouldn't happen!\n");
            exit(-1);
         }
         data[i][j] = sum / n;
      }
   }

   /* 
      Set up the axis coordinates 
      If helps to do this with thought so the line segments drawn
      by CONREC are in the most convenient coordinate system.
   */
   if ((xaxis = malloc(NX*SCALE*sizeof(double))) == NULL) {
      fprintf(stderr,"Failed to malloc space for the xaxis data\n");
      exit(-1);
   }
   for (i=0;i<SCALE*NX;i++) 
      xaxis[i] = i;
   if ((yaxis = malloc(NY*SCALE*sizeof(double))) == NULL) {
      fprintf(stderr,"Failed to malloc space for the yaxis data\n");
      exit(-1);
   }
   for (i=0;i<SCALE*NY;i++) 
      yaxis[i] = i;

   /* Set up the contour levels */
   contours[0] = zmin + (zmax - zmin) / 8;
   contours[1] = zmin + (zmax - zmin) / 4;
   contours[2] = (zmax + zmin) / 2;
   contours[3] = zmax - (zmax - zmin) / 4;
   contours[4] = zmax - (zmax - zmin) / 8;

   /* 
      Create the image in memory and draw a colour ramped version 
      GetColour() just performs a colour ramp from blue to red.
   */
   if ((image = Create_Bitmap(NX*SCALE,NY*SCALE)) == NULL) {
      fprintf(stderr,"Malloc of bitmap failed\n");
      exit(-1);
   }
   Erase_Bitmap(image,NX*SCALE,NY*SCALE,grey); /* Not strictly necessary */
   for (j=0;j<SCALE*NY;j++) {
      for (i=0;i<SCALE*NX;i++) {
         colour = GetColour(data[i][j],zmin,zmax);
         col.r = colour.r * 255;
         col.g = colour.g * 255;
         col.b = colour.b * 255;
         Draw_Pixel(image,SCALE*NX,SCALE*NY,(double)i,(double)j,col);
      }
   }

   /* Finally do the contouring */
   CONREC(data,0,SCALE*NX-1,0,SCALE*NY-1,
      xaxis,yaxis,NCONTOUR,contours,drawline);
   fprintf(stderr,"Drew %d vectors\n",vectorsdrawn);

   /* 
      Write the image as a TGA file 
      See bitmaplib.c for more details, or write "image"
      in your own prefered format.
   */
   if ((fptr = fopen("image.tga","w")) == NULL) {
      fprintf(stderr,"Failed to open output image\n");
      exit(-1);
   }
   Write_Bitmap(fptr,image,SCALE*NX,SCALE*NY,12);
   fclose(fptr);

	exit(0);
}

void drawline(double x1,double y1,double x2,double y2,double z)
{
   BITMAP4 black = {0,0,0,0};

   if (x1 < 0 || x1 >= NX*SCALE || x2 < 0 || x2 > NX*SCALE)
      fprintf(stderr,"Shouldn't get here, x out of bounds: %g %g\n",x1,x2);
   if (y1 < 0 || y1 >= NY*SCALE || y2 < 0 || y2 > NY*SCALE)
      fprintf(stderr,"Shouldn't get here, y out of bounds: %g %g\n",y1,y2);
   Draw_Line(image,SCALE*NX,SCALE*NY,(int)x1,(int)y1,(int)x2,(int)y2,black);
   vectorsdrawn++;
}

/*
   Derivation from CONREC
   d               ! matrix of data to contour
   ilb,iub,jlb,jub ! index bounds of data matrix
   x               ! data matrix column coordinates
   y               ! data matrix row coordinates
   nc              ! number of contour levels
   z               ! contour levels in increasing order
*/
void CONREC(double **d,int ilb,int iub,int jlb,int jub,
   double *x,double *y,int nc,double *z,
   void (*ConrecLine)(double,double,double,double,double))
{
#define xsect(p1,p2) (h[p2]*xh[p1]-h[p1]*xh[p2])/(h[p2]-h[p1])
#define ysect(p1,p2) (h[p2]*yh[p1]-h[p1]*yh[p2])/(h[p2]-h[p1])

   int m1,m2,m3,case_value;
   double dmin,dmax,x1,x2,y1,y2;
   int i,j,k,m;
   double h[5];
   int sh[5];
   double xh[5],yh[5];
   int im[4] = {0,1,1,0},jm[4]={0,0,1,1};
   int castab[3][3][3] = {
     { {0,0,8},{0,2,5},{7,6,9} },
     { {0,3,4},{1,3,1},{4,3,0} },
     { {9,6,7},{5,2,0},{8,0,0} }
   };
   double temp1,temp2;

   for (j=(jub-1);j>=jlb;j--) {
      for (i=ilb;i<=iub-1;i++) {
         temp1 = MIN(d[i][j],d[i][j+1]);
         temp2 = MIN(d[i+1][j],d[i+1][j+1]);
         dmin  = MIN(temp1,temp2);
         temp1 = MAX(d[i][j],d[i][j+1]);
         temp2 = MAX(d[i+1][j],d[i+1][j+1]);
         dmax  = MAX(temp1,temp2);
         if (dmax < z[0] || dmin > z[nc-1])
            continue;
         for (k=0;k<nc;k++) {
            if (z[k] < dmin || z[k] > dmax)
               continue;
            for (m=4;m>=0;m--) {
               if (m > 0) {
                  h[m]  = d[i+im[m-1]][j+jm[m-1]]-z[k];
                  xh[m] = x[i+im[m-1]];
                  yh[m] = y[j+jm[m-1]];
               } else {
                  h[0]  = 0.25 * (h[1]+h[2]+h[3]+h[4]);
                  xh[0] = 0.50 * (x[i]+x[i+1]);
                  yh[0] = 0.50 * (y[j]+y[j+1]);
               }
               if (h[m] > 0.0)
                  sh[m] = 1;
               else if (h[m] < 0.0)
                  sh[m] = -1;
               else
                  sh[m] = 0;
            }

            /*
               Note: at this stage the relative heights of the corners and the
               centre are in the h array, and the corresponding coordinates are
               in the xh and yh arrays. The centre of the box is indexed by 0
               and the 4 corners by 1 to 4 as shown below.
               Each triangle is then indexed by the parameter m, and the 3
               vertices of each triangle are indexed by parameters m1,m2,and m3.
               It is assumed that the centre of the box is always vertex 2
               though this isimportant only when all 3 vertices lie exactly on
               the same contour level, in which case only the side of the box
               is drawn.
                  vertex 4 +-------------------+ vertex 3
                           | \               / |
                           |   \    m-3    /   |
                           |     \       /     |
                           |       \   /       |
                           |  m=2    X   m=2   |       the centre is vertex 0
                           |       /   \       |
                           |     /       \     |
                           |   /    m=1    \   |
                           | /               \ |
                  vertex 1 +-------------------+ vertex 2
            */
            /* Scan each triangle in the box */
            for (m=1;m<=4;m++) {
               m1 = m;
               m2 = 0;
               if (m != 4)
                  m3 = m + 1;
               else
                  m3 = 1;
               if ((case_value = castab[sh[m1]+1][sh[m2]+1][sh[m3]+1]) == 0)
                  continue;
               switch (case_value) {
               case 1: /* Line between vertices 1 and 2 */
                   x1 = xh[m1];
                   y1 = yh[m1];
                   x2 = xh[m2];
                   y2 = yh[m2];
                   break;
               case 2: /* Line between vertices 2 and 3 */
                   x1 = xh[m2];
                   y1 = yh[m2];
                   x2 = xh[m3];
                   y2 = yh[m3];
                   break;
               case 3: /* Line between vertices 3 and 1 */
                   x1 = xh[m3];
                   y1 = yh[m3];
                   x2 = xh[m1];
                   y2 = yh[m1];
                   break;
               case 4: /* Line between vertex 1 and side 2-3 */
                   x1 = xh[m1];
                   y1 = yh[m1];
                   x2 = xsect(m2,m3);
                   y2 = ysect(m2,m3);
                   break;
               case 5: /* Line between vertex 2 and side 3-1 */
                   x1 = xh[m2];
                   y1 = yh[m2];
                   x2 = xsect(m3,m1);
                   y2 = ysect(m3,m1);
                   break;
               case 6: /* Line between vertex 3 and side 1-2 */
                   x1 = xh[m1];
                   y1 = yh[m1];
                   x2 = xsect(m1,m2);
                   y2 = ysect(m1,m2);
                   break;
               case 7: /* Line between sides 1-2 and 2-3 */
                   x1 = xsect(m1,m2);
                   y1 = ysect(m1,m2);
                   x2 = xsect(m2,m3);
                   y2 = ysect(m2,m3);
                   break;
               case 8: /* Line between sides 2-3 and 3-1 */
                   x1 = xsect(m2,m3);
                   y1 = ysect(m2,m3);
                   x2 = xsect(m3,m1);
                   y2 = ysect(m3,m1);
                   break;
               case 9: /* Line between sides 3-1 and 1-2 */
                   x1 = xsect(m3,m1);
                   y1 = ysect(m3,m1);
                   x2 = xsect(m1,m2);
                   y2 = ysect(m1,m2);
                   break;
               default:
                   break;
               }

               /* Finally draw the line */
               ConrecLine(x1,y1,x2,y2,z[k]);
            } /* m */
         } /* k - contour */
      } /* i */
   } /* j */
}

/*
   Create a bitmap structure
*/
BITMAP4 *Create_Bitmap(int nx,int ny)
{
   return((BITMAP4 *)malloc(nx*ny*sizeof(BITMAP4)));
}

/* 
   Clear the bitmap to a particular colour
*/
void Erase_Bitmap(BITMAP4 *bm, int nx, int ny, BITMAP4 col)
{  
   int i,j; 
   long index; 
   
   for (i=0;i<nx;i++) {
      for (j=0;j<ny;j++) {
         index = j * nx + i;
         bm[index] = col;
      }
   }
}    

/*
   Turn on a pixel of a bitmap
*/
int Draw_Pixel(BITMAP4 *bm, int nx, int ny, int x, int y, BITMAP4 col)
{
   long index;

   if (x < 0 || y < 0 || x >= nx || y >= ny)
      return(FALSE);
   index = y * nx + x;
   bm[index] = col;
   return(TRUE);
}

/*
   Draw a line from (x1,y1) to (x2,y2)
   Use colour col
*/
void Draw_Line(BITMAP4 *bm,int nx,int ny,int x1,int y1,int x2,int y2,BITMAP4 col)
{  
   int i,j; 
   long index; 
   double mu,dx,dy,dxy;
   
   dx = x2 - x1;
   dy = y2 - y1;
   dxy = sqrt(dx*dx + dy*dy);
   if (dxy <= 0) {
      Draw_Pixel(bm,nx,ny,x1,y1,col);
      return;
   }
   for (mu=0;mu<=2*dxy;mu++) {
      i = (int)(x1 + 0.5 * mu * dx / dxy);
      j = (int)(y1 + 0.5 * mu * dy / dxy);
      if (i < 0 || j < 0 || i >= nx || j >= ny)
         continue; 
      index = j * nx + i;
      bm[index] = col;
   } 
} 

/*
   Write a bitmap to a file
   The format is as follows
     1 == tga
    11 == tga with alpha
    12 == compressed tga
    13 == compressed tga with alpha
     2 == ppm
     3 == rgb
     4 == raw grey scale
     5 == tiff
     6 == EPS colour (Encapsulated PostScript)
     7 == EPS black and white
     8 == raw
     9 == BMP
   A negative format indicates a vertical flip
*/
void Write_Bitmap(FILE *fptr,BITMAP4 *bm,int nx,int ny,int format)
{
   int i,j,offset;
   long index,rowindex;
   int linelength = 0,size;
   char buffer[1024];

   /* Write the header */
   switch (ABS(format)) {
   case 1:
   case 11:
   case 12:
   case 13:
      putc(0,fptr);  /* Length of ID */
      putc(0,fptr);  /* No colour map */
      if (ABS(format) == 12 || ABS(format) == 13)
         putc(10,fptr); /* compressed RGB */
      else
         putc(2,fptr); /* uncompressed RGB  */
      putc(0,fptr); /* Index of colour map entry */
      putc(0,fptr);
      putc(0,fptr); /* Colour map length */
      putc(0,fptr);
      putc(0,fptr); /* Colour map size */
      putc(0,fptr); /* X origin */
      putc(0,fptr);
      putc(0,fptr); /* Y origin */
      putc(0,fptr);
      putc((nx & 0x00ff),fptr); /* X width */
      putc((nx & 0xff00) / 256,fptr);
      putc((ny & 0x00ff),fptr); /* Y width */
      putc((ny & 0xff00) / 256,fptr);
      if (ABS(format) == 11 || ABS(format) == 13) {
         putc(32,fptr);                      /* 32 bit bitmap     */
         putc(0x08,fptr);
      } else {
         putc(24,fptr);                       /* 24 bit bitmap       */
         putc(0x00,fptr);
      }
      break;
   case 2:
      fprintf(fptr,"P6\n# bitmaplib (Paul Bourke)\n%d %d\n255\n",nx,ny);
      break;
   case 3:
      putc(0x01,fptr);
      putc(0xda,fptr);
      putc(0x00,fptr);
      putc(0x01,fptr);
      putc(0x00,fptr);
      putc(0x03,fptr);
      putc((nx & 0xFF00) / 256,fptr);
      putc((nx & 0x00FF),fptr);
      putc((ny & 0xFF00) / 256,fptr);
      putc((ny & 0x00FF),fptr);
      BM_WriteHexString(fptr,"000300000000000000ff00000000");
      fprintf(fptr,"WriteBitmap, pdb");
      putc(0x00,fptr);
      putc(0x00,fptr);
      putc(0x00,fptr);
      putc(0x00,fptr);
      putc(0x00,fptr);
      putc(0x00,fptr);
      putc(0x00,fptr);
      putc(0x00,fptr);
      break;
   case 4:
      break;
   case 5:
      BM_WriteHexString(fptr,"4d4d002a");   /* Little endian & TIFF identifier */
      offset = nx * ny * 3 + 8;
      BM_WriteLongInt(fptr,buffer,offset);
      break;
   case 6:
      fprintf(fptr,"%%!PS-Adobe-3.0 EPSF-3.0\n");
      fprintf(fptr,"%%%%Creator: Created from bitmaplib by Paul Bourke\n");
      fprintf(fptr,"%%%%BoundingBox: %d %d %d %d\n",0,0,nx,ny);
      fprintf(fptr,"%%%%LanguageLevel: 2\n");
      fprintf(fptr,"%%%%Pages: 1\n");
      fprintf(fptr,"%%%%DocumentData: Clean7Bit\n");
      fprintf(fptr,"%d %d scale\n",nx,ny);
      fprintf(fptr,"%d %d 8 [%d 0 0 -%d 0 %d]\n",nx,ny,nx,ny,ny);
      fprintf(fptr,"{currentfile 3 %d mul string readhexstring pop} bind\n",nx);
      fprintf(fptr,"false 3 colorimage\n");
      break;
   case 7:
      fprintf(fptr,"%%!PS-Adobe-3.0 EPSF-3.0\n");
      fprintf(fptr,"%%%%Creator: Created from bitmaplib by Paul Bourke\n");
      fprintf(fptr,"%%%%BoundingBox: %d %d %d %d\n",0,0,nx,ny);
      fprintf(fptr,"%%%%LanguageLevel: 2\n");
      fprintf(fptr,"%%%%Pages: 1\n");
      fprintf(fptr,"%%%%DocumentData: Clean7Bit\n");
      fprintf(fptr,"%d %d scale\n",nx,ny);
      fprintf(fptr,"%d %d 8 [%d 0 0 -%d 0 %d]\n",nx,ny,nx,ny,ny);
      fprintf(fptr,"{currentfile %d string readhexstring pop} bind\n",nx);
      fprintf(fptr,"false 1 colorimage\n");
      break;
   case 8:
      break;
   case 9:
      /* Header 10 bytes */
      putc('B',fptr);
      putc('M',fptr);
      size = nx * ny * 3 + 14 + 40;
      putc((size) % 256,fptr);
      putc((size / 256) % 256,fptr);
      putc((size / 65536) % 256,fptr);
      putc((size / 16777216),fptr);
      putc(0,fptr); putc(0,fptr);
      putc(0,fptr); putc(0,fptr);
      /* Offset to image data */
      putc(14+40,fptr); putc(0,fptr); putc(0,fptr); putc(0,fptr);
      /* Information header 40 bytes */
      putc(0x28,fptr); putc(0,fptr); putc(0,fptr); putc(0,fptr);
      putc((nx) % 256,fptr);
      putc((nx / 256) % 256,fptr);
      putc((nx / 65536) % 256,fptr);
      putc((nx / 16777216),fptr);
      putc((ny) % 256,fptr);
      putc((ny / 256) % 256,fptr);
      putc((ny / 65536) % 256,fptr);
      putc((ny / 16777216),fptr);
      putc(1,fptr); putc(0,fptr); /* One plane */
      putc(24,fptr); putc(0,fptr); /* 24 bits */
      /* Compression type == 0 */
      putc(0,fptr); putc(0,fptr); putc(0,fptr); putc(0,fptr);
      size = nx * ny * 3;
      putc((size) % 256,fptr);
      putc((size / 256) % 256,fptr);
      putc((size / 65536) % 256,fptr);
      putc((size / 16777216),fptr);
      putc(1,fptr); putc(0,fptr); putc(0,fptr); putc(0,fptr);
      putc(1,fptr); putc(0,fptr); putc(0,fptr); putc(0,fptr);
      putc(0,fptr); putc(0,fptr); putc(0,fptr); putc(0,fptr); /* No palette */
      putc(0,fptr); putc(0,fptr); putc(0,fptr); putc(0,fptr);
      break;
   }


   // Write the binary data
   for (j=0;j<ny;j++) {
      if (format > 0)
         rowindex = j * nx;
      else
         rowindex = (ny - 1 - j) * nx;
      switch (ABS(format)) {
      case 12:
         WriteTGACompressedRow(fptr,&(bm[rowindex]),nx,3);
         break;
      case 13:
         WriteTGACompressedRow(fptr,&(bm[rowindex]),nx,4);
         break;
      }
      for (i=0;i<nx;i++) {
         if (format > 0)
            index = rowindex + i;
         else
            index = rowindex + i;
         switch (ABS(format)) {
         case 1:
         case 11:
         case 9:
            putc(bm[index].b,fptr);
            putc(bm[index].g,fptr);
            putc(bm[index].r,fptr);
            if (ABS(format) == 11)
               putc(bm[index].a,fptr);
            break;
         case 2:
         case 3:
         case 5:
         case 8:
            putc(bm[index].r,fptr);
            putc(bm[index].g,fptr);
            putc(bm[index].b,fptr);
            break;
         case 4:
            putc((bm[index].r+bm[index].g+bm[index].b)/3,fptr);
            break;
         case 6:
            fprintf(fptr,"%02x%02x%02x",bm[index].r,bm[index].g,bm[index].b);
            linelength += 6;
            if (linelength >= 72 || linelength >= nx) {
               fprintf(fptr,"\n");
               linelength = 0;
            }
            break;
         case 7:
            fprintf(fptr,"%02x",(bm[index].r+bm[index].g+bm[index].b)/3);
            linelength += 2;
            if (linelength >= 72 || linelength >= nx) {
               fprintf(fptr,"\n");
               linelength = 0;
            }
            break;
         }
      }
   }

   /* Write the footer */
   switch (ABS(format)) {
   case 1:
   case 11:
   case 12:
   case 13:
   case 2:
   case 3:
   case 4:
      break;
   case 5:
      putc(0x00,fptr); /* The number of directory entries (14) */
      putc(0x0e,fptr);

      /* Width tag, short int */
      BM_WriteHexString(fptr,"0100000300000001");
      putc((nx & 0xff00) / 256,fptr);      /* Image width */
      putc((nx & 0x00ff),fptr);
      putc(0x00,fptr);
      putc(0x00,fptr);

      /* Height tag, short int */
      BM_WriteHexString(fptr,"0101000300000001");
      putc((ny & 0xff00) / 256,fptr);    /* Image height */
      putc((ny & 0x00ff),fptr);
      putc(0x00,fptr);
      putc(0x00,fptr);

      /* bits per sample tag, short int */
      BM_WriteHexString(fptr,"0102000300000003");
      offset = nx * ny * 3 + 182;
      BM_WriteLongInt(fptr,buffer,offset);

      /* Compression flag, short int */
      BM_WriteHexString(fptr,"010300030000000100010000");

      /* Photometric interpolation tag, short int */
      BM_WriteHexString(fptr,"010600030000000100020000");

      /* Strip offset tag, long int */
      BM_WriteHexString(fptr,"011100040000000100000008");

      /* Orientation flag, short int */
      BM_WriteHexString(fptr,"011200030000000100010000");

      /* Sample per pixel tag, short int */
      BM_WriteHexString(fptr,"011500030000000100030000");

      /* Rows per strip tag, short int */
      BM_WriteHexString(fptr,"0116000300000001");
      putc((ny & 0xff00) / 256,fptr);
      putc((ny & 0x00ff),fptr);
      putc(0x00,fptr);
      putc(0x00,fptr);
     /* Strip byte count flag, long int */
      BM_WriteHexString(fptr,"0117000400000001");
      offset = nx * ny * 3;
      BM_WriteLongInt(fptr,buffer,offset);

      /* Minimum sample value flag, short int */
      BM_WriteHexString(fptr,"0118000300000003");
      offset = nx * ny * 3 + 188;
      BM_WriteLongInt(fptr,buffer,offset);

      /* Maximum sample value tag, short int */
      BM_WriteHexString(fptr,"0119000300000003");
      offset = nx * ny * 3 + 194;
      BM_WriteLongInt(fptr,buffer,offset);

      /* Planar configuration tag, short int */
      BM_WriteHexString(fptr,"011c00030000000100010000");

      /* Sample format tag, short int */
      BM_WriteHexString(fptr,"0153000300000003");
      offset = nx * ny * 3 + 200;
      BM_WriteLongInt(fptr,buffer,offset);

      /* End of the directory entry */
      BM_WriteHexString(fptr,"00000000");

      /* Bits for each colour channel */
      BM_WriteHexString(fptr,"000800080008");

      /* Minimum value for each component */
      BM_WriteHexString(fptr,"000000000000");

      /* Maximum value per channel */
      BM_WriteHexString(fptr,"00ff00ff00ff");

      /* Samples per pixel for each channel */
      BM_WriteHexString(fptr,"000100010001");

      break;
   case 6:
   case 7:
      fprintf(fptr,"\n%%%%EOF\n");
      break;
   case 8:
   case 9:
      break;
   }
}

void BM_WriteLongInt(FILE *fptr,char *s,long n)
{     
   int i;
      
   s[0] = (n & 0xff000000) / 16777216;
   s[1] = (n & 0x00ff0000) / 65536;
   s[2] = (n & 0x0000ff00) / 256;
   s[3] = (n & 0x000000ff);
      
   for (i=0;i<4;i++)
      putc(s[i],fptr);
}     
      
void BM_WriteHexString(FILE *fptr,char *s)
{
   unsigned int i,c;
   char hex[3];
      
   for (i=0;i<strlen(s);i+=2) {
      hex[0] = s[i];
      hex[1] = s[i+1];
      hex[2] = '\0';
      sscanf(hex,"%X",&c);
      putc(c,fptr);
   }  
}     

/*    
   Write a compressed TGA row 
   Depth is either 3 or 4
*/    
void WriteTGACompressedRow(FILE *fptr,BITMAP4 *bm,int width,int depth)
{        
   int i;
   int counter = 1;
   int pixelstart = 0;
   int packettype = 0;
   int readytowrite = FALSE;
   BITMAP4 currentpixel,nextpixel = {0,0,0,0};
            
   currentpixel = bm[0];
   for (;;) {
      if (pixelstart+counter >= width)
         readytowrite = TRUE;
      else
         nextpixel = bm[pixelstart+counter];
            
      if (!readytowrite) {
         if (Same_BitmapPixel(currentpixel,nextpixel)) {
            if (packettype == 0) {
               counter++;
               if (counter >= 128 || (pixelstart + counter) >= width)
                  readytowrite = TRUE;
            } else {
               counter--;
               readytowrite = TRUE;
            }
         } else {
            if (packettype == 1 || counter <= 1) {
               packettype = 1;
               currentpixel = nextpixel;
               counter++;
               if (counter >= 128 || (pixelstart + counter) >= width)
                  readytowrite = TRUE;
            } else {
               readytowrite = TRUE;
            }
         }
      }

      if (readytowrite) {
         if (pixelstart + counter > width)
            counter = width - pixelstart;
         if (packettype == 0) {
            putc(((counter-1) | 0x80),fptr);
            putc(currentpixel.b,fptr);
            putc(currentpixel.g,fptr);
            putc(currentpixel.r,fptr);
            if (depth == 4)
               putc(currentpixel.a,fptr);
            currentpixel = nextpixel;
         } else {
            putc(counter-1,fptr);
            for (i=0;i<counter;i++) {
               putc(bm[pixelstart+i].b,fptr);
               putc(bm[pixelstart+i].g,fptr);
               putc(bm[pixelstart+i].r,fptr);
               if (depth == 4)
                  putc(bm[pixelstart+i].a,fptr);
            }
         }
         if ((pixelstart = pixelstart + counter) >= width)
            break; /* From for (;;) */
         readytowrite = FALSE;
         packettype = 0;
         counter = 1;
      }
   }
}

/*
   Compare two pixels
*/
int Same_BitmapPixel(BITMAP4 p1,BITMAP4 p2)
{
   if (p1.r != p2.r) return(FALSE);
   if (p1.g != p2.g) return(FALSE);
   if (p1.b != p2.b) return(FALSE);
   if (p1.a != p2.a) return(FALSE);
   return(TRUE);
}

COLOUR GetColour(double v,double vmin,double vmax)
{   
   int iv;
   double dv,vmid;
   COLOUR c = {1.0,1.0,1.0};

   if (vmax < vmin) {
      dv = vmin;
      vmin = vmax;
      vmax = dv;
   }
   if (vmax - vmin < 0.000001) {
      vmin -= 1;
      vmax += 1;
   }
   if (v < vmin)
      v = vmin;
   if (v > vmax)
      v = vmax;
   dv = vmax - vmin; 
      
   if (v < (vmin + 0.25 * dv)) {
      c.r = 0;
      c.g = 4 * (v - vmin) / dv;
      c.b = 1;
   } else if (v < (vmin + 0.5 * dv)) {
      c.r = 0;
      c.g = 1;
      c.b = 1 + 4 * (vmin + 0.25 * dv - v) / dv;
   } else if (v < (vmin + 0.75 * dv)) {
      c.r = 4 * (v - vmin - 0.5 * dv) / dv;
      c.g = 1;
      c.b = 0;
   } else {
      c.r = 1;
      c.g = 1 + 4 * (vmin + 0.75 * dv - v) / dv;
      c.b = 0;
   }

	return(c);
}