积分图

引言

整理。

1.原

对于一幅灰度图像，积分图像中的任意一点(x,y)的值是指从图像的左上角到这个点的所构成的矩形区域内所有的点的灰度值之和：I表示积分图像，G表示原始图像。则 I(x,y)=sum(G(i,j))，其中 0<=i<=x,0<=j<=y.

2.改进版本1

Integral(i,j) = Integral(i,j-1) + Integral(i-1,j) - Integral(i-1,j-1) + Image(i,j);

code_1

void integral(unsigned char* inputMatrix, unsigned long* outputMatrix, int width, int height){
	// calculate integral of the first line
	for(int i=0;i<width;i++){
		outputMatrix[i] = inputMatrix[i];
		if(i>0){
			outputMatrix[i] += outputMatrix[i-1];
		}
	}
	for (int i=1;i<height;i++){
		int offset = i*width;
		// first column of each line
		outputMatrix[offset] = outputMatrix[offset-width]+inputMatrix[offset];
	    // other columns 
			outputMatrix[offset+j] = outputMatrix[offset+j-1] + outputMatrix[offset-width] + outputMatrix[offset-width-1] + inputMatrix[offset]; 
		}
	}
	return ;
}

3.改进版本2

Integral(i,j) = Integral(i,j-1) + ColumnSum(j);

code_2

void fastIntegral(unsigned char* inputMatrix, unsigned long* outputMatrix, int width, int height){
	unsigned long *columnSum = new unsigned long[width]; // sum of each column
	// calculate integral of the first line
	for(int i=0;i<width;i++){
		columnSum[i]=inputMatrix[i];
		outputMatrix[i] = inputMatrix[i];
		if(i>0){
			outputMatrix[i] += outputMatrix[i-1];
		}
	}
	for (int i=1;i<height;i++){
		int offset = i*width;
		// first column of each line
		columnSum[0] +=inputMatrix[offset];
		outputMatrix[offset] = columnSum[0];
		 // other columns 
		for(int j=1;j<width;j++){
			columnSum[j] += inputMatrix[offset+j];
			outputMatrix[offset+j] = outputMatrix[offset+j-1] + columnSum[j]; 
		}
	}
	return ;
}

4.改进版本3

Integral(i,j) = Integral(i-1,j) + ColumnSum(i);

（图片标注有误，应该是 Integral(i-1,j)）。

code_3

void GetGrayIntegralImage(unsigned char *Src, int *Integral, int Width, int Height, int Stride)
{
    memset(Integral, 0, (Width + 1) * sizeof(int));                    //    第一行都为0
    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        int *LinePL = Integral + Y * (Width + 1) + 1;                　　//    上一行位置            
        int *LinePD = Integral + (Y + 1) * (Width + 1) + 1;            //    当前位置，注意每行的第一列的值都为0
        LinePD[-1] = 0;                                                //    第一列的值为0
        for (int X = 0, Sum = 0; X < Width; X++)
        {
            Sum += LinePS[X];                                        //    行方向累加
            LinePD[X] = LinePL[X] + Sum;                            //    更新积分图
        }
    }
}

code_4

void GetGrayIntegralImage(unsigned char *Src, int *Integral, int Width, int Height, int Stride)
{
    int *ColSum = (int *)calloc(Width, sizeof(int));        //    用的calloc函数哦，自动内存清0
    memset(Integral, 0, (Width + 1) * sizeof(int));
    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        int *LinePL = Integral + Y * (Width + 1) + 1;
        int *LinePD = Integral + (Y + 1) * (Width + 1) + 1;
        LinePD[-1] = 0;
        for (int X = 0; X < Width; X++)
        {
            ColSum[X] += LinePS[X];
            LinePD[X] = LinePD[X - 1] + ColSum[X];
        }
    }
    free(ColSum);
}

5.应用

均值滤波

void BoxBlur(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Radius)
{
    int *Integral = (int *)malloc((Width + 1) * (Height + 1) * sizeof(int));
    GetGrayIntegralImage(Src, Integral, Width, Height, Stride);
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        int Y1 = max(Y - Radius, 0);
        int Y2 = min(Y + Radius + 1, Height - 1);
        /*int Y1 = Y - Radius;
        int Y2 = Y + Radius + 1;
        if (Y1 < 0) Y1 = 0;
        if (Y2 > Height) Y2 = Height;*/
        int *LineP1 = Integral + Y1 * (Width + 1);
        int *LineP2 = Integral + Y2 * (Width + 1);
        unsigned char *LinePD = Dest + Y * Stride;
        for (int X = 0; X < Width; X++)
        {
            int X1 = max(X - Radius, 0);
            int X2 = min(X + Radius + 1, Width);
            //    int X1 = X - Radius;
            //    if (X1 < 0) X1 = 0;
            //    int X2 = X + Radius + 1;
            //    if (X2 >= Width) X2 = Width - 1;
            int Sum = LineP2[X2] - LineP1[X2] - LineP2[X1] + LineP1[X1];
            int PixelCount = (X2 - X1) * (Y2 - Y1);
            LinePD[X] = (Sum + (PixelCount >> 1)) / PixelCount;
        }
    }
    free(Integral);
}

