#include <stdio.h>
#include <math.h>
FILE *fp;
char head[72][80];
char c256[512*512];
char b[512*512*2];
int  i,n1,n2,c1,m1,m2,mm,nn,k=0;
char cl[12],cb[12];
short  a[4096*4096],a2[4160];

main(ac,av)
int ac; char *av[];
{
  char f1[60],*b2;
  int j,kk;
  int da[8]={51,10,19,20,3,4,36,42};
  float white;
  if(ac<2){
    printf("\n\n\t *** output GIF&txt (512*512) (256*256) ***\n");
    printf("\n\t Ex:   gh8 in(d*.ccd)      2006,10\n\n");
    exit(0);
  }
  fp=fopen(av[1],"rb");
  if(fp==0)exit(0);
  fread(head,72,80,fp);
  sscanf(&head[3][20],"%d",&n1);
  sscanf(&head[4][20],"%d",&n2);
  sscanf(&head[5][20],"%d",&c1);
  m1=m2=512; mm=(n2-1)/m1+1; 
  kk=0;
  for(i=0;i<n2;i++){
    fread(a2,2,n1,fp);
    baseline(a2,n1,n2,c1);    
    for(j=0;j<n2;j++)a[kk++]=a2[j];
  }	    
  fclose(fp);
  if(n1>4000){
    nn=mm*m1*mm*m1*2;
    b2=(char *)malloc(nn);
    fp=fopen("/u/ccdev/flati2","rb");
    fread(b2,1,nn,fp);
    fread(&white,1,4,fp);
    fclose(fp);
    do_flat(a,b2,nn/2,white/2);
    free(b2);
    balence(a,b);
  }  
  shrink2(a,b); 

  fp=fopen("/u/ccdev/status/lastimage.txt","w");
  fprintf(fp,"%s\n",av[1]);
  for(i=0;i<8;i++){
    head[da[i]][79]=0;
    fprintf(fp,"%s\n",&head[da[i]][0]);
  }
  fclose(fp);

  comp4(b,c256);
  fp=fopen("last2.gif","wb");
  jcomp();
  fclose(fp);

  for(i=0;i<m1*m2;i++)b[i]=c256[i];
  shrink256(b,c256);
  m1/=2; m2/=2;		
  fp=fopen("last1.gif","wb");
  jcomp();
  fclose(fp);

  system("convert last1.gif /u/ccdev/status/last1.jpg");
  system("convert last2.gif /u/ccdev/status/last2.jpg");
}

do_flat(a,b,n,white)
short a[],b[];
int n;
float white;
{
  int i;
  float x,y;
  for(i=0;i<n;i++){ x=a[i]*white/b[i];
        if(x>32767.)x=32767.;
        if(x<0.)x=0.;
        a[i]=x;
  }
}

balence(a,b)
short a[],b[];
{
  int j;
  float white,black,sigma;  
  shrink1(a,b,1536,2560,1792,2048);
  white_(b,1024,256,&white,&sigma); 
//    printf("%f\n",white);
  shrink1(a,b,1536,2560,2048,2304);
  white_(b,1024,256,&black,&sigma); 
//    printf("%f\n",black);
  j=white-black+1.5;
  shrinkj(a,j);
}

shrink1(b,bb,i1,i2,j1,j2)
short b[4096][4096],bb[1024*256];
int i1,i2,j1,j2;
{
  int i,j,k=0;
  for(i=i1;i<i2;i++)for(j=j1;j<j2;j++)bb[k++]=b[i][j];
}

shrinkj(a,j)
short a[4096][4096];
int j;
{
  int i,k;
  for(i=0;i<4096;i++)for(k=0;k<2048;k++)a[i][k]-=j;
}

shrink256(a,b)
char a[m1][m2],b[m1/2][m2/2];
{
  int i,j;
  for(j=0;j<m2;j+=2)
  for(i=0;i<m1;i+=2)
  b[i/2][j/2]=a[i][j];
}

shrink2(a,b)
short a[],b[];
{
  int i,j,k=0,i1,j1,k1,x,y;
  x=mm*mm;
  for(i=0;i<m2*mm;i+=mm)for(j=0;j<m1*mm;j+=mm){
    y=0;
    for(i1=0;i1<mm;i1++){
      k1=(i1+i)*n2+j;
      for(j1=0;j1<mm;j1++)y+=a[k1+j1];
    }
    b[k++]=y/x; 
  }
}

comp4(b,c)
short b[]; unsigned char c[];
{
  float white,black,d,sigma;
  short bb[512*512];
  int i,j;
  if(m2>n2)movepic(b,bb,m2,n2);
  if(n2>=4000 || m2>n2){
    shrink0(b,bb);
    white_(bb,256,256,&white,&sigma); 
  } else white_(b,m1,m2,&white,&sigma);
  black=white+20.*sigma;
  white-=20.*sigma;
  d=250./(black-white);
  for(j=0;j<m1*m2;j++){
    i=(b[j]-white)*d+0.5;
    if(i<0)i=0;
    if(i>255)i=255;
    c[j]=i;
  }
}

