// define head function
#ifndef PS_ALGORITHM_H_INCLUDED
#define PS_ALGORITHM_H_INCLUDED#include <iostream>
#include <string>
#include "cv.h"
#include "highgui.h"
#include "cxmat.hpp"
#include "cxcore.hpp"using namespace std;
using namespace cv;void Show_Image(Mat&, const string &);#endif // PS_ALGORITHM_H_INCLUDED#include "PS_Algorithm.h"void Dark_Channel(Mat& src, Mat& dst, int Block_size);
double Atmosperic_Light(Mat& J_dark, Mat& Img);
void Recove_Img(Mat& src, Mat& dst, Mat& T, float Th, float A);int main(void)
{Mat Img;Img=imread("5.jpg");Mat D_Img(Img.size(), CV_32FC3);Img.convertTo(D_Img, CV_32FC3);Mat Dark_Img(D_Img.size(), CV_32FC1);imshow("Img", Img);int Block_size=3;Dark_Channel(D_Img, Dark_Img, Block_size);float A=0;A=Atmosperic_Light(Dark_Img, D_Img);float W=0.9;Mat T(D_Img.size(), CV_32FC1);T=1-W/A*Dark_Img;//imshow("Img", T);float Th=0.35;Mat Img_out(D_Img.size(), CV_32FC3);Recove_Img(D_Img, Img_out, T, Th, A);Img_out/=255;imshow("Out",Img_out);waitKey();cvDestroyAllWindows();cout<<"All is well."<<endl;
}void Dark_Channel(Mat& src, Mat& dst, int Block_size)
{Mat R(src.size(), CV_32FC1);Mat G(src.size(), CV_32FC1);Mat B(src.size(), CV_32FC1);Mat m_array[]={R,G,B};cv::split(src, m_array);int t=0;t=(Block_size-1)/2;Mat a1(Block_size, Block_size, CV_32FC1);Mat a2(Block_size, Block_size, CV_32FC1);Mat a3(Block_size, Block_size, CV_32FC1);double min_a1=0;double min_a2=0;double min_a3=0;double min_value=0;for(int i=t; i<dst.rows-t; i++){for(int j=t; j<dst.cols-t; j++){a1=R(Range(i-t,i+t+1), Range(j-t,j+t+1));a2=G(Range(i-t,i+t+1), Range(j-t,j+t+1));a3=B(Range(i-t,i+t+1), Range(j-t,j+t+1));cv::minMaxLoc(a1, &min_a1,NULL,NULL,NULL);cv::minMaxLoc(a2, &min_a2,NULL,NULL,NULL);cv::minMaxLoc(a3, &min_a3,NULL,NULL,NULL);min_value=min(min_a1, min_a2);min_value=min(min_a3, min_value);dst.at<float>(i,j)=(float)min_value;}}dst(Range(0,t), Range::all())=dst(Range(t,2*t), Range::all());dst(Range(dst.rows-t,dst.rows), Range::all())=dst(Range(dst.rows-(2*t),dst.rows-t), Range::all());dst(Range::all(), Range(0,t))=dst(Range::all(),Range(t,2*t));dst(Range::all(),Range(dst.cols-t,dst.cols))=dst(Range::all(), Range(dst.cols-2*t,dst.cols-t));}double Atmosperic_Light(Mat& J_dark, Mat& Img)
{Mat M1(J_dark.size(), CV_32FC1);M1=J_dark;M1.reshape(0,1);Mat M2(1,J_dark.rows*J_dark.cols, CV_32FC1);cv::sort(M1,M2,CV_SORT_ASCENDING);int Index=J_dark.rows*J_dark.cols*0.9999;float T_value=M2.at<float>(0, Index);float value=0;float Temp_value;float r_temp, g_temp, b_temp;for(int i=0; i<Img.rows; i++){for(int j=0; j<Img.cols; j++){Temp_value=J_dark.at<float>(i,j);if(Temp_value>T_value){r_temp=Img.at<Vec3f>(i,j)[0];g_temp=Img.at<Vec3f>(i,j)[1];b_temp=Img.at<Vec3f>(i,j)[2];Temp_value=(r_temp+g_temp+b_temp)/3.0;value=max(value, Temp_value);}}}return value;
}void Recove_Img(Mat& src, Mat& dst, Mat& T, float Th, float A)
{float value=0;for(int i=0; i<src.rows; i++){for(int j=0; j<src.cols; j++){value=max(Th, T.at<float>(i,j));dst.at<Vec3f>(i,j)[0]=(src.at<Vec3f>(i,j)[0]-A)/value+A;dst.at<Vec3f>(i,j)[1]=(src.at<Vec3f>(i,j)[1]-A)/value+A;dst.at<Vec3f>(i,j)[2]=(src.at<Vec3f>(i,j)[2]-A)/value+A;}}}