#include "includes.h"
 
 unsigned long *ML;
 long double *M;
 long double *CM;
 long double *NCM;
 long double *n;
 
 long double Pow2(long double d)
  {
   return (long double)(d*d);
  }
 
/**************************************|****************************************
Routine   : Classify
Input  par: BitMap *SubImage (pointer to structure which determines image to classify)
Output par: Int (Integer which determines the class of 'SubImage')
Function  : Determines the class of a given subimage. Possible classes: EDGE, SHADE.
***************************************|***************************************/
 
 float *Classify(BitMap *Image,int XStart,int YStart, int SubSize,
                   int Dim, float TSE)
  {
   register unsigned int x,y;
   unsigned long f,ff=0,xx,yy,xy,xxy,xyy,m=0;
   long double mx,my,iM2,iM3,temp;
   static int FirstCall=TRUE;
   float *Features=GimmeAFloatArray(Dim+PREDEFFEATURES);
   double Volume=SubSize*SubSize;
   
   if(FirstCall)
    {
     FirstCall=FALSE;
     ML =GimmeAULongArray(10);
     M  =GimmeADoubleArray(10);
     CM =GimmeADoubleArray(10);
     NCM=GimmeADoubleArray(7);
     n  =GimmeADoubleArray(7);
    }
   
   for(x=0;x<10;x++) ML[x]=0;
   
   for (x=0;x<SubSize;x++)
   for (y=0;y<SubSize;y++)
    {
      f=Image->Map[x+XStart][y+YStart];
      m += f;
      ff += f*f;
    }
   
   Features[MEAN]=(float)m/Volume;
   Features[VAR] =((float)ff-Volume*Features[MEAN]*Features[MEAN])/Volume;
   
   if (Features[VAR]>TSE) /* is this an edge ? */
    {
     for (x=0;x<SubSize;x++)
     for (y=0;y<SubSize;y++)
      {
       f=Image->Map[x+XStart][y+YStart];
       xx=(x+1)*(x+1);
       yy=(y+1)*(y+1);
       xy=(x+1)*(y+1);
       xxy=xx*(y+1);
       xyy=(x+1)*yy;
       ML[0] += f;          /* m_00 */
       ML[1] += (x+1)*f;    /* m_10 */
       ML[2] += (y+1)*f;    /* m_01 */
       ML[3] += xy*f;       /* m_11 */
       ML[4] += xx*f;       /* m_20 */
       ML[5] += yy*f;       /* m_02 */
       ML[6] += xxy*f;      /* m_21 */
       ML[7] += xyy*f;      /* m_12 */
       ML[8] += (x+1)*xx*f; /* m_30 */
       ML[9] += (y+1)*yy*f; /* m_03 */
      }
     
     Features[STDDEV]=sqrt(Features[VAR]);

     for(x=0;x<10;x++) M[x]=(double)(ML[x]);
     
     if(M[0]==0.0)
      {
       ErrorHandler(OUT_OF_RANGE2,"Null mean in Classify");
       M[0]=1E-100;
      }
     
     mx=M[1]/M[0];
     my=M[2]/M[0];
     
     CM[0] = M[0];              /* u_00 */
     CM[1] = 0.0;               /* u_10 */
     CM[2] = 0.0;               /* u_01 */
     CM[3] = M[3]-my*M[1];      /* u_11 */
     CM[4] = M[4]-mx*M[1];      /* u_20 */
     CM[5] = M[5]-my*M[2];      /* u_02 */
     CM[6] = M[6]-2.0*mx*M[3]-my*M[4]+2.0*mx*mx*M[2]; /* u_21 */
     CM[7] = M[7]-2.0*my*M[3]-mx*M[5]+2.0*my*my*M[1]; /* u_12 */
     CM[8] = M[8]-3.0*mx*M[4]+2.0*mx*mx*M[1];         /* u_30 */
     CM[9] = M[9]-3.0*my*M[5]+2.0*my*my*M[2];         /* u_03 */
     
     iM2=1.0/(M[0]*M[0]);
     iM3=1.0/pow(M[0],2.5);
     
     n[0]=CM[3]*iM2; /* n_11 */
     n[1]=CM[4]*iM2; /* n_20 */
     n[2]=CM[5]*iM2; /* n_02 */
     n[3]=CM[6]*iM3; /* n_21 */
     n[4]=CM[7]*iM3; /* n_12 */
     n[5]=CM[8]*iM3; /* n_30 */
     n[6]=CM[9]*iM3; /* n_03 */

     /* These features may give rise to numerical problems
        with floats which for some systems must be in 
        the range [1E-37;1E37]  */

     NCM[0]=n[1]+n[2];                               /* p_1 */
     if (Dim>1)  
      {
       NCM[1]=Pow2(n[1]-n[2])+                       /* p_2 */
       4.0*Pow2(n[0]);
       if (Dim>2) 
        {
         NCM[2]=Pow2(n[5]-3.0*n[4])+                 /* p_3 */
         Pow2(3.0*n[3]-n[6]);
         if (Dim>3) 
          {
           NCM[3]=Pow2(n[5]+n[4])+                   /* p_4 */
           Pow2(n[3]+n[6]);
           if (Dim>4) 
            {
             NCM[4]=(n[5]-3.0*n[4])*(n[5]+n[4])*     /* p_5 */
             (Pow2(n[5]+n[4])-
              3.0*Pow2(n[3]+n[6]))+
             (3.0*n[3]-n[6])*(n[3]+n[6])*
             (3.0*Pow2(n[5]+n[4])-
              Pow2(n[3]+n[6]));
             if (Dim>5) 
              {
               NCM[5]=(n[1]-n[2])*                   /* p_6 */
               (Pow2(n[5]+n[4])-      
                Pow2(n[3]+n[6]))+
               4.0*n[0]*(n[5]+n[4])*
               (n[3]+n[6]);
               if (Dim>6) 
                {
                 NCM[6]=(3.0*n[3]-n[5])*(n[5]+n[4])* /* p_7 */
                 (Pow2(n[5]+n[4])-      
                  3.0*Pow2(n[3]+n[6]))+
                 (3.0*n[4]-n[5])*(n[3]+n[6])*
                 (3.0*Pow2(n[5]+n[4])-
                  Pow2(n[3]+n[6]));
                }
              }
            }
          }
        }
      }
   
     for(x=0;x<Dim;x++)
      {
       temp=fabs(NCM[x]);
       if (NCM[x]>=0.0) Features[x+PREDEFFEATURES]=(float)log(temp);
       else             Features[x+PREDEFFEATURES]=(float)-log(temp);

       if (temp==0.0) 
        {
         vprintf(stderr,"Feature: %d",x); 
         ErrorHandler(OUT_OF_RANGE2,"Feature value out of range");
         Features[VAR]=-99999; /* replace with shade */
        }
      }
     
    }
   else 
    {
     for(x=0;x<Dim;x++) Features[x+PREDEFFEATURES]=0.0;
     Features[VAR] = -Features[VAR]; /* negative variance indicates a shade */
    }
   
   return Features;
  }