图像降噪------LLSURE

代码如下：

void Denoise_LLSURE(cv::Mat &image, int radius, float sigma) {
    int width = image.cols, height = image.rows;
    cv::Mat I1, I2;
    cv::integral(image, I1, I2, CV_64FC1, CV_64FC1);
    I1 = I1(cv::Rect(1, 1, width, height));
    I2 = I2(cv::Rect(1, 1, width, height));
    cv::Mat varInv = cv::Mat::zeros(height, width, CV_64FC1);
    cv::Mat AMat = cv::Mat::zeros(height, width, CV_64FC1);
    cv::Mat BMat = cv::Mat::zeros(height, width, CV_64FC1);
    double epsilon = 1e-6f, sigma2 = sigma * sigma;
    for (int i = 0; i < height; ++i) {
        auto aData = AMat.ptr<double>(i);
        auto bData = BMat.ptr<double>(i);
        auto varDataInv = varInv.ptr<double>(i);
        int yStart = std::max(i - radius, 0);
        int yEnd = std::min(i + radius, height - 1);
        auto IDataStart1 = I1.ptr<double>(yStart);
        auto IDataStart2 = I2.ptr<double>(yStart);
        auto IDataEnd1 = I1.ptr<double>(yEnd);
        auto IDataEnd2 = I2.ptr<double>(yEnd);
        int yLen = yEnd - yStart + 1;
        for (int j = 0; j < width; ++j) {
            int xStart = std::max(j - radius, 0);
            int xEnd = std::min(j + radius, width - 1);
            int xLen = xEnd - xStart + 1;
            double invNum = 1.0f / (yLen * xLen);
            double mean = (IDataEnd1[xEnd] - IDataEnd1[xStart] + IDataStart1[xStart] - IDataStart1[xEnd]) * invNum;
            double mean2 = (IDataEnd2[xEnd] - IDataEnd2[xStart] + IDataStart2[xStart] - IDataStart2[xEnd]) * invNum;
            double var = mean2 - mean * mean + epsilon;
            double invVar = 1.0 / (var);
            varDataInv[j] = invVar;
            double a = std::max(var - sigma2, 0.0) * invVar;
            double b = (1.0 - a) * mean;
            aData[j] = invVar * a;
            bData[j] = invVar * b;
        }
    }

    cv::Mat AIntegral, BIntegral, varInv_integral;
    cv::integral(AMat, AIntegral, CV_64FC1);
    AIntegral = AIntegral(cv::Rect(1, 1, width, height));
    cv::integral(BMat, BIntegral, CV_64FC1);
    BIntegral = BIntegral(cv::Rect(1, 1, width, height));
    cv::integral(varInv, varInv_integral, CV_64FC1);
    varInv_integral = varInv_integral(cv::Rect(1, 1, width, height));

    for (int i = 0; i < height; ++i) {
        int yStart = std::max(i - radius, 0);
        int yEnd = std::min(i + radius, height - 1);
        auto IDataStart1 = AIntegral.ptr<double>(yStart);
        auto IDataStart2 = BIntegral.ptr<double>(yStart);
        auto IDataEnd1 = AIntegral.ptr<double>(yEnd);
        auto IDataEnd2 = BIntegral.ptr<double>(yEnd);
        auto IDataStartInv = varInv_integral.ptr<double>(yStart);
        auto IDataEndInv = varInv_integral.ptr<double>(yEnd);
        auto imageData = image.ptr<uchar>(i);
        //auto aData = AMat.ptr<float>(i);
        //auto bData = BMat.ptr<float>(i);
        for (int j = 0; j < width; ++j) {
            int xStart = std::max(j - radius, 0);
            int xEnd = std::min(j + radius, width - 1);
            double mean = (IDataEnd1[xEnd] - IDataEnd1[xStart] + IDataStart1[xStart] - IDataStart1[xEnd]);
            double mean2 = (IDataEnd2[xEnd] - IDataEnd2[xStart] + IDataStart2[xStart] - IDataStart2[xEnd]);
            double sumW = (IDataEndInv[xEnd] - IDataEndInv[xStart] + IDataStartInv[xStart] - IDataStartInv[xEnd]);
            double a = mean / sumW;
            double b = mean2 / sumW;
            imageData[j] = cv::saturate_cast<uchar>(a * imageData[j] + b);
            //aData[j] = a;
            //bData[j] = b;
        }
    }
}