unit2.h

it will work it is recogntion code only

    int made1=sscanf(line,"%d %d %d",num,num+1,num+2); 
    int i=0; for(i=0;i<made1;i++) *(numbers+allmade+i)=*(num+i);
    allmade += made1;
  } 
  *pxSize=numbers[0];*pySize=numbers[1];int Ngray=numbers[2];
  if (Ngray>255) 
  {printf("\n\7 ???File '%s' is not 255 gray image\n",fullname);exit(1);}  

  int rSize=(*pxSize)*(*pySize);
  if (rSize>rSize0) // rSize0 should be >0 from the beginning.
  {delete[]retina2D; retina2D=new unsigned char [rSize];} 
  //unsigned char *ret=retina2D;
  int u;
  if(P EQ 2)for(int i=0;i<rSize;i++)
                {
                fscanf(pgmFile,"%d",&u);//fscanf(pgmFile,"%d",retina2D+i);
                retina2D[i]= (unsigned char) (u);
                }  
  if(P EQ 5) fread (retina2D, 1, rSize, pgmFile);        
  fclose(pgmFile); 
  return(retina2D);
};//END.
//%%%%%%%%%%%%%%%%%%%%%%% write and read integer text P2 (my style)

void pgmP2WriteInt(char *fullname,//Writes integer text P2 (my style)
                  int xSize, int ySize,
                  int *retina2D,int maxGray)
{  FILE *pgmPict;
pgmPict=fopen(fullname,"w"); 
fprintf(pgmPict,"P2\n");fprintf(pgmPict,"%d %d\n%d\n",xSize,ySize,maxGray);
int rSize=xSize*ySize;
int i;for(i=0;i<rSize;i++) {fprintf(pgmPict,"%d ",retina2D[i]);}
fclose(pgmPict);
};//END.

int *pgm2Readint(char *fullname,//Reads integer text P2 (my style) 
      /*RESULTS:*/int *pxSize,int *pySize)
  { 
  FILE *pg2Pict; char shouldBeP2[2]; int shouldBe255;
  pg2Pict=fopen(fullname,"r");
  fscanf(pg2Pict,"%2s\n",shouldBeP2);
  char line[100];                              
  fscanf(pg2Pict,"%[^\n]\n",line);              
  while(line[0]=='#')                           
  {fscanf(pg2Pict,"%[^\n]\n",line);}           
  sscanf(line,"%d %d",pxSize,pySize);            
  fscanf(pg2Pict,"%d\n",&shouldBe255);
  int rSize=(*pxSize)*(*pySize);
  int*retina2D=new int[rSize];
  int i;for(i=0;i<rSize;i++)fscanf(pg2Pict,"%d",&(retina2D[i]));
  fclose(pg2Pict); 
  return(retina2D);
  }//END.

unsigned char *Read_uchar(char *fullname,int rSize)//Reads uchar
  { 
  FILE *uFile; 
  if( (uFile  = fopen( fullname, "r" )) == NULL ) 
    {printf("\n\7 ???File '%s' was not opened\n",fullname );exit(1);}  
  
  unsigned char*retina2D=new unsigned char[rSize];int u;
  int i;for(i=0;i<rSize;i++)
     {fscanf(uFile,"%d",&u);//fscanf(uFile,"%d",&(retina2D[i]));
      retina2D[i]=u;}
  fclose(uFile); 
  return(retina2D);
  }//END. 
double *Read_double(char *fullname,int rSize)//Reads uchar
  { 
  FILE *dFile; 
  if( (dFile  = fopen( fullname, "r" )) == NULL ) 
    {printf("\n\7 ???File '%s' was not opened\n",fullname );exit(1);}  
  
  double*retina2D=new double[rSize]; float u;
  int i;for(i=0;i<rSize;i++)
      { fscanf(dFile,"%f",&u);retina2D[i]=u;}

  fclose(dFile); 
  return(retina2D);
  }//END. 