movepic(b,bb,m1,n1)
short b[512][512],bb[512][512];
int m1,n1;
{
  int i,j,k;
  for(i=0;i<512;i++)for(j=0;j<512;j++){ bb[i][j]=b[i][j]; b[i][j]=0; }
  k=(m1-n1)/2;
  for(i=0;i<n1;i++)for(j=0;j<n1;j++) b[i+k][j+k]=bb[i][j];
}

shrink0(b,bb)
short b[512][512],bb[256*256];
{
  int i,j,k=0;
  for(i=128;i<384;i++)for(j=128;j<384;j++)bb[k++]=b[i][j];
}

white_(map,n1,n2,peak,sigma)
short map[];
int n1,n2;
float *peak,*sigma;
{
  int *ihist;
  ihist=(int*)malloc(32767*4);
  white_1(map,n1,n2,peak,sigma,ihist);
  free(ihist);
}

white_1(map,n1,n2,xpeak,sigma,ihist)
short map[]; int n1,n2;
float *xpeak, *sigma;
long  ihist[];
{
  int i,ii;
  int imax=0;
  for(i=0;i<32767;i++)ihist[i]=0;
  for(i=0;i<n1*n2;i++){
    ii=map[i];
    ii+=101;
    if(ii<0 || ii >= 32767)continue;
    ihist[ii]++;
    if(ii>imax)imax=ii;
  }
  histt(ihist,xpeak,sigma,imax);
  *xpeak-=100.0;
}

histt(ihist,xpeak,fsigma,n)
long ihist[];
float *xpeak,*fsigma;
int   n;
{
  int ihtot=0, ihsum=0, iwid=0, maxh=0, mode=0;
  int i,ii,j,jj,k,m,nn, isigm, nfilt, noff, ixpeak, ilow, ihih,ilim,icount;
  float xpar[101],buf[201];
  float xmean=0., xmed=0, sigmsq=0., xcount=0;
  float rttpi, conv, xx, temp, xp,fs,gs,cogd,cogn,cog,xnumb,sd;

  rttpi=2.506628275;
  for(i=1;i<=n;i++){
    nn=n+1-i;
    if(ihist[nn]!=0)break;
  }
  for(i=1;i<=nn;i++)ihtot+=ihist[i];
  for(i=1;i<=nn;i++){
    ihsum+=ihist[i];
    if(ihsum < 0.1*ihtot || ihsum > 0.9*ihtot)continue;
    if(ihist[i] <= maxh)continue;
    maxh=ihist[i];
    mode=i;
  }
  for(i=1;i<=nn;i++){
    temp=ihist[i];
    xcount+=temp;
    xmean+=temp*i;
    sigmsq+=temp*i*i;
  }
  for(i=1;i<=nn;i++){
    xmed+=ihist[i];
   if(xmed > 0.5*xcount){ xmed=i; break; }
  }
  xmean/=xcount;
  gs=xcount/(maxh*rttpi);
  if(gs > 0.1*nn)gs=0.1*nn;
  isigm=gs+0.5;
  if(isigm<3)isigm=3;
  if(isigm>50)isigm=50;
  nfilt=isigm*4+1;
  noff=(nfilt+1)/2;
  conv=0.5/(float)(isigm*isigm);
  for(i=1;i<=noff;i++){
    xx=i-noff;
    buf[i]=exp(-conv*xx*xx);
    buf[nfilt+1-i]=buf[i];
  }
  xp=mode;
  for(k=1;k<=2;k++){
    ixpeak=xp+0.5;
    if(ixpeak > nn)ixpeak=nn;
    if(ixpeak <1 )ixpeak=1;
    ilow=ixpeak-isigm;
    if(ilow <1 )ilow=1;
    ihih=ixpeak+isigm;
    if(ihih >nn)ihih=nn;
    m=ihih-ilow+1;
    ii=1;  cogd=cogn=0.;
    for(i=ilow;i<=ihih;i++){
      cogd+=ihist[i];
      cogn+=ihist[i]*i;
      temp=xnumb=0.;
      for(j=1;j<=nfilt;j++){
	jj=j-noff;
	if(i+jj < 1 || i+jj > nn)continue;
	temp+=ihist[i+jj]*buf[j];
	xnumb+=buf[j];
      }
      temp = (xnumb > 0.) ? temp/=xnumb : 0. ;
      if(temp < 1.)temp=1.;
      xpar[ii++]=log((double)temp);
    }
    parbol(xpar,m,&xp,&sd);
    xx=sd*sd-isigm*isigm;
    if(xx < 0.)xx=0.;
    fs=sqrt((double)xx);
    xp+=ilow-1;
    cog=(cogd > .5) ? cogn/cogd : xp;
    xx=xp-mode; if(xx <0.)xx=-xx;
    if(xx > gs)xp=cog;
  }
  mode--;
  xp=xp-1.0;
  xmean=xmean-1.0;
  xmed=xmed-1.0;
  if(xp > nn-gs)xp=xmed;
  if(xp < gs)xp=xmed;
  ixpeak=xp+1.5;
  icount=ihist[ixpeak];
  ilim=0.68*xcount+0.5;
  for(i=1;i<=nn;i++){
    ii=ixpeak+i;
    if(ii <= nn){
      iwid++;
      icount+=ihist[ii];
      if(icount > ilim)break;
    }
    ii=ixpeak-i;
    if(ii >= 1){
      iwid++;
      icount+=ihist[ii];
      if(icount > ilim)break;
    }
  }
  gs=0.5*iwid;
  if(gs<fs)fs=gs;
  if(fs<1.)fs=1.;
  *xpeak=xp;
  *fsigma=fs;
}

