//  ____      __     
// |  __'____|  |__  
// |  _|  _  |  _  | 
// '_' '___._'_____' 
// Florian & Alain's Blender 
// 2009 januar - october 

// !!! APPAREMENT DES BUGS DANS LES DIVISIONS 


#define ALAIN

#include <stdio.h>  // file functions, etc... 
#include <stdlib.h> // exit(), malloc() 
#include <string.h> // string handling 
#include <math.h>   // mathematics
#include <time.h>

#define FALSE 0
#define TRUE 1
#define pi 3.1415926

/**** GLOBAL VARIABLES *****************************************************/
unsigned char ColorPalette[256][3];
FILE *filep,*fopen(const char *,const char *);

struct IMAGE {
  unsigned short *Adress;
  unsigned int Length;
};
typedef struct IMAGE IMAGE;
IMAGE Img0,Img1,Img2;

/* IMPORTANT: The number of bytes in one line MUST always be a multiple of 4,
   so we must sometimes append 0s to align. */
unsigned char *BmpAdress;
unsigned int BmpLength;
unsigned char BmpHeader[54]={
  'B','M',  /*  0: always "BM" */
  0,0,0,0,  /*  2: size of file */
  0,0,0,0,  /*  6:always set to 0, reserved (?) */
  54,0,0,0, /* 10: offset from start of file to bitmap data */
  40,0,0,0, /* 14: size of Infos (54-14) */
  0,0,0,0,  /* 18: width of image */
  0,0,0,0,  /* 22: height of image */
  1,0,      /* 26: nb of planes (?), always 1 */
  24,0,     /* 28: nb of bitplanes in {1;4;8;24} */
  0,0,0,0,  /* 30: compression mode (0=none) */
  0,0,0,0,  /* 34: uncompressed size of bitmap data (0 if no compression is ok) */
  0,0,0,0,  /* 38: horizontal pixels per meter */
  0,0,0,0,  /* 42: vertical pixels per meter */
  0,0,0,0,  /* 46: number of colors used in the bitmap (may be 0) */
  0,0,0,0   /* 50: number of colors important (0=>all ) */
};

/**** ROUTINES *************************************************************/


//**** rnd ****
// Returns an unsigned number in [rndmin;rndmax]. 
// Generator from D.E.Knuth's Art of Computer Programming, vol.2, p185. 
// Period of pseudorandom sequence is 2147483646. 
// Using modulo is mathematically incorrect but who cares? 
#define rndM 2147483647 // 2^31-1, a Mersenne prime 
#define rndA 48271      // this does well in the spectral test 
#define rndQ 44488      // rndM/rndA 
#define rndR 3399       // rndM%rndA, it's important that rndR<rndQ 
long rndX;              // should be initialised to 0<rndX<rndM 

unsigned int rnd(rndmin,rndmax) unsigned int rndmin,rndmax; {
  rndX=rndA*(rndX%rndQ)-rndR*(long)(rndX/rndQ); // X=AX mod M 
  if(rndX<0) {rndX+=rndM; }
  return( rndmin+((rndX)%(rndmax-rndmin+1)));
}          

/**** division with rounding ****/
/* Please note that 255:2 = 127.5 must be rounded to 127 in case we make 127*2 */
int divide(dividend,divisor) int dividend,divisor; {
  int quotient;
  if(divisor<0) { divisor=-divisor; dividend=-dividend; }
  quotient=dividend/divisor;
  if((dividend%divisor)>(divisor/2)) { quotient++; }
  return(quotient);
}

/**** power2 ****/
int power2(n) int n; { if(n>0) { return(2*power2(n-1)); } else { return(1); } }

void ErrorInFile(name,text) char *name,*text; {printf("Error: %s %s\n",name,text); exit(1);}

void CheckUsage(c,n,s1,s2) int c,n; char *s1; char *s2; {
  if(c<n) { printf("Usage: fab %s %s\n",s1,s2); exit(1); }
  printf("%s: ",s1);
}

/**** Word/LongWord Handling for BMP files ****/
void WriteLongWord(value,pos,adress) unsigned int value,pos; unsigned char *adress; {
  adress[pos]=value & 255; value=value>>8;
  adress[pos+1]=value & 255; value=value>>8;
  adress[pos+2]=value & 255; value=value>>8;
  adress[pos+3]=value & 255; value=value>>8;
}

void WriteWord(value,pos,adress) unsigned int value,pos; unsigned char *adress; {
  adress[pos]=value & 255; value=value>>8;
  adress[pos+1]=value & 255; value=value>>8;
}

int ReadLongWord(pos,adress) unsigned int pos; unsigned char *adress; {
  return(((adress[pos+3]*256+adress[pos+2])*256+adress[pos+1])*256+adress[pos]);
}

int ReadWord(pos,adress) unsigned int pos; unsigned char *adress; {
  return(adress[pos+1]*256+adress[pos]);
}

/**** Allocate Memory for w*h grey image ****/
IMAGE AllocMemForGreyImage(width,height) int width,height; {
  IMAGE Img;
  if((width<0) || (width>65535)) { ErrorInFile("image","too large"); }
  if((height<0) || (height>65535)) { ErrorInFile("image","too high"); }
  Img.Length=(2+width*height)*sizeof(unsigned short); /* size in bytes */
  Img.Adress=(unsigned short *) malloc(Img.Length);
  Img.Adress[0]=width;
  Img.Adress[1]=height;
  return(Img);
}

/**** Allocate Memory and Load Grey Image File ****/
IMAGE AllocMemAndLoadGreyImage(name) char *name; {
  IMAGE Img;
  filep=fopen(name,"rb"); /* rb = read binary */
  if(filep==0) { ErrorInFile(name,"not found"); }
  fseek(filep,0,2); /* set the file pointer at the end of file */
  Img.Length=(int) ftell(filep); /* returns the length of the file in bytes */
  fseek(filep,0,0); /* set the file pointer at the start of file */
  Img.Adress=(unsigned short *) malloc(Img.Length); /* allocate memory */
  fread(Img.Adress,1,Img.Length,filep); /* load file, 1=size of a byte */
  fclose(filep);
  return(Img);
}

/**** Save Image ****/
void SaveImage(name,Img) char *name; IMAGE Img; {
  unsigned int FileLength;
  filep=fopen(name,"wb"); /* write binary */
  FileLength=fwrite(Img.Adress,1,Img.Length,filep);
  fclose(filep);
  if(FileLength!=Img.Length) { ErrorInFile(name,"not saved"); }
}

/**** searchmincolinImg0 ****/
unsigned short searchmincolinImg0() {
  int x,y;
  unsigned short width0,height0,mincol,currentcol;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  mincol=Img0.Adress[2];
  for(y=0;y<height0;y++) {
    for(x=0;x<width0;x++) {
      currentcol=Img0.Adress[2+y*width0+x];
      if(currentcol<mincol) { mincol=currentcol; }
    }
  }
  return (mincol);  
}

/**** searchmaxcolinImg0 ****/
unsigned short searchmaxcolinImg0() {
  int x,y;
  unsigned short width0,height0,maxcol,currentcol;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  maxcol=Img0.Adress[2];
  for(y=0;y<height0;y++) {
    for(x=0;x<width0;x++) {
      currentcol=Img0.Adress[2+y*width0+x];
      if(currentcol>maxcol) { maxcol=currentcol; }
    }
  }
  return (maxcol);  
}

/******************************************************************************/
/************************** FAB FUNCTIONS *************************************/
/******************************************************************************/