//%%%%%%%%%%%%%%%%%%%%%%% weighted median:
int medianW(int Larr,unsigned char*arr,int Marr,double*weight,double wei_threshold,double*medtemp)
// all arrays at least Larr length
// arr considt of numbers >=0, <=Marr,(Mar<256)
{ int i,j;
  for (i=Marr-1;i>=0;i--) medtemp[i]=0;
  for (j=Larr-1;j>=0;j--) {i=arr[j];medtemp[i]+=weight[i];}
  double cumsum=0;
  for (i=0;i<Marr;i++) {cumsum+=medtemp[i];if (cumsum>=wei_threshold) break;}
  return(i);  
}//END
int medianW1(int Larr,unsigned char*arr,int Marr,double wei_threshold,double*medtemp)
// weight = 1 
// all arrays at least Larr length
// arr considt of numbers >=0, <=Marr,(Mar<256)
{ int i,j;
  for (i=Marr-1;i>=0;i--) medtemp[i]=0;
  for (j=Larr-1;j>=0;j--) {i=arr[j];medtemp[i]++;}
  double cumsum=0;
  for (i=0;i<Marr;i++) {cumsum+=medtemp[i];if (cumsum>=wei_threshold) break;}
  return(i);  
}//END

//%%%%%%%%%%%%%%%%%%%%%%%  sort
/*
const int MLA=1128;
int *c1H=new int[MLA],*c2H=new int[MLA];
//int c1H[MLA],c2H[MLA];
*/
int b0,b1,b2,b3;//- globals b0,b1,b2,b3 for output of sort4
void sort4(int a0,int a1,int a2,int a3)

{
if(a0<=a1)
     if(a2<=a3)
         if(a1<=a3)
             if(a1<=a2){b0=a0;b1=a1;b2=a2;b3=a3;}
             else if(a0<=a2){b0=a0;b1=a2;b2=a1;b3=a3;}
                  else {b0=a2;b1=a0;b2=a1;b3=a3;}
 
         else
             if(a0<=a2){b0=a0;b1=a2;b2=a3;b3=a1;}
             else if(a0<a3){b0=a2;b1=a0;b2=a3;b3=a1;}
                  else {b0=a2;b1=a3;b2=a0;b3=a1;}
     else//(a2<=a3)
         if(a1<=a2)
             if(a1<=a3){b0=a0;b1=a1;b2=a3;b3=a2;}
             else if(a0<a3){b0=a0;b1=a3;b2=a1;b3=a2;}
                  else {b0=a3;b1=a0;b2=a1;b3=a2;}
 
         else
             if(a0<=a3){b0=a0;b1=a3;b2=a2;b3=a1;}
             else if(a0<=a2){b0=a3;b1=a0;b2=a2;b3=a1;}
                  else {b0=a3;b1=a2;b2=a0;b3=a1;}

else//if(a0<=a1)
     if(a2<=a3)
         if(a0<=a3)
             if(a0<=a2){b0=a1;b1=a0;b2=a2;b3=a3;}
             else if(a1<=a2){b0=a1;b1=a2;b2=a0;b3=a3;}
                  else {b0=a2;b1=a1;b2=a0;b3=a3;}
 
         else
             if(a1<=a2){b0=a1;b1=a2;b2=a3;b3=a0;}
             else if(a1<=a3){b0=a2;b1=a1;b2=a3;b3=a0;}
                  else {b0=a2;b1=a3;b2=a1;b3=a0;}
     else//(a2<=a3)
         if(a0<=a2)
             if(a0<=a3){b0=a1;b1=a0;b2=a3;b3=a2;}
             else if(a1<=a3){b0=a1;b1=a3;b2=a0;b3=a2;}
                  else {b0=a3;b1=a1;b2=a0;b3=a2;}
 
         else
             if(a1<=a3){b0=a1;b1=a3;b2=a2;b3=a0;}
             else if(a1<=a2){b0=a3;b1=a1;b2=a2;b3=a0;}
                  else {b0=a3;b1=a2;b2=a1;b3=a0;}

}

