基于偏振差分成像的图像去雾算法

#include 
#include 
#include 
#include 
#include 
#include 

using namespace cv;
using namespace std;

typedef struct Pixel 
{
    int x, y;
    float data;
}Pixel;

bool structCmp(const Pixel &a, const Pixel &b) 
{
    return a.data > b.data;//descending降序
}

Mat getAdd(Mat& parimg, Mat& perimg);
Mat getSubstract(Mat& parimg, Mat& perimg);
Mat getA(Mat& parimg, Mat& perimg, float *array);
Mat recover(Mat& total, Mat& a, float *array);

int main() 
{
    string loc1 = "E:\\polarization\\firstopencv\\0degree.tiff";
    string loc2 = "E:\\polarization\\firstopencv\\90degree.tiff";
    string name = "forest";
    clock_t start, finish;
    double duration;

    cout <<"A defog program" << endl
        <<"----------------" << endl;

    Mat src1 = imread(loc1);
    Mat src2 = imread(loc2);
    Mat parimage = Mat(src1.size(), CV_16UC3);
    Mat perimage = Mat(src2.size(), CV_16UC3);
    src1.convertTo(parimage, CV_16UC3, 255.0 );
    src2.convertTo(perimage, CV_16UC3, 255.0 );
    imshow("0degree", parimage);
    imshow("90degree", perimage);
    cout <<"input hazy image" << endl;

    int rows = parimage.rows;
    int cols = parimage.cols;
    Mat newparimage; Mat newperimage;
    parimage.convertTo(newparimage, CV_32FC3, 1.0/65535.0, 0);
    perimage.convertTo(newperimage, CV_32FC3, 1.0/65535.0, 0);


    Mat totalimage(rows, cols, CV_32FC3);
    Mat diffimage(rows, cols, CV_32FC3);
    Mat A(rows, cols, CV_32FC3);
    Mat recoverimg(rows, cols, CV_32FC3);
    totalimage = getAdd(newparimage,newperimage);
    diffimage = getSubstract(newparimage, newperimage);
    int nr = rows;int nl = cols;
    start = clock();

    //estimate Airlight
    //开一个结构体数组存暗通道，再sort，取最大0.1%，利用结构体内存储的原始坐标在原图中取点
    cout <<"estimating airlight." << endl;
    vector planes;
    split(diffimage, planes);
    Mat rchannel = planes[2];
    int pixelTot = rows * cols * 0.001;
    float *Apar = new float[3]; float *Aper = new float[3];
    Pixel *toppixels, *allpixels;
    toppixels = new Pixel[pixelTot];
    allpixels = new Pixel[rows * cols];


    for (unsigned int r = 0; r ) 
    {
        const uchar *data = rchannel.ptr(r);
        for (unsigned int c = 0; c ) 
        {
            allpixels[r*cols + c].data = *data;
            allpixels[r*cols + c].x = r;
            allpixels[r*cols + c].y = c;
        }
    }
    
    std::sort(allpixels, allpixels + rows * cols, structCmp);
    memcpy(toppixels, allpixels, pixelTot * sizeof(Pixel));
    
    float tmp_par_A[3]; float tmp_per_A[3];
    float A_par_r, A_par_g, A_par_b, avg1, maximum1 = 0;
    int idx1, idy1, max_x1, max_y1;
    for (int i = 0; i ) 
    {
        float A_par_sum_r = 0, A_par_sum_g = 0, A_par_sum_b = 0;
        idx1 = allpixels[i].x; idy1 = allpixels[i].y;
        const uchar *pardata = newparimage.ptr(idx1);
        pardata += 3 * idy1;
        A_par_b = *pardata++;
        A_par_g = *pardata++;
        A_par_r = *pardata++;

        A_par_sum_b += A_par_b;
        A_par_sum_g += A_par_g;
        A_par_sum_r += A_par_r;
        
        tmp_par_A[0] = A_par_sum_b / pixelTot;
        tmp_par_A[1] = A_par_sum_g / pixelTot;
        tmp_par_A[2] = A_par_sum_r / pixelTot;
    }

    float A_per_r, A_per_g, A_per_b, avg2, maximum2 = 0;
    int idx2, idy2, max_x2, max_y2;
    for (int i = 0; i )
    {
        float A_per_sum_r = 0, A_per_sum_g = 0, A_per_sum_b = 0;
        idx2 = allpixels[i].x; idy2 = allpixels[i].y;
        const uchar *perdata = newperimage.ptr(idx2);
        perdata += 3 * idy2;
        A_per_b = *perdata++;
        A_per_g = *perdata++;
        A_per_r = *perdata++;
        
        A_per_sum_b += A_per_b;
        A_per_sum_g += A_per_g;
        A_per_sum_r += A_per_r;

        
        tmp_per_A[0] = A_per_sum_b / pixelTot;
        tmp_per_A[1] = A_per_sum_g / pixelTot;
        tmp_per_A[2] = A_per_sum_r / pixelTot;
    }

    delete[] toppixels;
    delete[] allpixels;

    float P[3]; float Ainf[3];
    for (int i = 0; i <3; i++)
    {
        P[i] = (tmp_per_A[i] - tmp_par_A[i]) / (tmp_per_A[i] + tmp_par_A[i]);
        Ainf[i] = tmp_per_A[i] + tmp_par_A[i];
    }

    A = getA(newparimage, newperimage, P);

    cout <<"start recovering." << endl;
    recoverimg = recover(totalimage, A, Ainf);
    cout <<"recovering finished." << endl;

    finish = clock();
    duration = (double)(finish - start);
    cout <<"defog used " <"ms time;" << endl;
    imshow("result", recoverimg);
    waitKey(0);

    
    imwrite(name + "_refined.png", recoverimg);
    destroyAllWindows();
    parimage.release();
    perimage.release();
    recoverimg.release();
    return 0;
}

