// 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;

        }

    }

}