#define MAXCODES 4096
#define TABLESIZES 4999
unsigned char SuffixTable[4999], ByteBuf[2592], BlockBuf[256];
int Encode, RunBits, MaxCodeSize, ByteCount, ShiftBits;
int Code, Width, Height, Dots, Rows;
int PrefixCode, SuffixCode, RandomIndex, Val;
int EncodeTable[4999], PrefixTable[4999];
short head1[7]={0x4947,0x3846,0x6137,0,0,0xf7,0};
short head2[5]={0,0,0,0,0x800};
unsigned long TempCode,jj;
jcomp()
{	
  head1[3]=head2[2]=Width=m1;
  head1[4]=head2[3]=Height=m2;
  fwrite(head1,1,13,fp);
  for(i=0;i<256;i++){
    ByteBuf[0]=ByteBuf[1]=ByteBuf[2]=i;
    fwrite(ByteBuf,1,3,fp);
  }
  fputc(0x2c,fp);
  fwrite(head2,2,5,fp);
  ByteCount = ShiftBits = TempCode = 0;
  ClearTable();
  Code = 256;
  FillBlockBuf();
  for(Rows = 0; Rows < Height; Rows++)  {
	  jj=(m1-Rows)*m1;
    for(i=0;i<m1;i++)ByteBuf[i]=c256[--jj];	  
    if(Rows == 0) { PrefixCode = ByteBuf[0];  Dots = 1; }
    else Dots = 0;
    while( Dots < Width)  {
      SuffixCode = ByteBuf[Dots++];
      RandomIndex = PrefixCode ^ (SuffixCode << 4);
      if(RandomIndex == 0)  Val = 1;
      else Val = TABLESIZES - RandomIndex;
      while(1)  {
        if(EncodeTable[RandomIndex] == 0)  {
          Code = PrefixCode;
          FillBlockBuf();
          if(Encode == MAXCODES)  {
            Code = 256;
            FillBlockBuf();
            ClearTable();
          }
          else  {
            if(Encode == MaxCodeSize)  {
              MaxCodeSize <<= 1;
              RunBits++;
            }
            PrefixTable[RandomIndex] = PrefixCode;
            SuffixTable[RandomIndex] = SuffixCode;
            EncodeTable[RandomIndex] = Encode++;
          }
          PrefixCode = SuffixCode;
          break;
        }
        if(PrefixTable[RandomIndex] == PrefixCode &&
           SuffixTable[RandomIndex] == SuffixCode)  {
          PrefixCode = EncodeTable[RandomIndex];
          break;
        }
        else  {
          RandomIndex -= Val;
          if(RandomIndex < 0)  RandomIndex += TABLESIZES;
        }
      }
    }
  }
  Code = PrefixCode;
  FillBlockBuf();
  Code = 257;
  FillBlockBuf();
  if(ShiftBits > 0 || ByteCount > 0)   {
    BlockBuf[++ByteCount] = TempCode & 0x00FF;
    SaveCode();
  }
  fputc(0, fp);  fputc(';', fp);
}

ClearTable()
{
  Encode = 258;  RunBits = 9;  MaxCodeSize =512;
  for(i = 0; i < TABLESIZES; i++) EncodeTable[i] = 0;
}

FillBlockBuf()
{
  TempCode |= (unsigned long)Code << ShiftBits;
  ShiftBits += RunBits;
  while(ShiftBits >= 8)  {
    BlockBuf[++ByteCount] = TempCode & 0x00FF;
    if(ByteCount == 255)  SaveCode();
    TempCode >>= 8;
    ShiftBits -= 8;
  }
}

SaveCode()
{
  BlockBuf[0] = ByteCount;
  fwrite(BlockBuf, 1, ByteCount+1, fp);
  ByteCount = 0;
}

baseline(a,n1,n2,c1)
short a[];
int n1,n2,c1;
{
  int i,base1,base2;
  swap2(a,n1*2);
  for(i=0;i<n1;i++)a[i]=(a[i]+32768)>>1;
// neednot concer of BZERO, a[]'s bzero =0
  base2=median(&a[n1-32],32);
  if(n1-n2==64){
    base1=median(&a[n1-64],32);
    for(i=0;i<n2;i++)if((i+c1)<2048) a[i]-=base1; else a[i]-=base2;
  } else for(i=0;i<n2;i++)a[i]-=base2;
} 

median(a,n)
short a[];
int   n;
{
  int i,j;
  short k;
  for(i=0;i<n-1;i++)for(j=i+1;j<n;j++)if(a[i]>a[j]){
    k=a[j]; a[j]=a[i]; a[i]=k;
  }
  i=a[n/2];
  return(i);
}