int* sortup(int *A,int LA,int*sorttemp)//sorttemp is twice as LA
                 // A is not damaged. Output=pointer to a part of sorttemp
                 // globals b0,b1,b2,b3 for output of sort4
{
int *c1,*c2,*c3;//int *Cout;
c1=A;c2=sorttemp; 
register int  *a,*b,*c;
c=c2;int k;int *Ak=A-1,*ck=c;int LA4=LA/4; 
for (k=LA4;k>0;k--)
{
//sort4(A[k],A[k+1],A[k+2],A[k+3]);c[k]=b0;c[k+1]=b1;c[k+2]=b2;c[k+3]=b3;
  sort4(*++Ak,*++Ak,*++Ak,*++Ak);*ck++=b0;*ck++=b1;*ck++=b2;*ck++=b3;
}
int rest=LA-LA4*4,a0,a1,a2;
if (rest>0)
{
   if (rest EQ 1) *ck++=*++Ak;
   else
     if (rest EQ 2){a0=*++Ak;a1=*++Ak;
        if(a0<a1){*ck++=a0;*ck++=a1;}else{*ck++=a1;*ck++=a0;}}
     else //rest EQ 3
     {
         {a0=*++Ak;a1=*++Ak;a2=*++Ak;}
         if (a0<a1) 
             if (a1<a2){*ck++=a0;*ck++=a1;*ck++=a2;}
             else
                if (a0<a2){*ck++=a0;*ck++=a2;*ck++=a1;}
                else {*ck++=a2;*ck++=a0;*ck++=a1;}
         else//(a0<a1) 
             if (a0<a2){*ck++=a1;*ck++=a0;*ck++=a2;}
             else
                if (a1<a2){*ck++=a1;*ck++=a2;*ck++=a0;}
                else {*ck++=a2;*ck++=a1;*ck++=a0;}
     }
}

if (LA<5) return (c2);
c1=sorttemp;c2=sorttemp+LA;

int N=4,N2=8,k1,k2,k22;
for(;;)
{       k=0,k1=k+N,k2=k+N2;if (k2>LA) k2=LA;
    for (; k1<k2; )
    {  int *a0=c1+k,*b0=c1+k1; 
       int i,j=k2-k1-1;c=c2+k+N+j;a=a0+N-1;b=b0+j;
       int free=N; if (j+1<N) free=j+1;
   
            {register int u=*a,v=*b;int it=free;
               if(u>v)  do{*c--=u;u=*--a;it--;}while(it && (u>=v));
                 else   do{*c--=v;v=*--b;it--;}while(it && (u<=v));
               for (;it>0;it--) 
                  if(u>v){*c--=u;u=*--a;}else{*c--=v;v=*--b;} 
            }
            
            {int ir=a-a0,jr=b-b0;register int u=a0[ir],v=b0[jr];

#define SEVENSTYLE 
#ifdef  SEVENSTYLE          
             while (ir>=0 && jr>=0) 
                      if(u>v){*c--=u;u=a0[--ir];}else{*c--=v;v=b0[--jr];};
#else
             if (ir<jr)
               {while (ir>=0 && jr>=0) 
                        if(u>v){*c--=u;u=a0[--ir];}else{*c--=v;v=b0[--jr];};}
             else
               {while (jr>=0 && ir>=0) 
                        if(u>v){*c--=u;u=a0[--ir];}else{*c--=v;v=b0[--jr];};} 
#endif
             i=ir;j=jr;
            }
    //if (i<0) while(j>=0) *c--=b0[j--];
    //else if (j<0) while(i>=0) *c--=a0[i--];
    if (i<0) { register int jr=j; while(jr>=0) *c--=b0[jr--]; }
    else if (j<0) { register int ir=i; while(ir>=0) *c--=a0[ir--]; }
    k22=k2;
    k+=N2,k1+=N2,k2+=N2; if (k2>LA) k2=LA;
    }

//Elements *(c1+k2),..,*(c1+LA-1)  copy to adresses  c2+k2,..,c2+LA-1
int numtail=LA-k22;
if (numtail>0){c=c2+LA-1;a=c1+LA;register int t=numtail;while(t-->0)*c--=*--a;}

    N = N2; N2=N2<<1;
    c3=c2;c2=c1;c1=c3;
    if (N>=LA) return(c1);  
}//for(;;;)
}//END.
//////////////////////////////////////////////////////////
int goodfft(int N)
{ 
 int i,n=round(ceil(log(N)/log(2.)+0.00001)),Ngood=1<<n,Ntry; 
 if (n>0) if ( N EQ (1<<(n-1)) ) return(N);
 int N10=round(ceil(N*1.1));
 for (i=0;i<=n;i++) 
 { Ntry=(1<<i)*3; 
   if (Ntry>=N)
     if (Ntry<Ngood) {Ngood=Ntry;break;}
     else break;     
 } if (Ngood<N10) return(Ngood);
 for (i=0;i<=n;i++) 
 { Ntry=(1<<i)*9; 
   if (Ntry>=N)
     if (Ntry<Ngood) {Ngood=Ntry;break;}
     else break;     
 } if (Ngood<N10) return(Ngood);
 for (i=0;i<=n;i++) 
 { Ntry=(1<<i)*27; 
   if (Ntry>=N)
     if (Ntry<Ngood) {Ngood=Ntry;break;}
     else break;     
 } if (Ngood<N10) return(Ngood);
 for (i=0;i<=n;i++) 
 { Ntry=(1<<i)*81; 
   if (Ntry>=N)
     if (Ntry<Ngood) {Ngood=Ntry;break;}
     else break;     
 } if (Ngood<N10) return(Ngood);
 for (i=0;i<=n;i++) 
 { Ntry=(1<<i)*5; 
   if (Ntry>=N)
     if (Ntry<Ngood) {Ngood=Ntry;break;}
     else break;     
 } if (Ngood<N10) return(Ngood);
 for (i=0;i<=n;i++) 
 { Ntry=(1<<i)*15; 
   if (Ntry>=N)
     if (Ntry<Ngood) {Ngood=Ntry;break;}
     else break;     
 } if (Ngood<N10) return(Ngood);

 for (i=0;i<=n;i++) 
 { Ntry=(1<<i)*25; 
   if (Ntry>=N)
     if (Ntry<Ngood) {Ngood=Ntry;break;}
     else break;     
 } if (Ngood<N10) return(Ngood);
for (i=0;i<=n;i++) 
 { Ntry=(1<<i)*45; 
   if (Ntry>=N)
     if (Ntry<Ngood) {Ngood=Ntry;break;}
     else break;     
 } if (Ngood<N10) return(Ngood);

  for (i=0;i<=n;i++) 
 { Ntry=(1<<i)*125; 
   if (Ntry>=N)
     if (Ntry<Ngood) {Ngood=Ntry;break;}
     else break;     
 }
return Ngood;
}//END.

