home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 5 / Amiga Tools 5.iso / grafik / converter / limbo4.0 / src / 2d / classify.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-07-16  |  5.2 KB  |  185 lines
  1. #include "includes.h"
  2.  
  3.  unsigned long *ML;
  4.  long double *M;
  5.  long double *CM;
  6.  long double *NCM;
  7.  long double *n;
  8.  
  9.  long double Pow2(long double d)
  10.   {
  11.    return (long double)(d*d);
  12.   }
  13.  
  14. /**************************************|****************************************
  15. Routine   : Classify
  16. Input  par: BitMap *SubImage (pointer to structure which determines image to classify)
  17. Output par: Int (Integer which determines the class of 'SubImage')
  18. Function  : Determines the class of a given subimage. Possible classes: EDGE, SHADE.
  19. ***************************************|***************************************/
  20.  
  21.  float *Classify(BitMap *Image,int XStart,int YStart, int SubSize,
  22.                    int Dim, float TSE)
  23.   {
  24.    register unsigned int x,y;
  25.    unsigned long f,ff=0,xx,yy,xy,xxy,xyy,m=0;
  26.    long double mx,my,iM2,iM3,temp;
  27.    static int FirstCall=TRUE;
  28.    float *Features=GimmeAFloatArray(Dim+PREDEFFEATURES);
  29.    double Volume=SubSize*SubSize;
  30.    
  31.    if(FirstCall)
  32.     {
  33.      FirstCall=FALSE;
  34.      ML =GimmeAULongArray(10);
  35.      M  =GimmeADoubleArray(10);
  36.      CM =GimmeADoubleArray(10);
  37.      NCM=GimmeADoubleArray(7);
  38.      n  =GimmeADoubleArray(7);
  39.     }
  40.    
  41.    for(x=0;x<10;x++) ML[x]=0;
  42.    
  43.    for (x=0;x<SubSize;x++)
  44.    for (y=0;y<SubSize;y++)
  45.     {
  46.       f=Image->Map[x+XStart][y+YStart];
  47.       m += f;
  48.       ff += f*f;
  49.     }
  50.    
  51.    Features[MEAN]=(float)m/Volume;
  52.    Features[VAR] =((float)ff-Volume*Features[MEAN]*Features[MEAN])/Volume;
  53.    
  54.    if (Features[VAR]>TSE) /* is this an edge ? */
  55.     {
  56.      for (x=0;x<SubSize;x++)
  57.      for (y=0;y<SubSize;y++)
  58.       {
  59.        f=Image->Map[x+XStart][y+YStart];
  60.        xx=(x+1)*(x+1);
  61.        yy=(y+1)*(y+1);
  62.        xy=(x+1)*(y+1);
  63.        xxy=xx*(y+1);
  64.        xyy=(x+1)*yy;
  65.        ML[0] += f;          /* m_00 */
  66.        ML[1] += (x+1)*f;    /* m_10 */
  67.        ML[2] += (y+1)*f;    /* m_01 */
  68.        ML[3] += xy*f;       /* m_11 */
  69.        ML[4] += xx*f;       /* m_20 */
  70.        ML[5] += yy*f;       /* m_02 */
  71.        ML[6] += xxy*f;      /* m_21 */
  72.        ML[7] += xyy*f;      /* m_12 */
  73.        ML[8] += (x+1)*xx*f; /* m_30 */
  74.        ML[9] += (y+1)*yy*f; /* m_03 */
  75.       }
  76.      
  77.      Features[STDDEV]=sqrt(Features[VAR]);
  78.  
  79.      for(x=0;x<10;x++) M[x]=(double)(ML[x]);
  80.      
  81.      if(M[0]==0.0)
  82.       {
  83.        ErrorHandler(OUT_OF_RANGE2,"Null mean in Classify");
  84.        M[0]=1E-100;
  85.       }
  86.      
  87.      mx=M[1]/M[0];
  88.      my=M[2]/M[0];
  89.      
  90.      CM[0] = M[0];              /* u_00 */
  91.      CM[1] = 0.0;               /* u_10 */
  92.      CM[2] = 0.0;               /* u_01 */
  93.      CM[3] = M[3]-my*M[1];      /* u_11 */
  94.      CM[4] = M[4]-mx*M[1];      /* u_20 */
  95.      CM[5] = M[5]-my*M[2];      /* u_02 */
  96.      CM[6] = M[6]-2.0*mx*M[3]-my*M[4]+2.0*mx*mx*M[2]; /* u_21 */
  97.      CM[7] = M[7]-2.0*my*M[3]-mx*M[5]+2.0*my*my*M[1]; /* u_12 */
  98.      CM[8] = M[8]-3.0*mx*M[4]+2.0*mx*mx*M[1];         /* u_30 */
  99.      CM[9] = M[9]-3.0*my*M[5]+2.0*my*my*M[2];         /* u_03 */
  100.      
  101.      iM2=1.0/(M[0]*M[0]);
  102.      iM3=1.0/pow(M[0],2.5);
  103.      
  104.      n[0]=CM[3]*iM2; /* n_11 */
  105.      n[1]=CM[4]*iM2; /* n_20 */
  106.      n[2]=CM[5]*iM2; /* n_02 */
  107.      n[3]=CM[6]*iM3; /* n_21 */
  108.      n[4]=CM[7]*iM3; /* n_12 */
  109.      n[5]=CM[8]*iM3; /* n_30 */
  110.      n[6]=CM[9]*iM3; /* n_03 */
  111.  
  112.      /* These features may give rise to numerical problems
  113.         with floats which for some systems must be in 
  114.         the range [1E-37;1E37]  */
  115.  
  116.      NCM[0]=n[1]+n[2];                               /* p_1 */
  117.      if (Dim>1)  
  118.       {
  119.        NCM[1]=Pow2(n[1]-n[2])+                       /* p_2 */
  120.        4.0*Pow2(n[0]);
  121.        if (Dim>2) 
  122.         {
  123.          NCM[2]=Pow2(n[5]-3.0*n[4])+                 /* p_3 */
  124.          Pow2(3.0*n[3]-n[6]);
  125.          if (Dim>3) 
  126.           {
  127.            NCM[3]=Pow2(n[5]+n[4])+                   /* p_4 */
  128.            Pow2(n[3]+n[6]);
  129.            if (Dim>4) 
  130.             {
  131.              NCM[4]=(n[5]-3.0*n[4])*(n[5]+n[4])*     /* p_5 */
  132.              (Pow2(n[5]+n[4])-
  133.               3.0*Pow2(n[3]+n[6]))+
  134.              (3.0*n[3]-n[6])*(n[3]+n[6])*
  135.              (3.0*Pow2(n[5]+n[4])-
  136.               Pow2(n[3]+n[6]));
  137.              if (Dim>5) 
  138.               {
  139.                NCM[5]=(n[1]-n[2])*                   /* p_6 */
  140.                (Pow2(n[5]+n[4])-      
  141.                 Pow2(n[3]+n[6]))+
  142.                4.0*n[0]*(n[5]+n[4])*
  143.                (n[3]+n[6]);
  144.                if (Dim>6) 
  145.                 {
  146.                  NCM[6]=(3.0*n[3]-n[5])*(n[5]+n[4])* /* p_7 */
  147.                  (Pow2(n[5]+n[4])-      
  148.                   3.0*Pow2(n[3]+n[6]))+
  149.                  (3.0*n[4]-n[5])*(n[3]+n[6])*
  150.                  (3.0*Pow2(n[5]+n[4])-
  151.                   Pow2(n[3]+n[6]));
  152.                 }
  153.               }
  154.             }
  155.           }
  156.         }
  157.       }
  158.    
  159.      for(x=0;x<Dim;x++)
  160.       {
  161.        temp=fabs(NCM[x]);
  162.        if (NCM[x]>=0.0) Features[x+PREDEFFEATURES]=(float)log(temp);
  163.        else             Features[x+PREDEFFEATURES]=(float)-log(temp);
  164.  
  165.        if (temp==0.0) 
  166.         {
  167.          vprintf(stderr,"Feature: %d",x); 
  168.          ErrorHandler(OUT_OF_RANGE2,"Feature value out of range");
  169.          Features[VAR]=-99999; /* replace with shade */
  170.         }
  171.       }
  172.      
  173.     }
  174.    else 
  175.     {
  176.      for(x=0;x<Dim;x++) Features[x+PREDEFFEATURES]=0.0;
  177.      Features[VAR] = -Features[VAR]; /* negative variance indicates a shade */
  178.     }
  179.    
  180.    return Features;
  181.   }
  182.  
  183.  
  184.  
  185.