/**** bmp2rgb ****/
void bmp2rgb(name) char *name; {
  unsigned int width,realbmpwidth,height;
  int pos,x,y;
  /* load bmp file */
  filep=fopen(name,"rb"); /* rb = read binary */
  if(filep==0) { ErrorInFile(name,"not found"); }
  fseek(filep,0,2); /* set the file pointer at the end of file */
  BmpLength=(int) ftell(filep); /* returns the length of the file in bytes */
  fseek(filep,0,0); /* set the file pointer at the start of file */
  BmpAdress=(unsigned char *) malloc(BmpLength); /* allocate memory */
  fread(BmpAdress,1,BmpLength,filep); /* load file, 1=size of a byte */
  fclose(filep);
  /* get width & height and check image validity */
  if(BmpLength<=54) { ErrorInFile(name,"must be more than 54 bytes"); }
  if((BmpAdress[0]!='B') || (BmpAdress[1]!='M')) { ErrorInFile(name,"isn't a .bmp file"); }
  if(BmpLength!=ReadLongWord(2,BmpAdress)) { ErrorInFile(name,"file corrupted"); }
  width=ReadLongWord(18,BmpAdress);
  height=ReadLongWord(22,BmpAdress);
  if((width<0) || (width>65535)) { ErrorInFile(name,"too large"); }
  if((height<0) || (height>65535)) { ErrorInFile(name,"too high"); }
  printf("%d x %d",width,height);
  realbmpwidth=3*width;
  while((realbmpwidth & 3)!=0) { realbmpwidth++; }
  pos=ReadLongWord(10,BmpAdress);
  if(pos<54) { ErrorInFile(name,"file corrupted"); }
  Img0=AllocMemForGreyImage(width,height); /* allocate memory for rgb channels  */
  Img1=AllocMemForGreyImage(width,height);
  Img2=AllocMemForGreyImage(width,height);
  /* convert image */
  switch(ReadWord(28,BmpAdress)) {
    case 24:
      /* convert 24 bpp to 3 channels red.g green.g blue.g in Img0 Img1 Img2 */
      if(pos+3*width*height>BmpLength){ ErrorInFile(name,"file corrupted"); }
      for(y=0;y<height;y++) {
        for(x=0;x<width;x++) {
          Img0.Adress[2+y*width+x]=(unsigned short) BmpAdress[pos+y*realbmpwidth+x*3+2]; /* red */
          Img1.Adress[2+y*width+x]=(unsigned short) BmpAdress[pos+y*realbmpwidth+x*3+1]; /* green */
          Img2.Adress[2+y*width+x]=(unsigned short) BmpAdress[pos+y*realbmpwidth+x*3+0]; /* blue */
        }
      }
      break;
    case 8:
      /* convert 8 bpp paletted to 3 channels red.g green.g blue.g in Img0 Img1 Img2 */
      ErrorInFile(name,"is a 8 bpp file (not yet supported)");
      break;
    case 4: ErrorInFile(name,"is a 4 bpp file (not yet supported)"); break;
    case 1: ErrorInFile(name,"is a 1 bpp file (not yet supported)"); break;
    default: ErrorInFile(name,"file corrupted"); break;
  }
}

/**** rgb2bmp ****/
void rgb2bmp(name) char *name; {
  unsigned int width,realbmpwidth,height,width0,height0,width1,height1,width2,height2;
  unsigned int FileLength;
  int pos,x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  width1=Img1.Adress[0];
  height1=Img1.Adress[1];
  width2=Img2.Adress[0];
  height2=Img2.Adress[1];
  /* get common size */
  width=width0;
  if(width1<width) { width=width1; }
  if(width2<width) { width=width2; }
  height=height0;
  if(height1<height) { height=height1; }
  if(height2<height) { height=height2; }
  printf("%d x %d",width,height);
  realbmpwidth=width*3;
  while((realbmpwidth & 3)!=0) { realbmpwidth++; }
  BmpLength=54+realbmpwidth*height;
  BmpAdress=(unsigned char *) malloc(BmpLength); /* allocate memory */
  for(pos=0;pos<54;pos++) { BmpAdress[pos]=BmpHeader[pos]; } /* copy header */
  WriteLongWord(54+realbmpwidth*height,2,BmpAdress); /* size of file */
  WriteLongWord(width,18,BmpAdress); /* width */
  WriteLongWord(height,22,BmpAdress); /* height */
  for(x=0;x<realbmpwidth*height;x++) { BmpAdress[pos+x]=0; }
  for(y=0;y<height;y++) {
    for(x=0;x<width;x++) {
      BmpAdress[pos+x*3+y*realbmpwidth+0]=Img2.Adress[2+x+y*width2]; /* blue */
      BmpAdress[pos+x*3+y*realbmpwidth+1]=Img1.Adress[2+x+y*width1]; /* green */
      BmpAdress[pos+x*3+y*realbmpwidth+2]=Img0.Adress[2+x+y*width0]; /* red */
    }
  }
  filep=fopen(name,"wb"); /* write binary */
  FileLength=fwrite(BmpAdress,1,BmpLength,filep);
  fclose(filep);
  if(FileLength!=BmpLength) { ErrorInFile(name,"not saved"); }
}

/**** rgb2grey ****/
void rgb2grey() {
  unsigned int width,height,width0,height0,width1,height1,width2,height2;
  int x,y,t;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  width1=Img1.Adress[0];
  height1=Img1.Adress[1];
  width2=Img2.Adress[0];
  height2=Img2.Adress[1];
  /* get common size */
  width=width0;
  if(width1<width) { width=width1; }
  if(width2<width) { width=width2; }
  height=height0;
  if(height1<height) { height=height1; }
  if(height2<height) { height=height2; }
  for(y=0;y<height;y++) {
    for(x=0;x<width;x++) {
      t=Img0.Adress[2+x+y*width0]+Img1.Adress[2+x+y*width1]+Img2.Adress[2+x+y*width2]; /* red+green+blue */
      t=divide(t,3);
      Img0.Adress[2+x+y*width]=t;
    }
  }
}

/**** negative ****/
void negative() {
  unsigned int c,width0,height0;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  for(c=0;c<width0*height0;c++) {
    Img0.Adress[2+c]=65535-Img0.Adress[2+c];
  }
}

/**** and ****/
void and() {
  unsigned int width0,height0,width1,height1,width2,height2;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  width1=Img1.Adress[0];
  height1=Img1.Adress[1];
  if (width0<width1){width2=width0;}
  else {width2=width1;}
  if (height0<height1){height2=height0;}
  else {height2=height1;}
  Img2=AllocMemForGreyImage(width2,height2);
  for (y=0;y<height2;y++) {
    for (x=0;x<width2;x++) {
      Img2.Adress[2+width2*y+x]=Img0.Adress[2+width0*y+x]&Img1.Adress[2+width1*y+x];
    }
  }
}
/**** andI ****/
void andI(value) unsigned int value; {
  unsigned int width0,height0;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img2=AllocMemForGreyImage(width0,height0);
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      Img2.Adress[2+width0*y+x]=Img0.Adress[2+width0*y+x]&value;
    }
  }
}
/**** or ****/
void or() {
  unsigned int width0,height0,width1,height1,width2,height2;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  width1=Img1.Adress[0];
  height1=Img1.Adress[1];
  if (width0<width1){width2=width0;}
  else {width2=width1;}
  if (height0<height1){height2=height0;}
  else {height2=height1;}
  Img2=AllocMemForGreyImage(width2,height2);
  for (y=0;y<height2;y++) {
    for (x=0;x<width2;x++) {
      Img2.Adress[2+width2*y+x]=Img0.Adress[2+width0*y+x]|Img1.Adress[2+width1*y+x];
    }
  }
}
/**** orI ****/
void orI(value) unsigned int value; {
  unsigned int width0,height0;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img2=AllocMemForGreyImage(width0,height0);
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      Img2.Adress[2+width0*y+x]=Img0.Adress[2+width0*y+x]|value;
    }
  }
}
/**** xor ****/
void xor() {
  unsigned int width0,height0,width1,height1,width2,height2;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  width1=Img1.Adress[0];
  height1=Img1.Adress[1];
  if (width0<width1){width2=width0;}
  else {width2=width1;}
  if (height0<height1){height2=height0;}
  else {height2=height1;}
  Img2=AllocMemForGreyImage(width2,height2);
  for (y=0;y<height2;y++) {
    for (x=0;x<width2;x++) {
      Img2.Adress[2+width2*y+x]=Img0.Adress[2+width0*y+x]^Img1.Adress[2+width1*y+x];
    }
  }
}
/**** xorI ****/
void xorI(value) unsigned int value; {
  unsigned int width0,height0;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img2=AllocMemForGreyImage(width0,height0);
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      Img2.Adress[2+width0*y+x]=Img0.Adress[2+width0*y+x]^value;
    }
  }
}  
  
