跳转至

前言

之前写过这篇论文的复现,https://blog.csdn.net/just_sort/article/details/84110518 ,但是当时估算透射率的时候没有用softmating或者导向滤波,是用何博士的公式粗略的计算的,所幸,在了解了导向滤波后,我完成了使用导向滤波估算 更加精细的透射率的任务,也确实取得了更好的去雾效果,所以打算分享一下这个源码。我所有开源的论文实现均在github地址里:https://github.com/BBuf/-Image-processing-algorithm

算法原理

这些就不说了,可以查看之前的博客,直接公布代码吧。

代码

#include <opencv2/opencv.hpp>
#include <iostream>
#include <algorithm>
#include <vector>
using namespace cv;
using namespace std;

cv::Mat guidedFilter(cv::Mat I, cv::Mat p, int r, double eps)
{
    /*
    % GUIDEDFILTER   O(1) time implementation of guided filter.
    %
    %   - guidance image: I (should be a gray-scale/single channel image)
    %   - filtering input image: p (should be a gray-scale/single channel image)
    %   - local window radius: r
    %   - regularization parameter: eps
    */

    cv::Mat _I;
    I.convertTo(_I, CV_64FC1);
    I = _I;

    cv::Mat _p;
    p.convertTo(_p, CV_64FC1);
    p = _p;

    //[hei, wid] = size(I);
    int hei = I.rows;
    int wid = I.cols;

    //N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.
    cv::Mat N;
    cv::boxFilter(cv::Mat::ones(hei, wid, I.type()), N, CV_64FC1, cv::Size(r, r));

    //mean_I = boxfilter(I, r) ./ N;
    cv::Mat mean_I;
    cv::boxFilter(I, mean_I, CV_64FC1, cv::Size(r, r));

    //mean_p = boxfilter(p, r) ./ N;
    cv::Mat mean_p;
    cv::boxFilter(p, mean_p, CV_64FC1, cv::Size(r, r));

    //mean_Ip = boxfilter(I.*p, r) ./ N;
    cv::Mat mean_Ip;
    cv::boxFilter(I.mul(p), mean_Ip, CV_64FC1, cv::Size(r, r));

    //cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.
    cv::Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);

    //mean_II = boxfilter(I.*I, r) ./ N;
    cv::Mat mean_II;
    cv::boxFilter(I.mul(I), mean_II, CV_64FC1, cv::Size(r, r));

    //var_I = mean_II - mean_I .* mean_I;
    cv::Mat var_I = mean_II - mean_I.mul(mean_I);

    //a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper; 
    cv::Mat a = cov_Ip / (var_I + eps);

    //b = mean_p - a .* mean_I; % Eqn. (6) in the paper;
    cv::Mat b = mean_p - a.mul(mean_I);

    //mean_a = boxfilter(a, r) ./ N;
    cv::Mat mean_a;
    cv::boxFilter(a, mean_a, CV_64FC1, cv::Size(r, r));
    mean_a = mean_a / N;

    //mean_b = boxfilter(b, r) ./ N;
    cv::Mat mean_b;
    cv::boxFilter(b, mean_b, CV_64FC1, cv::Size(r, r));
    mean_b = mean_b / N;

    //q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
    cv::Mat q = mean_a.mul(I) + mean_b;

    return q;
}


int rows, cols;
//获取最小值矩阵
int **getMinChannel(cv::Mat img) {
    rows = img.rows;
    cols = img.cols;
    if (img.channels() != 3) {
        fprintf(stderr, "Input Error!");
        exit(-1);
    }
    int **imgGray;
    imgGray = new int *[rows];
    for (int i = 0; i < rows; i++) {
        imgGray[i] = new int[cols];
    }
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            int loacalMin = 255;
            for (int k = 0; k < 3; k++) {
                if (img.at<Vec3b>(i, j)[k] < loacalMin) {
                    loacalMin = img.at<Vec3b>(i, j)[k];
                }
            }
            imgGray[i][j] = loacalMin;
        }
    }
    return imgGray;
}

//求暗通道
int **getDarkChannel(int **img, int blockSize = 3) {
    if (blockSize % 2 == 0 || blockSize < 3) {
        fprintf(stderr, "blockSize is not odd or too small!");
        exit(-1);
    }
    //计算pool Size
    int poolSize = (blockSize - 1) / 2;
    int newHeight = rows + poolSize - 1;
    int newWidth = cols + poolSize - 1;
    int **imgMiddle;
    imgMiddle = new int *[newHeight];
    for (int i = 0; i < newHeight; i++) {
        imgMiddle[i] = new int[newWidth];
    }
    for (int i = 0; i < newHeight; i++) {
        for (int j = 0; j < newWidth; j++) {
            if (i < rows && j < cols) {
                imgMiddle[i][j] = img[i][j];
            }
            else {
                imgMiddle[i][j] = 255;
            }
        }
    }
    int **imgDark;
    imgDark = new int *[rows];
    for (int i = 0; i < rows; i++) {
        imgDark[i] = new int[cols];
    }
    int localMin = 255;
    for (int i = poolSize; i < newHeight - poolSize; i++) {
        for (int j = poolSize; j < newWidth - poolSize; j++) {
            for (int k = i - poolSize; k < i + poolSize + 1; k++) {
                for (int l = j - poolSize; l < j + poolSize + 1; l++) {
                    if (imgMiddle[k][l] < localMin) {
                        localMin = imgMiddle[k][l];
                    }
                }
            }
            imgDark[i - poolSize][j - poolSize] = localMin;
        }
    }
    return imgDark;
}