void normalising(double maxAngle,double beta,int N2pi,double*c,double uomega,
                 double *cN)// normalises c into cN; cN should come allocated
                            // c is a circus on N2pi length
                            // uomega is of c.
                            //for associated circus c uomega=-1 (it is also clled "order")
{ double dfi=2*pi/N2pi;

if (beta EQ 0) { for(int i=0;i<N2pi;i++) cN[i]=0; return; }

for(int i=0;i<N2pi;i++)
        {double fi=i*dfi;  
         double fi1=fi-maxAngle, co=cos(fi1)/beta,si=sin(fi1)*beta,
         fiinci=maxAngle+atan2(si,co), K=1.0/sqrt(co*co+si*si);
         fiinci-=floor(fiinci/pi2)*pi2;// go to [0,2pi)                  
         double jfiinci=fiinci/pi2*N2pi;int j=int(floor(jfiinci));
         double x=jfiinci-j;
         //double *c=circ[nf]+cf*N2pi;//- pointer to a given circus NoNo=(nf,cf) 
                // double *c=circ..;// or circA..
       
         double ciinj; // a little bit unefficient:
         if (0<j && j<N2pi-2) ciinj=spline4(x,c[j-1],c[j],c[j+1],c[j+2]);
           else
           {if (j EQ 0) ciinj=spline4(x,c[N2pi-1],c[0],c[1],c[2]);
            if (j EQ N2pi-1) ciinj=spline4(x,c[N2pi-2],c[N2pi-1],c[0],c[1]);
            if (j EQ N2pi-2) ciinj=spline4(x,c[N2pi-3],c[N2pi-2],c[N2pi-1],c[0]);
           }
         double K1=K;//uomega=-1, it is as well case of associated circus
          if (uomega != -1.0) K1=exp(log(K)* -uomega);//K is always positive
         cN[i]=K1 *ciinj;
        }          
}//END