/**** addborders ****/
void addborders(left,right,bottom,top,color) int left,right,bottom,top,color; {
  unsigned int width0,height0,width1,height1;
  int x,y;
  printf("%d %d %d %d",left,right,bottom,top);
  if(left<0) { printf("Error: left should be positive"); exit(0); }
  if(right<0) { printf("Error: right should be positive"); exit(0); }
  if(bottom<0) { printf("Error: bottom should be positive"); exit(0); }
  if(top<0) { printf("Error: top should be positive"); exit(0); }
  if((color<0) || (color>65535)) { printf("Error: color should be in [0;65535]"); exit(0); }
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  width1=left+width0+right;
  height1=bottom+height0+top;
  printf(" -> %d %d",width1,height1);
  Img1=AllocMemForGreyImage(width1,height1);
  for (x=2;x<2+width1*height1;x++) {Img1.Adress[x]=color;}
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      Img1.Adress[2+(bottom+y)*width1+left+x]=Img0.Adress[2+y*width0+x];
    }
  }
}
/**** removeborders ****/
void removeborders(left,right,bottom,top) int left,right,bottom,top; {
  unsigned int width0,height0,width1,height1;
  int x,y;
  printf("%d %d %d %d",left,right,bottom,top);
  if(left<0) { printf("Error: left should be positive"); exit(0); }
  if(right<0) { printf("Error: right should be positive"); exit(0); }
  if(bottom<0) { printf("Error: bottom should be positive"); exit(0); }
  if(top<0) { printf("Error: top should be positive"); exit(0); }
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  if(left+right>width0) { printf("Error: width do not exist"); exit(0); }
  if(top+bottom>height0) { printf("Error: height do not exist"); exit(0); }
  width1=width0-(left+right);
  height1=height0-(bottom+top);
  printf(" -> %d %d",width1,height1);
  Img1=AllocMemForGreyImage(width1,height1);
  for (y=0;y<height1;y++) {
    for (x=0;x<width1;x++) {
      Img1.Adress[2+y*width1+x]=Img0.Adress[2+(y+top)*width0+x+left];
    }
  }
}
/**** add ****/
void add() {
  unsigned int width0,height0,width1,height1,width2,height2;
  int x,y,t;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  width1=Img1.Adress[0];
  height1=Img1.Adress[1];
  if (width0<width1){width2=width0;} else {width2=width1;}
  if (height0<height1){height2=height0;} else {height2=height1;}
  Img2=AllocMemForGreyImage(width2,height2);
  for (y=0;y<height2;y++) {
    for (x=0;x<width2;x++) {
      t=Img0.Adress[2+width0*y+x]+Img1.Adress[2+width1*y+x];
      if (t>65535) {
        Img2.Adress[2+width2*y+x]=65535;
      } else {
        Img2.Adress[2+width2*y+x]=t;
      }
    }
  }
}

/**** addI ****/
void addI(value) unsigned int value; {
  unsigned int width0,height0;
  int x,y,t;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img2=AllocMemForGreyImage(width0,height0);
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      t=Img0.Adress[2+width0*y+x]+value;
      if(t>65535) { t=65535; }
      Img2.Adress[2+width0*y+x]=t;
    }
  }
} 

/**** sub ****/
void sub() {
  unsigned int width0,height0,width1,height1,width2,height2;
  int x,y,t;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  width1=Img1.Adress[0];
  height1=Img1.Adress[1];
  if (width0<width1){width2=width0;} else {width2=width1;}
  if (height0<height1){height2=height0;} else {height2=height1;}
  Img2=AllocMemForGreyImage(width2,height2);
  for (y=0;y<height2;y++) {
    for (x=0;x<width2;x++) {
      t=Img0.Adress[2+width0*y+x]-Img1.Adress[2+width1*y+x];
      if (t<0) {
        Img2.Adress[2+width2*y+x]=0;
      } else {
        Img2.Adress[2+width2*y+x]=t;
      }
    }
  }
} 

/**** subI ****/
void subI(value) unsigned int value; {
  unsigned int width0,height0;
  int x,y,t;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img2=AllocMemForGreyImage(width0,height0);
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      t=Img0.Adress[2+width0*y+x]-value;
      if(t<0) { t=0; }
      Img2.Adress[2+width0*y+x]=t;
    }
  }
} 

/**** rsbI ****/
void rsbI(value) unsigned int value; {
  unsigned int width0,height0;
  int x,y,t;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img2=AllocMemForGreyImage(width0,height0);
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      t=value-Img0.Adress[2+width0*y+x];
      if(t<0) { t=0; }
      Img2.Adress[2+width0*y+x]=t;
    }
  }
} 

/**** mul ****/
void mul() {
  unsigned int width0,height0,width1,height1,width2,height2;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  width1=Img1.Adress[0];
  height1=Img1.Adress[1];
  if (width0<width1){width2=width0;} else {width2=width1;}
  if (height0<height1){height2=height0;} else {height2=height1;}
  Img2=AllocMemForGreyImage(width2,height2);
  for (y=0;y<height2;y++) {
    for (x=0;x<width2;x++) {
      Img2.Adress[2+width2*y+x]=Img0.Adress[2+width0*y+x]*Img1.Adress[2+width1*y+x];
    }
  }
} 

/**** mulI ****/
void mulI(factor) unsigned int factor; {
  unsigned int width0,height0;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img2=AllocMemForGreyImage(width0,height0);
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      Img2.Adress[2+width0*y+x]=Img0.Adress[2+width0*y+x]*factor;
    }
  }
} 

/**** div ****/
void division() {
  unsigned int width0,height0,width1,height1,width2,height2;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  width1=Img1.Adress[0];
  height1=Img1.Adress[1];
  if (width0<width1){width2=width0;} else {width2=width1;}
  if (height0<height1){height2=height0;} else {height2=height1;}
  Img2=AllocMemForGreyImage(width2,height2);
  for (y=0;y<height2;y++) {
    for (x=0;x<width2;x++) {
      Img2.Adress[2+width2*y+x]=divide(Img0.Adress[2+width0*y+x],Img1.Adress[2+width1*y+x]);
    }
  }
} 

/**** divI ****/
void divI(divisor) unsigned int divisor; {
  unsigned int width0,height0;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img2=AllocMemForGreyImage(width0,height0);
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      Img2.Adress[2+width0*y+x]=divide(Img0.Adress[2+width0*y+x],divisor);
    }
  }
} 

/**** filter ****/
/* IN  Img0 */
/* OUT Img1 */
void filter(a,b,c,d,e,f,g,h,i,n,offset) int a,b,c,d,e,f,g,h,i,n,offset; {
  unsigned int width,height;
  int x0,x1,x2,y,y0,y1,y2,pixel;
  unsigned short *Adress;
  width=Img0.Adress[0];
  height=Img0.Adress[1];
  Adress=Img0.Adress+2;
  Img1=AllocMemForGreyImage(width,height);
  for(y=0;y<height;y++) {
    y0=((y+height-1)%height)*width; /* if y-1 then modulo doesn't act the intended way */
    y1=y*width;
    y2=((y+1)%height)*width;
    for(x1=0;x1<width;x1++) {
      x0=(x1+width-1)%width;
      x2=(x1+1)%width;
      pixel=a*Adress[y0+x0]+b*Adress[y0+x1]+c*Adress[y0+x2];
      pixel+=d*Adress[y1+x0]+e*Adress[y1+x1]+f*Adress[y1+x2];
      pixel+=g*Adress[y2+x0]+h*Adress[y2+x1]+i*Adress[y2+x2];
      pixel=divide(pixel,n)+offset;
      if(pixel<0) { pixel=0; }
      if(pixel>65535) { pixel=65535; }
      Img1.Adress[2+y1+x1]=pixel;
    }
  }
}

