【短道速滑十】从单幅图像中评估加性噪音的均方差。

　　在Halcon中有这样一个函数：

　　estimate_noise estimate_noise — Estimate the image noise from a single image.

　　Signature

　　　　estimate_noise(Image : : Method, Percent : Sigma)

　　Description

　　The operator estimate_noise estimates the standard deviation of additive noise within the domain of the image that is passed in Image. The standard deviation is returned in Sigma.

　　　　The operator is useful in the following use cases:

　　　　determination of MinContrast for matching,

　　　　determination of the amplitude for edge filters,

　　　　camera evaluation,

　　　　monitoring errors in camera operation (e.g., user overdrives camera gain).

　　即从单幅图像中评估图像噪音的均方差，这个算子可以用于计算匹配时的最小对比度（发现新大陆了，原路模板匹配还可以用这个做自动化）、边缘检测滤波器的幅度、摄像机评估、控相机操作中的错误（例如用户过度调节相机增益）。

       我觉得还可以把他作为自动去噪的一个参考指标。

　　Halcon里提供了四个评估噪音的方法：: 'foerstner', 'immerkaer', 'least_squares', 'mean'，其本身最推荐的方法是immerkaer，如其帮助文档里所说：

　　Use the method 'immerkaer', instead of the methods 'foerstner', 'least_squares', or 'mean'. The method 'immerkaer' does not rely on the existence of homogeneous image regions, and hence is almost always applicable.

　　关于immerkaer方法，开放的Halcon基本上提供了完整的算法思路：

　　　　'immerkaer': If Method is set to 'immerkaer', first the following filter mask is applied to the input image:  

　　　　The advantage of this method is that M is almost insensitive to image structure but only depends on the noise in the image. Assuming a Gaussian distributed noise, its standard deviation is finally obtained as

            where N is the number of image pixels to which M is applied. Note that the result obtained by this method is independent of the value passed in Percent.

　　这个M算子明显就是类似一个边缘检测的算子，然后把所有这个算子的结果相加，再求某个意义下的平均值，Halcon说这个方法的好处是对图像的结构不敏感，而只完全依赖于图像的噪音本身。 

　　我想有了这个提示，要实现这个功能应该就是很简单的过程了。 

　　我的一个实现如下所示：

//    模拟实现halcon的estimate_noise函数

int IM_EstimateNoise(unsigned char *Src, int Width, int Height, int Stride, float &Sigma)
{
    int Channel = Stride / Width;
    if (Src == NULL)                                return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                return IM_STATUS_INVALIDPARAMETER;
    if (Channel != 1)                                return IM_STATUS_NOTSUPPORTED;
    unsigned int Sum = 0;
    for (int Y = 1; Y < Height - 1; Y++)
    {
        unsigned char *LinePL = Src + (Y - 1) * Stride;
        unsigned char *LinePC = Src + (Y - 0) * Stride;
        unsigned char *LinePN = Src + (Y + 1) * Stride;
        for (int X = 1; X < Width - 1; X++)
        {
            int L = LinePL[X - 1] - 2 * LinePL[X] + LinePL[X + 1];
            int C = -2 * LinePC[X - 1] + 4 * LinePC[X] - 2 * LinePC[X + 1];
            int N = LinePN[X - 1] - 2 * LinePN[X] + LinePN[X + 1];
            Sum += IM_Abs(L + C + N);
        }
    }
    Sigma = sqrtf(IM_PI / 2) / (6 * Width * Height) * Sum;
    return IM_STATUS_OK;
}

　　为了简化代码，没有考虑图像周边单位像素的信息了，如果要严格意义的实现，也应该不是很难吧。 

　　我们比较下halcon的结果和上面这段代码的结果，使用Halcon自带的测试代码和图片：

dev_update_off ()
dev_close_window ()
dev_open_window (0, 0, 512, 512, 'black', WindowHandle)
set_display_font (WindowHandle, 14, 'mono', 'true', 'false')
TestImages := ['for6','pumpe','die/die_02','clip','ic0','horses','board/board-01','combine']
NumImages := |TestImages|
for I := 0 to NumImages - 1 by 1
    read_image (Image, TestImages[I])
    dev_resize_window_fit_image (Image, 0, 0, -1, -1)
    dev_display (Image)
    for AddedNoise := 0 to 15 by 5
        gauss_distribution (AddedNoise + 1e-2, Distribution)
        add_noise_distribution (Image, ImageNoise, Distribution)
        write_image (ImageNoise, 'bmp', 0, 'C:/Users/Administrator/Desktop/1.bmp')
        estimate_noise (ImageNoise, 'foerstner', 20, SigmaFoerstner)
        estimate_noise (ImageNoise, 'immerkaer', 20, SigmaImmerkaer)
        estimate_noise (ImageNoise, 'least_squares', 20, SigmaLeastSquares)
        estimate_noise (ImageNoise, 'mean', 20, SigmaMean)
        dev_display (ImageNoise)
        disp_message (WindowHandle, 'Added Gaussian noise: Sigma = ' + AddedNoise, 'window', 12, 12, 'black', 'true')
        Message := 'Estimated image noise (Sigma):'
        Message[1] := 'Method \'foerstner\':     ' + SigmaFoerstner$'5.2f'
        Message[2] := 'Method \'immerkaer\':     ' + SigmaImmerkaer$'5.2f'
        Message[3] := 'Method \'least_squares\': ' + SigmaLeastSquares$'5.2f'
        Message[4] := 'Method \'mean\':          ' + SigmaMean$'5.2f'
        disp_message (WindowHandle, Message, 'windowe', 40, 12, 'black', 'true')
        disp_continue_message (WindowHandle, 'black', 'true')
        stop ()
    endfor
endfor  

　　　　　　　　　　　　　　 噪音图像　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　Halcon的结果

 　使用上述C的代码获取的结果为： 5.240565，和Halcon的结果基本一致。

　我们再找一些正常的图片看看这个噪音值是否合理：

　　　　　　　　　　噪音值 0.7763　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　噪音值 2.6604

　　基本啥都比较小。 

　　　　　　　　　　　        　        噪音值  7.2155　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　   噪音值  20.04

　　对于高斯噪音，如上所示，还是能明显的区别出来的。

　　不过测试也表面，有些图的噪音虽然视觉看起来比较明显，但是用这参数去衡量时，确是很小，这个可能是因为他针对的是加性噪音做的评估吧。

　　参考资料：

　　　　W. F?rstner: “Image Preprocessing for Feature Extraction in Digital Intensity, Color and Range Images“, Springer Lecture Notes on Earth Sciences, Summer School on Data Analysis and the Statistical Foundations of Geomatics, 1999
　　　　J. Immerkaer: “Fast Noise Variance Estimation“, Computer Vision and Image Understanding, Vol. 64, No. 2, pp. 300-302, 1996