struct node {
    int x, y, val;
    node() {}
    node(int _x, int _y, int _val) :x(_x), y(_y), val(_val) {}
    bool operator<(const node &rhs) {
        return val > rhs.val;
    }
};

//估算全局大气光值
int getGlobalAtmosphericLightValue(int **darkChannel, cv::Mat img, bool meanMode = false, float percent = 0.001) {
    int size = rows * cols;
    std::vector <node> nodes;
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            node tmp;
            tmp.x = i, tmp.y = j, tmp.val = darkChannel[i][j];
            nodes.push_back(tmp);
        }
    }
    sort(nodes.begin(), nodes.end());
    int atmosphericLight = 0;
    if (int(percent*size) == 0) {
        for (int i = 0; i < 3; i++) {
            if (img.at<Vec3b>(nodes[0].x, nodes[0].y)[i] > atmosphericLight) {
                atmosphericLight = img.at<Vec3b>(nodes[0].x, nodes[0].y)[i];
            }
        }
    }
    //开启均值模式
    if (meanMode == true) {
        int sum = 0;
        for (int i = 0; i < int(percent*size); i++) {
            for (int j = 0; j < 3; j++) {
                sum = sum + img.at<Vec3b>(nodes[i].x, nodes[i].y)[j];
            }
        }
    }
    //获取暗通道在前0.1%的位置的像素点在原图像中的最高亮度值
    for (int i = 0; i < int(percent*size); i++) {
        for (int j = 0; j < 3; j++) {
            if (img.at<Vec3b>(nodes[i].x, nodes[i].y)[j] > atmosphericLight) {
                atmosphericLight = img.at<Vec3b>(nodes[i].x, nodes[i].y)[j];
            }
        }
    }
    return atmosphericLight;
}

//恢复原图像
// Omega 去雾比例 参数
//t0 最小透射率值
cv::Mat getRecoverScene(cv::Mat img, float omega = 0.95, float t0 = 0.1, int blockSize = 15, bool meanModel = false, float percent = 0.001) {
    int** imgGray = getMinChannel(img);
    int **imgDark = getDarkChannel(imgGray, blockSize = blockSize);
    int atmosphericLight = getGlobalAtmosphericLightValue(imgDark, img, meanModel = meanModel, percent = percent);
    float **imgDark2, **transmission;
    imgDark2 = new float *[rows];
    for (int i = 0; i < rows; i++) {
        imgDark2[i] = new float[cols];
    }
    Mat B(rows, cols, CV_8UC1);
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            B.at<uchar>(i, j) = img.at<Vec3b>(i, j)[0];
        }
    }
    Mat p(rows, cols, CV_64FC1);
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            p.at<double>(i, j) = 1 - omega * imgDark[i][j] / atmosphericLight;
        }
    }
    Mat transmission_filter = guidedFilter(B, p, 80, 1e-3);
    imshow("transmission", transmission_filter);
    imwrite("F:\\res3.jpg", transmission_filter);
    waitKey(0);
    transmission = new float *[rows];
    for (int i = 0; i < rows; i++) {
        transmission[i] = new float[cols];
    }
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            imgDark2[i][j] = float(imgDark[i][j]);
            //transmission[i][j] = 1 - omega * imgDark[i][j] / atmosphericLight;
            transmission[i][j] = transmission_filter.at<double>(i, j);
            if (transmission[i][j] < 0.1) {
                transmission[i][j] = 0.1;
            }
        }
    }
    cv::Mat dst(img.rows, img.cols, CV_8UC3);
    for (int channel = 0; channel < 3; channel++) {
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                int temp = (img.at<Vec3b>(i, j)[channel] - atmosphericLight) / transmission[i][j] + atmosphericLight;
                if (temp > 255) {
                    temp = 255;
                }
                if (temp < 0) {
                    temp = 0;
                }
                dst.at<Vec3b>(i, j)[channel] = temp;
            }
        }
    }
    return dst;
}

int main() {
    Mat src = imread("F:\\fog\\3.jpg");
    //cv::Mat src = cv::imread("/home/zxy/CLionProjects/Acmtest/3.jpg");
    rows = src.rows;
    cols = src.cols;
    cv::Mat dst = getRecoverScene(src);
    cv::imshow("origin", src);
    cv::imshow("result", dst);
    cv::imwrite("F:\\res2.jpg", dst);
    waitKey(0);
}

效果

原图

不使用导向滤波的结果

使用导向滤波的透射率图

使用导向滤波的结果图

结论

可以看到效果确实好一些,关于导向滤波可以看我的这篇博客:https://blog.csdn.net/just_sort/article/details/84324239