/**** mknoise ****/
void mknoise(min,max) short min,max; {
  unsigned int width0,height0;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      Img0.Adress[2+y*width0+x]=rnd(min,max);
    }
  }  
}

/***** hmirror *****/
void hmirror() {
  unsigned int width0,height0;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img2=AllocMemForGreyImage(width0,height0);
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      Img2.Adress[2+(height0-1-y)*width0+x]=Img0.Adress[2+y*width0+x];
    }
  }
}

/***** vmirror *****/
void vmirror() {
  unsigned int width0,height0;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img2=AllocMemForGreyImage(width0,height0);
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      Img2.Adress[2+y*width0+width0-1-x]=Img0.Adress[2+y*width0+x];
    }
  }
}

/***** rotate90 *****/
void rotate90() {
  unsigned int width0,height0;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img2=AllocMemForGreyImage(height0,width0);
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      Img2.Adress[2+x*height0+height0-1-y]=Img0.Adress[2+y*width0+x];
    }
  } 
}

/***** rotate180 *****/
void rotate180() {
  unsigned int width0,height0;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img2=AllocMemForGreyImage(width0,height0);
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      Img2.Adress[2+(height0-1-y)*width0+width0-1-x]=Img0.Adress[2+y*width0+x];
    }
  } 
}

/***** rotate 270 *****/
void rotate270() {
  unsigned int width0,height0;
  int x,y;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img2=AllocMemForGreyImage(height0,width0);
  for (y=0;y<height0;y++) {
    for (x=0;x<width0;x++) {
      Img2.Adress[2+(width0-1-x)*height0+y]=Img0.Adress[2+y*width0+x];
    }
  } 
}

/**** edge ****/
void edge(size,back_color,edge_color) unsigned int size,back_color,edge_color; {
  unsigned int width0,height0,dist,dist2,size2;
  int x,y,yy,found;
  size2=size*size;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img1=AllocMemForGreyImage(height0,width0); /* hdist2 */
  /* Compute and store squares of horizontal distances */
  for(y=0;y<height0;y++) {
    /* from left to right */    
    dist=size;
    dist2=size2;
    for(x=0;x<width0;x++) {
      if(Img0.Adress[2+y*width0+x]==back_color) {
        if(dist<255) { dist2+=2*dist+1; } /* (x+1)^2=x^2+2x+1 */
        dist+=1; /* we won't increase dist2 if dist>=254 */
      } else {
        dist2=0;
        dist=0;
      }
      Img1.Adress[y*width0+x]=dist2;
    }
    /* from right to left */    
    dist=size;
    dist2=size2;
    for(x=width0-1;x>=0;x--) {
      if(Img0.Adress[2+y*width0+x]==back_color) {
        if(dist<255) { dist2+=2*dist+1; }
        dist+=1;
      } else {
        dist2=0;
        dist=0;
      }
      if(dist2<Img1.Adress[y*width0+x]) { Img1.Adress[y*width0+x]=dist2; }
    }
  }
  /* Compute dist2 using previous table and Pythagora's theorem */
  for(y=0;y<height0;y++) {
    for(x=0;x<width0;x++) {
      if(Img0.Adress[2+y*width0+x]==back_color) {
        found=FALSE;
        dist=0;
        dist2=0;
        while((dist<size) && (found==FALSE)) {
          yy=y+dist;
          if((yy<height0) && (dist2+Img1.Adress[yy*width0+x]<size2)) { found=TRUE; }
          yy=y-dist;
          if((yy>=0) && (dist2+Img1.Adress[yy*width0+x]<size2)) { found=TRUE; }
          if(dist<255) { dist2+=2*dist+1; }
          dist++;
        }
        if(found==TRUE) { Img0.Adress[2+y*width0+x]=edge_color; }
      }
    }
  }
}

/***** explode *****/
void explode(NbIt) unsigned int NbIt; {
  unsigned int width0,height0,x,y,c,t;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  for(c=0;c<NbIt;c++) {
    x=rnd(0,width0-2);
    y=rnd(0,height0-2);
    t=Img0.Adress[2+y*width0+x];
    if(rnd(0,1)==0) {
      Img0.Adress[2+y*width0+x]=Img0.Adress[2+y*width0+x+1];
      Img0.Adress[2+y*width0+x+1]=t;
    } else {
      Img0.Adress[2+y*width0+x]=Img0.Adress[2+(y+1)*width0+x];
      Img0.Adress[2+(y+1)*width0+x]=t;
    }
  }
}

/**** pal2rgb ****/
void pal2rgb(name) char *name; {
  unsigned int width0,height0,NbColors;
  int x,y,PaletteLength;
  unsigned short t;
  unsigned char *Palette;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  Img1=AllocMemForGreyImage(width0,height0);
  Img2=AllocMemForGreyImage(width0,height0);
  filep=fopen(name,"rb"); /* rb = read binary */
  if(filep==0) { ErrorInFile(name,"not found"); }
  fseek(filep,0,2); /* set the file pointer at the end of file */
  PaletteLength=(int) ftell(filep); /* returns the length of the file in bytes */
  fseek(filep,0,0); /* set the file pointer at the start of file */
  Palette=(unsigned char *) malloc(PaletteLength); /* allocate memory */
  fread(Palette,1,PaletteLength,filep); /* load file, 1=size of a byte */
  fclose(filep);
  NbColors=PaletteLength/3;
  for(y=0;y<height0;y++) {
    for(x=0;x<width0;x++) {
      t=Img0.Adress[2+y*width0+x];
      if(t>=NbColors) { t=NbColors-1; }
      Img0.Adress[2+y*width0+x]=Palette[t*3];
      Img1.Adress[2+y*width0+x]=Palette[t*3+1];
      Img2.Adress[2+y*width0+x]=Palette[t*3+2];
    }
  }
}

/**** setmincol ****/
void setmincol(mincol) unsigned short mincol; {
  int x,y;
  unsigned short width0,height0,currentmincol;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  currentmincol=searchmincolinImg0();
  printf("oldmincol=%d",currentmincol);
  for(y=0;y<height0;y++) {
    for(x=0;x<width0;x++) {
       Img0.Adress[2+y*width0+x]=Img0.Adress[2+y*width0+x]+mincol-currentmincol;
    }
  }
}

/**** setmaxcol ****/
void setmaxcol(maxcol) unsigned short maxcol; {
  int x,y;
  unsigned short width0,height0,currentmaxcol;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  currentmaxcol=searchmaxcolinImg0();
  printf("oldmaxcol=%d",currentmaxcol);
  for(y=0;y<height0;y++) {
    for(x=0;x<width0;x++) {
       Img0.Adress[2+y*width0+x]=Img0.Adress[2+y*width0+x]+maxcol-currentmaxcol;
    }
  }
}

