OpenCV GaussianBlur 结果不一致

TL;DR

OpenCV高版本的GaussianBlur，某些条件下会用定点化计算替代浮点计算，导致跨版本的结果不一致：

void GaussianBlurFixedPoint(const Mat& src, /*const*/ Mat& dst,
                            const uint16_t/*ufixedpoint16*/* fkx, int fkx_size,
                            const uint16_t/*ufixedpoint16*/* fky, int fky_size,
                            int borderType)
{
    CV_INSTRUMENT_REGION();

    CV_Assert(src.depth() == CV_8U && ((borderType & BORDER_ISOLATED) || !src.isSubmatrix()));
    fixedSmoothInvoker<uint8_t, ufixedpoint16> invoker(
            src.ptr<uint8_t>(), src.step1(),
            dst.ptr<uint8_t>(), dst.step1(), dst.cols, dst.rows, dst.channels(),
            (const ufixedpoint16*)fkx, fkx_size, (const ufixedpoint16*)fky, fky_size,
            borderType & ~BORDER_ISOLATED);
    {
        // TODO AVX guard (external call)
        parallel_for_(Range(0, dst.rows), invoker, std::max(1, std::min(getNumThreads(), getNumberOfCPUs())));
    }
}

问题描述

在移植OpenCV的GaussianBlur相关函数，测试阶段发现和OpenCV结果不一致，但是肉眼看不出差异。于是对比了多个版本的OpenCV，发现它们之间结果也不一致。

样例代码

namespace cv
{
	static void test_GaussianBlur()
	{
		std::string im_pth = "../../imgs/IU.bmp";
		Mat src = imread(im_pth);
		Mat dst;

		Size size(3, 3);
		GaussianBlur(src, dst, size, 0, 0);
		imwrite("IU-gaussian-blur.bmp", dst);
	}
}

测试结果

第一种结果： OpenCV 2.4.13 OpenCV 3.1.0

第二种结果： OpenCV 3.4.5 OpenCV 3.4.9 以及我的移植版本

说明：所用OpenCV都是PC CPU版本，关闭IPP优化。

调试排查

以OpenCV3.4.9为例，会把满足一定条件的图像，执行定点化版本的高斯模糊，而不是浮点数版本的计算。这是相当于OpenCV3.1.0 / 2.4.13版本增加的内容。

看完整源码：

