#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <cmath>
#include <algorithm>
#include <opencv2\opencv.hpp>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<Mat> 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 < rows; r++) {const uchar *data = rchannel.ptr<uchar>(r);for (unsigned int c = 0; c < cols; 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 < pixelTot; 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<uchar>(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 < pixelTot; 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<uchar>(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 " << duration << "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 < nr; 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 < nl; 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 < nr; 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 < nl; 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 < nr; 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 < nl; 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 < nr; 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 < nl; c++){for (int i=0;i<3;i++){*recoverPtr = (*totalPtr - *aPtr) /(1-(*aPtr / *array));recoverPtr++;totalPtr++;aPtr++;}}}return recoverimage;
}