/***** rectangle *****/
void rectangle(top,bottom,left,right,color) unsigned int top,bottom,left,right,color; {
  int x,y;
  unsigned short width0,height0;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  printf("%d %d %d %d",left,right,bottom,top);
  if(left<0) { printf("Error: left should be positive"); exit(0); }
  if(right<0) { printf("Error: right should be positive"); exit(0); }
  if(bottom<0) { printf("Error: bottom should be positive"); exit(0); }
  if(top<0) { printf("Error: top should be positive"); exit(0); }
  if((color<0) || (color>65535)) { printf("Error: color should be in [0;65535]"); exit(0); }
  if (top+bottom>height0) {printf("Error: height too small"); exit(0); }
  if (left+right>width0) {printf("Error: width too small"); exit(0); }
  for (y=bottom;y<height0-top;y++) {
    for (x=left;x<width0-right;x++) {
      Img0.Adress[2+y*width0+x]=color;
    }
  }
}

/***** rectangle2 *****/
void rectangle2(bottom,left,length,height,color) unsigned int bottom,left,length,height,color; {
  int x,y;
  unsigned short width0,height0;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  printf("%d %d %d %d",bottom,left,length,height);
  if(left<0) { printf("Error: left should be positive"); exit(0); }
  if(length<0) { printf("Error: length should be positive"); exit(0); }
  if(height<0) { printf("Error: height should be positive"); exit(0); }
  if(bottom<0) { printf("Error: top should be positive"); exit(0); }
  if((color<0) || (color>65535)) { printf("Error: color should be in [0;65535]"); exit(0); }
  if (bottom+height>height0) {printf("Error: height too small"); exit(0); }
  if (left+length>width0) {printf("Error: width too small"); exit(0); }
  for (y=bottom;y<height+bottom;y++) {
    for (x=left;x<left+length;x++) {
      Img0.Adress[2+y*width0+x]=color;
    }
  }
}

/**** landscape ****/
void landscape(nbiter,fractalfactor,absorption) int nbiter,fractalfactor,absorption; {
  int x,y,xx,yy,t,i;
  unsigned short w,h;
  unsigned short a,b,c,d;
  unsigned short *temp;
  w=Img0.Adress[0];
  h=Img0.Adress[1];
  t=power2(nbiter);
  Img1=AllocMemForGreyImage(w*t,h*t);
  Img2=AllocMemForGreyImage(w*t,h*t);
  
  for(x=2;x<2+w*h;x++) {Img1.Adress[x]=Img0.Adress[x];}
  for(i=0;i<nbiter;i++) {
    fractalfactor=power2(nbiter-i);
    
    for(y=0;y<h;y++) {
      yy=y+1; if(yy==h) { yy=0; }
      /*printf("yy=%d\n",yy);*/
      for(x=0;x<w;x++) {
        xx=x+1; if(xx==w) { xx=0; }
        a=Img1.Adress[2+y*w+x];
        b=Img1.Adress[2+y*w+xx];
        c=Img1.Adress[2+yy*w+x];
        d=Img1.Adress[2+yy*w+xx];
        
        /* printf("%d %d\n",2+yy*w+xx,d); */
        
        Img2.Adress[2+2*(y*2*w+x)]=a;
        t=(a+b)/2+rnd(-fractalfactor,fractalfactor);
        if(t<0) { t=0;} if(t>255) { t=255; }
        Img2.Adress[2+2*(y*2*w+x)+1]=t;
        t=(a+c)/2+rnd(-fractalfactor,fractalfactor);
        if(t<0) { t=0;} if(t>255) { t=255; }
        Img2.Adress[2+2*(y*2*w+x)+2*w]=t;
        t=(a+b+c+d)/4+rnd(-fractalfactor,fractalfactor);
        if(t<0) { t=0;} if(t>255) { t=255; }
        Img2.Adress[2+2*(y*2*w+x)+2*w+1]=t;
        
        /*printf("(%d+%d)/2=%d\n",a,b,Img2.Adress[2+2*(y*2*w+x)+1]);
        printf("(%d+%d)/2=%d\n",a,c,Img2.Adress[2+2*(y*2*w+x)+2*w]);
        printf("(%d+%d+%d+%d)/2=%d\n",a,b,c,d,Img2.Adress[2+2*(y*2*w+x)+2*w+1]);*/
        
        
      }
    }
    /* fractalfactor=fractalfactor*100/absorption;*/
    
    temp=Img2.Adress;
    Img2.Adress=Img1.Adress;
    Img1.Adress=temp;
    w=2*w;
    h=2*h;
  }
}

/**** between ****/
void between(colmin,colmax) unsigned short colmin,colmax; {
  int x,y;
  unsigned short width0,height0,col;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  for(y=0;y<height0;y++) {
    for(x=0;x<width0;x++) {
      col=Img0.Adress[2+y*width0+x];
      if(col<colmin) { Img0.Adress[2+y*width0+x]=colmin; }
      if(col>colmax) { Img0.Adress[2+y*width0+x]=colmax; }
    }
  }
}

/**** range ****/
void range(colmin,colmax) unsigned short colmin,colmax; {
  int x,y;
  unsigned short width0,height0,col;
  unsigned short currentmincol,currentmaxcol;
  int delta;
  currentmincol=searchmincolinImg0();
  currentmaxcol=searchmaxcolinImg0();
  delta=currentmaxcol-currentmincol;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  for(y=0;y<height0;y++) {
    for(x=0;x<width0;x++) {
      col=Img0.Adress[2+y*width0+x];
      Img0.Adress[2+y*width0+x]=colmin+divide((colmax-colmin)*(col-currentmincol),delta);
    }
  }
}

/**** disc ****/
/* returns nb of pixels in disc */
int disc(cx,cy,radius,col) unsigned short cx,cy,radius,col; {
  int x,y,dx,dy,dist2,dist2max;
  int NbOfPixels=0;
  unsigned short width0,height0;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  if(radius!=0) {
    // print first line
    //!!!SHOULD BE if((cy>=0) && (cy<height0)) { BUT CY>=0 ALWAYS TRUE
    if(cy<height0) {
      for(x=cx-radius;x<=cx+radius;x++) {
        if((x>=0) && (x<width0)) { Img0.Adress[2+cy*width0+x]=col; }
      }
    }
    NbOfPixels+=2*radius+1;
    // now process two lines at a time
    dx=radius;
    dy=1;
    dist2=radius*radius;
    dist2max=dist2+radius;
    dist2+=1;
    while(dy<=radius) {
      while(dist2>dist2max) {
        dist2+=-2*dx+1; // dist2=(dx-1)^2+dy^2=old_dist2-2*dx+1
        dx-=1;
      }
      // draw two lines
      y=cy+dy;
      if((y>=0) && (y<height0)) {
        y=y*width0;
        for(x=cx-dx;x<=cx+dx;x++) {
          if((x>=0) && (x<width0)) { Img0.Adress[2+y+x]=col; }
        }
      }
      y=cy-dy;
      if((y>=0) && (y<height0)) {
        y=y*width0;
        for(x=cx-dx;x<=cx+dx;x++) {
          if((x>=0) && (x<width0)) { Img0.Adress[2+y+x]=col; }
        }
      }
      NbOfPixels+=2*(2*dx+1);
      dist2+=2*dy+1; // dist2=dx^2+(dy+1)^2=old_dist2+2*dy+1
      dy+=1;
    }
  }
  return(NbOfPixels);
}

/**** fillwithOs ****/
unsigned short FWC_width0,FWC_height0;
int FWC_cx,FWC_cy,FWC_radius;