Mat getAdd(Mat& parimg, Mat& perimg)
{
    int nr = parimg.rows; int nl = perimg.cols;
    Mat addimg = Mat::zeros(nr, nl, CV_32FC3);
    for (unsigned int r = 0; r )
    {
        const float* parPtr = parimg.ptr<float>(r);
        const float* perPtr = perimg.ptr<float>(r);
        float* addPtr = addimg.ptr<float>(r);
        for (unsigned int c = 0; c )
        {
            for (int i = 0; i <3; i++)
            {
                *addPtr++ = *parPtr++ + *perPtr++;
            }
        }
    }
    return addimg;
}

Mat getSubstract(Mat& parimg, Mat& perimg)
{
    int nr = parimg.rows; int nl = perimg.cols;
    Mat subimg = Mat::zeros(nr, nl, CV_32FC3);
    for (unsigned int r = 0; r )
    {
        const float* parPtr = parimg.ptr<float>(r);
        const float* perPtr = perimg.ptr<float>(r);
        float* subPtr = subimg.ptr<float>(r);
        for (unsigned int c = 0; c )
        {
            for (int i = 0; i <3; i++)
            {
                *subPtr++ = *perPtr++ - *parPtr++;
            }
        }
    }
    return subimg;
}

Mat getA(Mat& parimg, Mat& perimg, float *array)
{
    int nr = parimg.rows; int nl = perimg.cols;
    Mat A= Mat::zeros(nr, nl, CV_32FC3);
    for (unsigned int r = 0; r )
    {
        const float* parPtr = parimg.ptr<float>(r);
        const float* perPtr = perimg.ptr<float>(r);
        float* APtr = A.ptr<float>(r);
        for (unsigned int c = 0; c )
        {
            for (int i = 0; i <3; i++)
            {
                *APtr++ = (*perPtr++ - *parPtr++) / *array;
            }
        }
    }
    return A;
}

Mat recover(Mat& total,Mat& a,float *array)
{
    int nr = total.rows;int nl=total.cols;
    Mat recoverimage = Mat::zeros(nr, nl, CV_32FC3);
    for (unsigned int r = 0; r )
    {
        const float* totalPtr = total.ptr<float>(r);
        const float* aPtr = a.ptr<float>(r);
        float* recoverPtr = recoverimage.ptr<float>(r);
        for (unsigned int c = 0; c )
        {
            for (int i=0;i<3;i++)
            {
                *recoverPtr = (*totalPtr - *aPtr) /(1-(*aPtr / *array));
                recoverPtr++;
                totalPtr++;
                aPtr++;
            }
        }
    }
    return recoverimage;
}