/*

-----------------------------------------------------------------------------------------------
Fast Euler Numbers thresholding
-----------------------------------------------------------------------------------------------

Computes optimal threshold according to fast Euer numbers method
("Real-time thresholding with Euler numbers", Pattern Recognition Letters).

Input:	diff is a 8 bit difference image (e.g. |Background-Current|) as OpenCV IplImage

Output:	ts is the optimal threshold in [0,255] range according to the Euler numbers techniques


Notes:
- Code uses OpenCV
- Code is not optimized for speed (e.g. cvGet2D is slow)
- Many thanks to Prof. Paul Rosin for his precious help in debugging the code


Version 1.0

Lauro Snidaro 2013


*/



#include <algorithm>

unsigned int FastEuler(const IplImage* diff)
{
    /* Takes the difference map and returns the threshold for a 8bit image*/

    unsigned int* q1start= new unsigned int[256];
    unsigned int* q1stop=  new unsigned int[256];
    unsigned int* q3start= new unsigned int[256];
    unsigned int* q3stop=  new unsigned int[256];
    unsigned int* qdstart= new unsigned int[256];
    unsigned int* qdstop=  new unsigned int[256];
    int* q1=new int[256];
    int* q3=new int[256];
    int* qd=new int[256];
    int* euler=new int[256];
    float* de=new float[256];


    for (register unsigned int i=0; i<256; i++)
    {
        q1start[i]=0;
        q1stop[i]=0;
        q3start[i]=0;
        q3stop[i]=0;
        qdstart[i]=0;
        qdstop[i]=0;
    }


	//Check if ROI has been set
	CvRect roi=cvGetImageROI(diff);

	const unsigned int rx1=roi.x;
	const unsigned int ry1=roi.y;
	const unsigned int rx2=rx1+roi.width;
	const unsigned int ry2=ry1+roi.height;
	const unsigned int width= roi.width;

    for (unsigned int y=ry1; y<ry2-1; y+=1)
    {
        for (unsigned int x=rx1; x<rx2-1; x+=1)
        {
			uchar* temp_ptr = (uchar*)(diff->imageData);
            unsigned int p1=(y*diff->widthStep)+x;
            unsigned int p2=p1+width;
            unsigned int p3=p1+1;
            unsigned int p4=p2+1;

            int min, min2, max2, max;
            int sortingVector[4];
			
			int p1_v=min=max=max2=min2=sortingVector[0]=(cvGet2D(diff, y,x)).val[0];
			int p2_v=sortingVector[1]=(cvGet2D(diff, y+1,x)).val[0];
			int p3_v=sortingVector[2]=(cvGet2D(diff, y,x+1)).val[0];
			int p4_v=sortingVector[3]=(cvGet2D(diff, y+1,x+1)).val[0];



	// SORTING
		
		   int elements = sizeof(sortingVector) / sizeof(sortingVector[0]); 
			std::sort(sortingVector, sortingVector + elements);
	
			min=sortingVector[0];
			min2=sortingVector[1];
			max2=sortingVector[2];
			max=sortingVector[3];
			

		   if(max2 <max) // determining the q1 quad contributing interval
           {
                q1start[max2]++;
				q1stop[max]++;
            }

			if (min<min2) //determining the q3 quad contributing interval

            {
                q3start[min]++;
                q3stop[min2]++;
            }

            if  ( (min2<max2) && (max2<=max) && (((max2<=p1_v && p1_v<=max) && (max2<=p4_v && p4_v<=max)) ||
                 ((max2<=p2_v && p2_v<=max) && (max2<=p3_v && p3_v<=max))))


            {   //determining the qd quad contributing interval
                qdstart[min2]++;
				qdstop[max2]++;
            }

   
        }
    }

    //constructing the euler vector using q1,q3,qd vectors

    q1[0]= q1start[0];// -q1stop[0];
    q3[0]= q3start[0];// -q3stop[0];
    qd[0]= qdstart[0];// -qdstop[0];
    euler[0] = (q1[0]- q3[0] - 2*qd[0])/4;

    de[0]=0;
    unsigned int izero=0;
    if (euler[0])  izero++;

	{
    for (register int i=1; i<256; i++)
    {
       q1[i]= q1[i-1] + (q1start[i]-q1stop[i]);
	   q3[i]= q3[i-1] + (q3start[i]-q3stop[i]);
	   qd[i]= qd[i-1] + (qdstart[i]-qdstop[i]);
	   euler[i] = (q1[i]- q3[i] - 2*qd[i])/4;

	   //Stopping condition when the euler number gets to zero
       //if (!euler[i])
       //{
       //     izero=i;
       //     break;
       //}
    }
	}


    /* finding maximum */
	unsigned int maxe = 0;
	for (register unsigned int k = 0; k < izero; k++)
    {
       if (euler[maxe] <= euler[k]) maxe = k;
    }


/*
	Calculating coefficients of the line passing through maximum and minimum of the Euler numbers curve
	(x1,y1) maxmimum
	(x2, y2) minimum
	
	Equation:
	y-y1     x-x1
	----- =  -----
	y2-y1    y2-y1

	in implicit form:

	a=1/(x2-x1)
	b= -1/(y2-y1)
	c= [y1(x2-x1)-x1(y2-y1)]/ [(y2-y1)(x2-x1)]
*/	
	
	float x1= float(maxe);
	float x2= float(izero);
	float y1= float(euler[maxe]);
	float y2= float(euler[izero]);

	float a = 1 / (x2-x1);
	float b = -1 / (y2-y1);
	float c = (y1*(x2-x1)-x1*(y2-y1))/ ((y2-y1)*(x2-x1));
		
	
	
	/* distance line-point
       d = (a*x0 + b*y0 + c) / sqrt(aČ+bČ)*/

   float denominator = sqrt(a*a + b*b);

   /* optimal threshold at maximum distance */
   unsigned int ts = maxe;
   float max = 0;
   {
   for (unsigned register int i=maxe; i<=izero; i++)
   {
     float max_t =  abs((float)(a*i + b*euler[i] + c)) / denominator;
     if (max_t > max) {ts = i; max = max_t;}
   }
   }


   // "Zero-Crossing" technique

    for (unsigned register int i = 1; i < 255; i++)
    {
       de[i]=(euler[i]-euler[i-1]);
       if (  (de[i]*de[i-1]<=0) && (i > ts))
       {
            ts=i;
            break;
       }
    }





   //Cleaning up
   delete[] q1start;
   delete[] q1stop;
   delete[] q3start;
   delete[] q3stop;
   delete[] qdstart;
   delete[] qdstop;
   delete[] q1;
   delete[] q3;
   delete[] qd;
   delete[] euler;
   delete[] de;

   return  ts;
}
//--------------------------------------------------------------------------