// returns -1 if circle contains one (or more) non empty pixel
int FWC_CountEmptyPixelsInCircle() {
  int x,y1,y2,dx,dy,dist2,dist2max;
  int NbPixelsInDisc=0;
  // first line
  y1=FWC_cy*FWC_width0;
  for(x=FWC_cx-FWC_radius;x<=FWC_cx+FWC_radius;x++) {
    if(Img0.Adress[2+y1+x]!=0) { return(-1); }
  }
  NbPixelsInDisc+=2*FWC_radius+1;
  // now process two lines at a time
  dx=FWC_radius;
  dy=1;
  dist2=FWC_radius*FWC_radius;
  dist2max=dist2+FWC_radius;
  dist2+=1;
  while(dy<=FWC_radius) {
    while(dist2>dist2max) {
      dist2+=-2*dx+1; // dist2=(dx-1)^2+dy^2=old_dist2-2*dx+1
      dx-=1;
    }
    // check both hlines
    y1=(FWC_cy+dy)*FWC_width0;
    y2=(FWC_cy-dy)*FWC_width0;
    for(x=FWC_cx-FWC_radius;x<=FWC_cx+FWC_radius;x++) {
      if(Img0.Adress[2+y1+x]!=0) { return(-1); }
      if(Img0.Adress[2+y2+x]!=0) { return(-1); }
    }
    NbPixelsInDisc+=2*(2*dx+1);
    dist2+=2*dy+1; // dist2=dx^2+(dy+1)^2=old_dist2+2*dy+1
    dy+=1;
  }
  return(NbPixelsInDisc);
}

void fillwithOs(OutRMin,OutRMax,InRMin,InRMax,Border,ClearColor,CircleColor)
  unsigned short OutRMin,OutRMax,InRMin,InRMax,Border,ClearColor,CircleColor; {
    
  unsigned int c,NbFruitlessTries;
  int NbEmptyPixels;
  int x,y;
  FWC_width0=Img0.Adress[0];
  FWC_height0=Img0.Adress[1];
  // Count the number of empty pixels
  NbEmptyPixels=0;
  for(y=0;y<FWC_height0;y++) {
    for(x=0;x<FWC_width0;x++) {
      if(Img0.Adress[2+y*FWC_width0+x]==0) { NbEmptyPixels++; }
    }
  }
  
  // As long as we can, try to add a circle in black parts
  NbFruitlessTries=0;
  while((NbFruitlessTries<1000) && (NbEmptyPixels>0)) {
    // Choose the center amongst empty pixels
    c=rnd(1,NbEmptyPixels);
    FWC_cx=-1; FWC_cy=0;
    do {
      do {
        FWC_cx++;
        if(FWC_cx==FWC_width0) {
          FWC_cx=0;
          FWC_cy+=1;
        }
      } while(Img0.Adress[2+FWC_cy*FWC_width0+FWC_cx]!=0);
      c-=1;
    } while(c>0);
    // Here FWC_cx,FWC_cy is the center of the circle
    
    // Set up radius as rmax at first, clip to borders
    FWC_radius=OutRMax+Border;
    if(FWC_cx-FWC_radius<0) { FWC_radius=FWC_cx; }
    if(FWC_cx+FWC_radius>FWC_width0-1) { FWC_radius=FWC_width0-1-FWC_cx; }
    if(FWC_cy-FWC_radius<0) { FWC_radius=FWC_cy; }
    if(FWC_cy+FWC_radius>FWC_height0-1) { FWC_radius=FWC_height0-1-FWC_cy; }
    
    // See if (non clipped) circle contains only empty pixels, else decrease radius
    c=FWC_CountEmptyPixelsInCircle();
    while((c==-1) && (FWC_radius>=OutRMin)) {
      FWC_radius-=1;
      c=FWC_CountEmptyPixelsInCircle();
    }
    if(FWC_radius-Border>=OutRMin) {
      NbEmptyPixels-=c; // remove pixels from all circles
      c=disc(FWC_cx,FWC_cy,FWC_radius,ClearColor); // clear
      c-=disc(FWC_cx,FWC_cy,FWC_radius-Border,CircleColor); // draw circle
      if(FWC_radius-Border-1>=InRMin) {
        x=InRMax;
        if(FWC_radius-Border-1<x) { x=FWC_radius-Border-1; }
        FWC_radius=rnd(InRMin,x);
        c+=disc(FWC_cx,FWC_cy,FWC_radius,ClearColor); // clear inside
      }
      if(ClearColor==0) { NbEmptyPixels+=c; } // pixels==0 can be accessed again
      NbFruitlessTries=0;
    } else {
      NbFruitlessTries+=1;
    }
  }
}

/**** floyd ****/
//error propagated to other pixels
//    . 7
//  3 5 1
void floyd(NbColors) unsigned short NbColors; {
  int x,y;
  unsigned short width0,height0,col,newcol;
  unsigned short currentmincol,currentmaxcol;
  float amplitude,intensity,newerror,nbcols;
  
  float errorcurrentline[4096];
  float errornextline[4096];
  
  nbcols=(float) (NbColors);
  currentmincol=searchmincolinImg0();
  currentmaxcol=searchmaxcolinImg0();
  amplitude=(float) (currentmaxcol-currentmincol);
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  for(x=0;x<width0;x++) { errornextline[x]=0.0; }
  for(y=0;y<height0;y++) {
    for(x=0;x<width0;x++) {
      errorcurrentline[x]=errornextline[x];
      errornextline[x]=0.0;
    }
    for(x=0;x<width0;x++) {
      col=Img0.Adress[2+y*width0+x]-currentmincol;
      intensity=((float) col)/amplitude;
      intensity+=errorcurrentline[x];
      if(intensity>1.0) { intensity=1.0; }
      if(intensity<0.0) { intensity=0.0; }
      newcol=(unsigned short) (intensity*(nbcols-1));
//!!!      newerror=intensity- ((float) newcol)/nbcols;
//!!!      newerror=intensity-((float) newcol)/(nbcols-1.0);
      newerror=intensity-(0.5+(float) newcol)/nbcols;
      Img0.Adress[2+y*width0+x]=newcol;
      // propagate error
      if(x!=width0-1) {
        errorcurrentline[x+1]+=7*newerror/16;
        errornextline[x+1]+=1*newerror/16;
      }
      errornextline[x]+=5*newerror/16;
      if(x!=0) { errornextline[x-1]+=3*newerror/16; }
    }
  }
}

/**** subsample ****/
void subsample(xratio,yratio) short xratio,yratio; {
  unsigned int width0,height0,width1,height1;
  int x,y,xsrc,ysrc,xsrc2,ysrc2,nbpixelsinrectangle,sumpixels;
  nbpixelsinrectangle=xratio*yratio;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  if(width0%xratio!=0) { printf("Error: %d is not a multiple of %d",width0,xratio); exit(0); }
  if(height0%yratio!=0) { printf("Error: %d is not a multiple of %d",height0,yratio); exit(0); }
  width1=width0/xratio;
  height1=height0/yratio;
  Img1=AllocMemForGreyImage(width1,height1);
  for (y=0;y<height1;y++) {
    ysrc=y*yratio;
    for (x=0;x<width1;x++) {
      xsrc=x*xratio;
      sumpixels=0;
      for(ysrc2=ysrc;ysrc2<ysrc+yratio;ysrc2++) {
        for(xsrc2=xsrc;xsrc2<xsrc+xratio;xsrc2++) {
          sumpixels+=Img0.Adress[2+ysrc2*width0+xsrc2];
        }
      }
      Img1.Adress[2+y*width1+x]=sumpixels/nbpixelsinrectangle;
    }
  }
}

/**** mkstrippedtube ****/
void mkstrippedtube(diameter,nbsteps,nbstripes,col0,col1) short diameter,nbsteps,nbstripes; {
  unsigned int width0,height0;
  int x,y,t;
  double startangle,angle;
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  printf("%d x %d",width0,height0);
  for(y=0;y<nbsteps;y++) {
    startangle=2*pi*((double) y)/((double) nbstripes*nbsteps);
    for(x=0;x<diameter;x++) {
      angle=acos( ((double) (2*x-diameter))/((double) diameter));
      t=(short) (pi+(angle-startangle)*((double) nbstripes/pi)); //+pi to avoid negative value
      if(t%2==0) { t=col0; } else { t=col1; }
      Img0.Adress[2+y*width0+x]=t;
    }
  }
}