均值2

利用积分图进行均值滤波 



void integral_blur()
{
    Mat src = imread("src960_720.bmp", 1);
    int height = src.rows;
    int width = src.cols;
    int Rmax = 5;
    int sh2 = Rmax / 2;
    int sw2 = Rmax / 2;
    int Rsize = Rmax;
    int RsizeArea = 1.0 / (Rsize * Rsize) * 2048;
    cout << RsizeArea << endl;
    Mat srcborder;
    MakeBorder(src, srcborder, sw2, sh2);
    int ybegin = sh2;
    //int yend = height - sh2;
    int xbegin = sw2;
    //int xend = width - sw2;
    int nchannel = src.channels();
    Mat blurimg(height, width, CV_8UC3);
    for (int i = 0; i < height; ++i)
    {
        uchar *blurimgRow = blurimg.ptr<uchar>(i);
        for (int j = 0; j < width; ++j)
        {
            uchar *blurimgCol = blurimgRow + j * nchannel;
            int nCenter = Rsize >> 1;
            int sumR = 0;
            int sumG = 0;
            int sumB = 0;
            int jxbegin = j + xbegin - nCenter;
            int iybegin = i + ybegin - nCenter;
            for (int y = 0; y < Rsize; ++y)
            {
                uchar *srcRow = srcborder.ptr<uchar>(iybegin + y);
                for (int x = 0; x < Rsize; ++x)
                {
                    uchar *srcCol = srcRow + (jxbegin + x) *nchannel;
                    sumR += srcCol[2];
                    sumG += srcCol[1];
                    sumB += srcCol[0];
                }
            }
            blurimgCol[0] = (sumB * RsizeArea) >> 11;
            blurimgCol[1] = (sumG * RsizeArea) >> 11;
            blurimgCol[2] = (sumR * RsizeArea) >> 11;
        }
    }

 


    Mat srcborder_integral;
    integral(srcborder, srcborder_integral, CV_32S);
    Mat blurimg_integral(height, width, CV_8UC3);
    for (int i = 0; i < height; ++i)
    {
        uchar *blurimgRow = blurimg_integral.ptr<uchar>(i);
        for (int j = 0; j < width; ++j)
        {
            uchar *blurimgCol = blurimgRow + j * nchannel;
            int nCenter = Rsize >> 1;
            int jxbegin = j + xbegin - nCenter;
            int iybegin = i + ybegin - nCenter;
            int jxend = jxbegin + Rsize;
            int iyend = iybegin + Rsize;

            int *p1 = srcborder_integral.ptr<int>(iybegin) +(jxbegin)*nchannel;
            int *p2 = srcborder_integral.ptr<int>(iybegin) +(jxend)*nchannel;
            int *p3 = srcborder_integral.ptr<int>(iyend) +(jxbegin)*nchannel;
            int *p4 = srcborder_integral.ptr<int>(iyend) +(jxend)*nchannel;

            int sumR = p1[2] +  p4[2] -  p2[2] - p3[2];
            int sumG = p1[1] + p4[1] - p2[1] - p3[1];
            int sumB = p1[0] + p4[0] - p2[0] - p3[0];
            blurimgCol[0] = (sumB * RsizeArea) >> 11;
            blurimgCol[1] = (sumG * RsizeArea) >> 11;
            blurimgCol[2] = (sumR * RsizeArea) >> 11;
        }
    }


    Mat blurimg_opencv;
    blur(src, blurimg_opencv, Size(Rmax, Rmax));

    namedWindow("src", 0);
    imshow("src", src);
    namedWindow("dst0", 0);
    imshow("dst0", blurimg);
    namedWindow("dst1", 0);
    imshow("dst1", blurimg_opencv);
    namedWindow("dst2", 0);
    imshow("dst2", blurimg_integral);
    waitKey(0);
}
//https://www.cnblogs.com/adong7639/p/5832327.html

参考1：快速计算积分图
参考2：快速高效的积分图像算法
参考3
参考4：opencv积分图计算
参考5：基于积分图的均值滤波算法
参考6：积分图应用–局部标准差
参考7：积分图应用–非局部均值去噪（NL-means）
参考8：OpenCV中积分图介绍与应用.