void GaussianBlur(InputArray _src, OutputArray _dst, Size ksize,
                  double sigma1, double sigma2,
                  int borderType)
{
    CV_INSTRUMENT_REGION();

    int type = _src.type();
    Size size = _src.size();
    _dst.create( size, type );

    if( (borderType & ~BORDER_ISOLATED) != BORDER_CONSTANT &&
        ((borderType & BORDER_ISOLATED) != 0 || !_src.getMat().isSubmatrix()) )
    {
        if( size.height == 1 )
            ksize.height = 1;
        if( size.width == 1 )
            ksize.width = 1;
    }

    if( ksize.width == 1 && ksize.height == 1 )
    {
        _src.copyTo(_dst);
        return;
    }

    bool useOpenCL = (ocl::isOpenCLActivated() && _dst.isUMat() && _src.dims() <= 2 &&
               ((ksize.width == 3 && ksize.height == 3) ||
               (ksize.width == 5 && ksize.height == 5)) &&
               _src.rows() > ksize.height && _src.cols() > ksize.width);
    CV_UNUSED(useOpenCL);

    int sdepth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);

    Mat kx, ky;
    createGaussianKernels(kx, ky, type, ksize, sigma1, sigma2);

    CV_OCL_RUN(useOpenCL, ocl_GaussianBlur_8UC1(_src, _dst, ksize, CV_MAT_DEPTH(type), kx, ky, borderType));

    CV_OCL_RUN(_dst.isUMat() && _src.dims() <= 2 && (size_t)_src.rows() > kx.total() && (size_t)_src.cols() > kx.total(),
               ocl_sepFilter2D(_src, _dst, sdepth, kx, ky, Point(-1, -1), 0, borderType))

    Mat src = _src.getMat();
    Mat dst = _dst.getMat();

    Point ofs;
    Size wsz(src.cols, src.rows);
    if(!(borderType & BORDER_ISOLATED))
        src.locateROI( wsz, ofs );

    CALL_HAL(gaussianBlur, cv_hal_gaussianBlur, src.ptr(), src.step, dst.ptr(), dst.step, src.cols, src.rows, sdepth, cn,
             ofs.x, ofs.y, wsz.width - src.cols - ofs.x, wsz.height - src.rows - ofs.y, ksize.width, ksize.height,
             sigma1, sigma2, borderType&~BORDER_ISOLATED);

    CV_OVX_RUN(true,
               openvx_gaussianBlur(src, dst, ksize, sigma1, sigma2, borderType))

    //CV_IPP_RUN_FAST(ipp_GaussianBlur(src, dst, ksize, sigma1, sigma2, borderType));

    if(sdepth == CV_8U && ((borderType & BORDER_ISOLATED) || !_src.getMat().isSubmatrix()))
    {
        std::vector<ufixedpoint16> fkx, fky;
        createGaussianKernels(fkx, fky, type, ksize, sigma1, sigma2);

        static bool param_check_gaussian_blur_bitexact_kernels = utils::getConfigurationParameterBool("OPENCV_GAUSSIANBLUR_CHECK_BITEXACT_KERNELS", false);
        if (param_check_gaussian_blur_bitexact_kernels && !validateGaussianBlurKernel(fkx))
        {
            CV_LOG_INFO(NULL, "GaussianBlur: bit-exact fx kernel can't be applied: ksize=" << ksize << " sigma=" << Size2d(sigma1, sigma2));
        }
        else if (param_check_gaussian_blur_bitexact_kernels && !validateGaussianBlurKernel(fky))
        {
            CV_LOG_INFO(NULL, "GaussianBlur: bit-exact fy kernel can't be applied: ksize=" << ksize << " sigma=" << Size2d(sigma1, sigma2));
        }
        else
        {
            if (src.data == dst.data)
                src = src.clone();
				
			//-------------!! 注意这里，dispatch到fixedpoint这一版本的实现上
            CV_CPU_DISPATCH(GaussianBlurFixedPoint, (src, dst, (const uint16_t*)&fkx[0], (int)fkx.size(), (const uint16_t*)&fky[0], (int)fky.size(), borderType),
                CV_CPU_DISPATCH_MODES_ALL);
            return;
        }
    }

	//-------！！先前几行的dispatch分支算好后直接return，不会fall back到sepFilter2D
    sepFilter2D(src, dst, sdepth, kx, ky, Point(-1, -1), 0, borderType);
}

而OpenCV 3.1.0的实现则是这样的：

void cv::GaussianBlur( InputArray _src, OutputArray _dst, Size ksize,
                   double sigma1, double sigma2,
                   int borderType )
{
    int type = _src.type();
    Size size = _src.size();
    _dst.create( size, type );

    if( borderType != BORDER_CONSTANT && (borderType & BORDER_ISOLATED) != 0 )
    {
        if( size.height == 1 )
            ksize.height = 1;
        if( size.width == 1 )
            ksize.width = 1;
    }

    if( ksize.width == 1 && ksize.height == 1 )
    {
        _src.copyTo(_dst);
        return;
    }

#ifdef HAVE_TEGRA_OPTIMIZATION
    Mat src = _src.getMat();
    Mat dst = _dst.getMat();
    if(sigma1 == 0 && sigma2 == 0 && tegra::useTegra() && tegra::gaussian(src, dst, ksize, borderType))
        return;
#endif

    CV_IPP_RUN(true, ipp_GaussianBlur( _src,  _dst,  ksize, sigma1,  sigma2, borderType));

    Mat kx, ky;
    createGaussianKernels(kx, ky, type, ksize, sigma1, sigma2);
    sepFilter2D(_src, _dst, CV_MAT_DEPTH(type), kx, ky, Point(-1,-1), 0, borderType );
}

显然它是始终调用sepFilter2D的实现的。