/**** grey2bitplanes ****/
void grey2bitplanes(nbbitplanes,name)  short nbbitplanes; char *name; {
  unsigned int FileLength;
  unsigned int width0,height0,nbwords;
  unsigned short int bitplane[nbbitplanes];
  unsigned char *File;
  FILE *filep;
  int x,y,word,pixel,t,FilePos;
  
  width0=Img0.Adress[0];
  height0=Img0.Adress[1];
  if(width0%16!=0) { printf("Error: width=%d is not a multiple of 16",width0); exit(0); }
  nbwords=width0/16;
  FileLength=2*nbbitplanes*nbwords*height0;
  File=(unsigned char *) malloc(FileLength);
  
  FilePos=0;
  for(y=0;y<height0;y++) {
    for(word=0;word<nbwords;word++) {
      for(t=0;t<nbbitplanes;t++) { bitplane[t]=0; }
      for(x=word*16;x<word*16+16;x++) {
        pixel=Img0.Adress[2+y*width0+x];
        for(t=0;t<nbbitplanes;t++) {
          bitplane[t]=2*bitplane[t]+(pixel & 1);
          pixel=pixel>>1;
        }
      }
      for(t=0;t<nbbitplanes;t++) {
        //Big endian vs small
        File[FilePos]=(unsigned char) (bitplane[t]>>8);
        FilePos++;
        File[FilePos]=(unsigned char) (bitplane[t] & 255);
        FilePos++;
      }
    }
  }
  
  filep=fopen(name,"wb"); /* write binary */
  t=fwrite(File,1,FileLength,filep);
  fclose(filep);
  if(FileLength!=t) { ErrorInFile(name,"not saved"); }
}
  

///**** Save Image ****/
//void SaveImage(name,Img) char *name; IMAGE Img; {
//  unsigned int FileLength;
//  filep=fopen(name,"wb"); /* write binary */
//  FileLength=fwrite(Img.Adress,1,Img.Length,filep);
//  fclose(filep);
//  if(FileLength!=Img.Length) { ErrorInFile(name,"not saved"); }
//}
    
    
/**** ReadNumber ****/
int ReadNumber(s) char *s; {
  int n,t,sign;
  unsigned char c;
  n=0;
  sign=1;
  t=0;
  c=s[t];
  if(c=='-') { sign=-1; t++; c=s[t]; }
  while(c!=0) {
    if((c<'0') || (c>'9')) {
      printf("Error: %s should be a number",s);
      exit(1);
    }
    n=n*10+c-'0';
    t++;
    c=s[t];
  }
  return(n*sign);
}

