[구현(c)/2008] Efficient Implementation of Local Adaptive Thresholding Techniques Using Integral Images

link: https://staffhome.ecm.uwa.edu.au/~00082689/papers/Shafait-efficient-binarization-SPIE08.pdf

 

#define OBJECT						255
#define BACKGROUND					0

#define SAUVOLA_WINDOW_SIZE			31
#define SAUVOLA_K_VAL				0.85f
#define SAUVOLA_R_VAL				128

#define FIND_MAX(a,b)						(((a) > (b)) ? (a) : (b))
#define FIND_MIN(a,b)						(((a) < (b)) ? (a) : (b))

#define MAKE_INTEGRAL_IMAGE(x, y)			(integralImage[(y)*imageWidth+(x)]) 
#define MAKE_INTEGRAL_SQUARE_IMAGE(x, y)	(integralSquareImage[(y)*imageWidth+(x)])  

void RunIntegralSauvola(unsigned char* srcImage, unsigned char* dstImage, const int imageWidth, const int imageHeight)
{
	int halfWindow = SAUVOLA_WINDOW_SIZE>>1; 
	int differece, squareDifference;  

	std::vector<int> integralImage(imageWidth*imageHeight);  
	std::vector<int> integralSquareImage(imageWidth*imageHeight);  

	integralImage[0] = srcImage[0];
 
	for(int x=1; x<imageWidth; ++x) 
	{
		int pixelVal = srcImage[x];   

		integralImage[x] = integralImage[x-1]+pixelVal;   
		integralSquareImage[x] = integralSquareImage[x-1]+pixelVal*pixelVal; 

	}

	for(int y=1; y<imageHeight; ++y) 
	{
		int lineSum =0, lineSquareSum = 0; 

		for(int x=0; x<imageWidth; ++x)  
		{
			int pixelVal = srcImage[y*imageWidth+x]; 

			lineSum += pixelVal;  
			lineSquareSum += pixelVal*pixelVal; 

			__int64	a = integralSquareImage[(y-1)*imageWidth+x];
			__int64	b = lineSquareSum;

			integralImage[y*imageWidth+x] = integralImage[(y-1)*imageWidth+x] + lineSum; 
			integralSquareImage[y*imageWidth+x] = integralSquareImage[(y-1)*imageWidth+x] + lineSquareSum; 
		}

	}

	for(int j=0; j<imageHeight; j++) 
	{
		for(int i=0; i<imageWidth; i++) 
		{ 
			int xMin = FIND_MAX(0, i-halfWindow);
			int yMin = FIND_MAX(0, j-halfWindow);
			int xMax = FIND_MIN(imageWidth-1, i+halfWindow);
			int yMax = FIND_MIN(imageHeight-1, j+halfWindow);
            
			int area = (xMax-xMin+1)*(yMax-yMin+1);

			if(!xMin && !yMin)
			{
				differece = MAKE_INTEGRAL_IMAGE(xMax, yMax);
				squareDifference = MAKE_INTEGRAL_SQUARE_IMAGE(xMax, yMax);
			}
            
			else if(!xMin && yMin)
			{
				differece = MAKE_INTEGRAL_IMAGE(xMax, yMax) - MAKE_INTEGRAL_IMAGE(xMax, yMin-1);
				squareDifference = MAKE_INTEGRAL_SQUARE_IMAGE(xMax, yMax) - MAKE_INTEGRAL_SQUARE_IMAGE(xMax, yMin-1);
			}

			else if(xMin && !yMin)
			{
				differece = MAKE_INTEGRAL_IMAGE(xMax, yMax) - MAKE_INTEGRAL_IMAGE(xMin-1, yMax);
				squareDifference = MAKE_INTEGRAL_SQUARE_IMAGE(xMax, yMax) - MAKE_INTEGRAL_SQUARE_IMAGE(xMin-1, yMax);

			}

			else 
			{
				int diagonalSum1 = MAKE_INTEGRAL_IMAGE(xMax, yMax) + MAKE_INTEGRAL_IMAGE(xMin-1, yMin-1);
				int diagonalSum2 = MAKE_INTEGRAL_IMAGE(xMax, yMin-1) + MAKE_INTEGRAL_IMAGE(xMin-1, yMax);

				differece = diagonalSum1 - diagonalSum2;
                
				int squareDiagonalSum1 = MAKE_INTEGRAL_SQUARE_IMAGE(xMax, yMax) + MAKE_INTEGRAL_SQUARE_IMAGE(xMin-1, yMin-1);
				int squareiDiagonalSum2 = MAKE_INTEGRAL_SQUARE_IMAGE(xMax, yMin-1) + MAKE_INTEGRAL_SQUARE_IMAGE(xMin-1, yMax);
				
                squareDifference = squareDiagonalSum1 - squareiDiagonalSum2;

			}

			float mean = (float)(differece)/(float)(area); 

			float	a = ( squareDifference-(float)(differece)*(float)(differece) / (float)(area) ) / (float)(area-1);
			float	b_mean = differece / (float)area;
			float	b = squareDifference / (float)area - b_mean * b_mean;

			if( abs( a - b ) > 0.00001 )
			{
				int a = 0;
			}

			float stdDev = sqrt( ( squareDifference-(float)(differece)*(float)(differece) / (float)(area) ) / (float)(area-1) );

			int findThreshold = (int)(( mean * ( 1.0f + SAUVOLA_K_VAL*( (stdDev/SAUVOLA_R_VAL) - 1.0f) ) )); 
            
			dstImage[j*imageWidth+i] = ( ( srcImage[j*imageWidth+i] > findThreshold ) ? OBJECT : BACKGROUND );
			
		}

	}

}