/**** MAIN *****************************************************************/
int main(argc,argv) int argc; char *argv[]; {
  setbuf(stdout, NULL); /* no buffering for output! */
  
  /* rndX=1; */ /* initialise random seed */
#ifdef ALAIN
  rndX=rawclock(); /* initialise random seed */
#else  
  rndX=time(0); /* initialise random seed */
#endif    

  if(argc<=2) { printf("Usage: fab function parameter ... parameter\nSee AllFunctions.txt"); exit(1); }
  
  if(strcmp(argv[1],"bmp2rgb")==0) {
    CheckUsage(argc,6,"bmp2rgb","in.bmp red.g green.g blue.g");
    bmp2rgb(argv[2]);
    SaveImage(argv[3],Img0);
    SaveImage(argv[4],Img1);
    SaveImage(argv[5],Img2);
    
  } else if (strcmp(argv[1],"rgb2bmp")==0) {
    CheckUsage(argc,6,"rgb2bmp","red.g green.g blue.g out.bmp");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    Img1=AllocMemAndLoadGreyImage(argv[3]);
    Img2=AllocMemAndLoadGreyImage(argv[4]);
    rgb2bmp(argv[5]);
    
  } else if (strcmp(argv[1],"rgb2grey")==0) {
    CheckUsage(argc,6,"rgb2grey","red.g green.g blue.g grey.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    Img1=AllocMemAndLoadGreyImage(argv[3]);
    Img2=AllocMemAndLoadGreyImage(argv[4]);
    rgb2grey();
    SaveImage(argv[5],Img0);
   
  } else if (strcmp(argv[1],"negative")==0) {
    CheckUsage(argc,4,"negative","in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    negative();
    SaveImage(argv[3],Img0);

  } else if (strcmp(argv[1],"filter")==0) {
    CheckUsage(argc,15,"filter","a b c d e f g h i divisor offset in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[13]);
    filter(ReadNumber(argv[2]),ReadNumber(argv[3]),ReadNumber(argv[4]),ReadNumber(argv[5]),ReadNumber(argv[6]),ReadNumber(argv[7]),ReadNumber(argv[8]),ReadNumber(argv[9]),ReadNumber(argv[10]),ReadNumber(argv[11]),ReadNumber(argv[12]));
    SaveImage(argv[14],Img1);
    
  } else if (strcmp(argv[1],"and")==0) {
    CheckUsage(argc,5,"and","in1.g in2.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    Img1=AllocMemAndLoadGreyImage(argv[3]);
    and();
    SaveImage(argv[4],Img2);
    
  } else if (strcmp(argv[1],"andI")==0) {
    CheckUsage(argc,5,"andI","value in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[3]);
    andI(ReadNumber(argv[2]));
    SaveImage(argv[4],Img2);
    
  } else if (strcmp(argv[1],"or")==0) {
    CheckUsage(argc,5,"or","in1.g in2.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    Img1=AllocMemAndLoadGreyImage(argv[3]);
    or();
    SaveImage(argv[4],Img2);
    
  } else if (strcmp(argv[1],"orI")==0) {
    CheckUsage(argc,5,"orI","value in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[3]);
    orI(ReadNumber(argv[2]));
    SaveImage(argv[4],Img2);
    
  } else if (strcmp(argv[1],"xor")==0) {
    CheckUsage(argc,5,"xor","in1.g in2.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    Img1=AllocMemAndLoadGreyImage(argv[3]);
    xor();
    SaveImage(argv[4],Img2);
    
  } else if (strcmp(argv[1],"xorI")==0) {
    CheckUsage(argc,5,"xorI","value in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[3]);
    xorI(ReadNumber(argv[2]));
    SaveImage(argv[4],Img2);
      
  } else if (strcmp(argv[1],"addborders")==0) {
    CheckUsage(argc,9,"addborders","left right bottom top color in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[7]);
    addborders(ReadNumber(argv[2]),ReadNumber(argv[3]),ReadNumber(argv[4]),ReadNumber(argv[5]),ReadNumber(argv[6]));
    SaveImage(argv[8],Img1);
    
  } else if (strcmp(argv[1],"removeborders")==0) {
    CheckUsage(argc,8,"removeborders","left right bottom top in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[6]);
    removeborders(ReadNumber(argv[2]),ReadNumber(argv[3]),ReadNumber(argv[4]),ReadNumber(argv[5]));
    SaveImage(argv[7],Img1);
    
  } else if (strcmp(argv[1],"add")==0) {
    CheckUsage(argc,5,"add","in1.g in2.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    Img1=AllocMemAndLoadGreyImage(argv[3]);
    add();
    SaveImage(argv[4],Img2);    
  
   } else if (strcmp(argv[1],"addI")==0) {
    CheckUsage(argc,5,"addI","value in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[3]);
    addI(ReadNumber(argv[2]));
    SaveImage(argv[4],Img2);    
    
   } else if (strcmp(argv[1],"sub")==0) {
    CheckUsage(argc,5,"sub","in1.g in2.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    Img1=AllocMemAndLoadGreyImage(argv[3]);
    sub();
    SaveImage(argv[4],Img2);    
  
   } else if (strcmp(argv[1],"subI")==0) {
    CheckUsage(argc,5,"subI","in.g value out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    subI(ReadNumber(argv[3]));
    SaveImage(argv[4],Img2);    
    
   } else if (strcmp(argv[1],"rsbI")==0) {
    CheckUsage(argc,5,"rsbI","value in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[3]);
    rsbI(ReadNumber(argv[2]));
    SaveImage(argv[4],Img2);    
    
  } else if (strcmp(argv[1],"mul")==0) {
    CheckUsage(argc,5,"mul","in1.g in2.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    Img1=AllocMemAndLoadGreyImage(argv[3]);
    mul();
    SaveImage(argv[4],Img2);     
    
  } else if (strcmp(argv[1],"mulI")==0) {
    CheckUsage(argc,5,"mulI","factor in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[3]);
    mulI(ReadNumber(argv[2]));
    SaveImage(argv[4],Img2);     
    
  } else if (strcmp(argv[1],"div")==0) {
    CheckUsage(argc,5,"div","in1.g in2.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    Img1=AllocMemAndLoadGreyImage(argv[3]);
    division();
    SaveImage(argv[4],Img2);
    
  } else if (strcmp(argv[1],"divI")==0) {
    CheckUsage(argc,5,"divI","in.g divisor out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    divI(ReadNumber(argv[3]));
    SaveImage(argv[4],Img2);

  } else if (strcmp(argv[1],"mknoise")==0) {    
    CheckUsage(argc,7,"mknoise","min max width height out.g");
    Img0=AllocMemForGreyImage(ReadNumber(argv[4]),ReadNumber(argv[5]));
    mknoise(ReadNumber(argv[2]),ReadNumber(argv[3]));
    SaveImage(argv[6],Img0);

  } else if (strcmp(argv[1],"hmirror")==0) {    
    CheckUsage(argc,4,"hmirror","in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    hmirror();
    SaveImage(argv[3],Img2);

  } else if (strcmp(argv[1],"vmirror")==0) {    
    CheckUsage(argc,4,"vmirror","in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    vmirror();
    SaveImage(argv[3],Img2);

  } else if (strcmp(argv[1],"rotate90")==0) {    
    CheckUsage(argc,4,"rotate90","in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    rotate90();
    SaveImage(argv[3],Img2);
    
  } else if (strcmp(argv[1],"rotate180")==0) {    
    CheckUsage(argc,4,"rotate180","in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    rotate180();
    SaveImage(argv[3],Img2);
    
  } else if (strcmp(argv[1],"rotate270")==0) {    
    CheckUsage(argc,4,"rotate270","in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[2]);
    rotate270();
    SaveImage(argv[3],Img2);

  } else if (strcmp(argv[1],"edge")==0) {    
    CheckUsage(argc,7,"edge","size emptycol edgecol in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[5]);
    edge(ReadNumber(argv[2]),ReadNumber(argv[3]),ReadNumber(argv[4]));
    SaveImage(argv[6],Img0);

  } else if (strcmp(argv[1],"explode")==0) {
    CheckUsage(argc,5,"explode","nbpoints in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[3]);
    explode(ReadNumber(argv[2]));
    SaveImage(argv[4],Img0);
             
  } else if (strcmp(argv[1],"pal2rgb")==0) {
    CheckUsage(argc,7,"pal2rgb","in.pal in.g red.g green.g blue.g");
    Img0=AllocMemAndLoadGreyImage(argv[3]);
    pal2rgb(argv[2]);
    SaveImage(argv[4],Img0);
    SaveImage(argv[5],Img1);
    SaveImage(argv[6],Img2);

  } else if (strcmp(argv[1],"setmincol")==0) {
    CheckUsage(argc,5,"setmincol","mincol in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[3]);
    setmincol(ReadNumber(argv[2]));
    SaveImage(argv[4],Img0);
        
  } else if (strcmp(argv[1],"setmaxcol")==0) {
    CheckUsage(argc,5,"setmaxcol","maxcol in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[3]);
    setmaxcol(ReadNumber(argv[2]));
    SaveImage(argv[4],Img0);
    
  } else if (strcmp(argv[1],"rectangle")==0) {
    CheckUsage(argc,9,"rectangle","top bottom left right color in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[7]);
    rectangle(ReadNumber(argv[2]),ReadNumber(argv[3]),ReadNumber(argv[4]),ReadNumber(argv[5]),ReadNumber(argv[6]));
    SaveImage(argv[8],Img0);
    
  } else if (strcmp(argv[1],"rectangle2")==0) {
    CheckUsage(argc,9,"rectangle2","top left length height color in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[7]);
    rectangle2(ReadNumber(argv[2]),ReadNumber(argv[3]),ReadNumber(argv[4]),ReadNumber(argv[5]),ReadNumber(argv[6]));
    SaveImage(argv[8],Img0);
    
  } else if (strcmp(argv[1],"landscape")==0) {
    CheckUsage(argc,7,"landscape","nbiterations fractalfactor*100 absorption*100 in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[5]);
    landscape(ReadNumber(argv[2]),ReadNumber(argv[3]),ReadNumber(argv[4]));
    SaveImage(argv[6],Img1);

  } else if (strcmp(argv[1],"between")==0) {
    CheckUsage(argc,6,"between","colmin colmax in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[4]);
    between(ReadNumber(argv[2]),ReadNumber(argv[3]));
    SaveImage(argv[5],Img0);

  } else if (strcmp(argv[1],"range")==0) {
    CheckUsage(argc,6,"range","colmin colmax in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[4]);
    range(ReadNumber(argv[2]),ReadNumber(argv[3]));
    SaveImage(argv[5],Img0);

  } else if (strcmp(argv[1],"disc")==0) {
    CheckUsage(argc,8,"disc","x y r color in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[6]);
    disc(ReadNumber(argv[2]),ReadNumber(argv[3]),ReadNumber(argv[4]),ReadNumber(argv[5]));
    SaveImage(argv[7],Img0);

  } else if (strcmp(argv[1],"fillwithOs")==0) {
    CheckUsage(argc,11,"fillwithOs","OutRMin OutRMax InRMmin InRMmax Border ClearColor CircleColor in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[9]);
    fillwithOs(ReadNumber(argv[2]),ReadNumber(argv[3]),ReadNumber(argv[4]),ReadNumber(argv[5]),ReadNumber(argv[6]),ReadNumber(argv[7]),ReadNumber(argv[8]));
    SaveImage(argv[10],Img0);
        
  } else if (strcmp(argv[1],"floyd")==0) {
    CheckUsage(argc,5,"floyd","nb_colors in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[3]);
    floyd(ReadNumber(argv[2]));
    SaveImage(argv[4],Img0);
    
  } else if (strcmp(argv[1],"subsample")==0) {
    CheckUsage(argc,6,"subsample","xratio yratio in.g out.g");
    Img0=AllocMemAndLoadGreyImage(argv[4]);
    subsample(ReadNumber(argv[2]),ReadNumber(argv[3]));
    SaveImage(argv[5],Img1);

  } else if (strcmp(argv[1],"mkstrippedtube")==0) {
    CheckUsage(argc,6,"mkstrippedtube","diameter nbsteps nbstripes col0 col1 out.g");
    Img0=AllocMemForGreyImage(ReadNumber(argv[2]),ReadNumber(argv[3]));
    mkstrippedtube(ReadNumber(argv[2]),ReadNumber(argv[3]),ReadNumber(argv[4]),ReadNumber(argv[5]),ReadNumber(argv[6]));
    SaveImage(argv[7],Img0);

  } else if (strcmp(argv[1],"grey2bitplanes")==0) {
    CheckUsage(argc,3,"grey2bitplanes","nbbitplanes in.g out.xxx");
    Img0=AllocMemAndLoadGreyImage(argv[3]);
    grey2bitplanes(ReadNumber(argv[2]),argv[4]);

  } else {
    printf("Error: %s is not a valid function\nSee AllFunctions.txt\n",argv[1]); exit(1);
  }
  
  printf(" OK");
  exit(0);
}