
#include <iostream>
#include "opencv2/opencv.hpp"
#include <stdio.h>
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <stdlib.h>

using namespace cv;
using namespace std;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
float cvMax(cv::Mat& mat) //get the max value inside a matrix
{
	mat.convertTo(mat,CV_32F);
	float max=mat.at<float>(0,0);
	float* pData=(float*)mat.data;
	for(int i=0;i<mat.rows;i++)
	{
		for(int j=0;j<mat.cols;j++)
		{
			float value = pData[j+i*mat.cols];
			if(value>max) max=value;
		}
	}
	return max;
}

float cvMin(cv::Mat& mat) //get the min value inside a matrix
{
	mat.convertTo(mat,CV_32F);
	float min=mat.at<float>(0,0);
	float* pData=(float*)mat.data;
	for(int i=0;i<mat.rows;i++)
	{
		for(int j=0;j<mat.cols;j++)
		{
			float value = pData[j+i*mat.cols];
			if(value<min) min=value;
		}
	}
	return min;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void scale01(cv::Mat& in) //scale matrix values in [0,1]
{
	in.convertTo(in,CV_32F);
	float min,max; min = cvMin(in); max = cvMax(in);

	for(int i=0; i<in.rows; i++)
	{
		for(int j=0; j<in.cols; j++)
		{
			in.at<float>(i,j)=(in.at<float>(i,j)-min)/(max-min);
		}
	}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void preprocessing(cv::Mat& image, cv::Mat& dest1, cv::Mat& dest2) //preprocessing phase
{
	imwrite("00image.tif",image);
	Mat channel[3]; split(image, channel); image = channel[1];

	//---------------------------------I-gamma: //central vessel light removal
	Mat gamma(image.rows,image.cols,CV_8UC1);
	Mat element = getStructuringElement(1, Size(3,3), Point( -1, -1 ) );
	morphologyEx(image,gamma,2,element );
	//medianBlur(gamma,gamma,3);
	//GaussianBlur(gamma, gamma,Size(9,9),3.24, 0); //implemented in the original paper
	imwrite("01gamma.tif",gamma);

	//---------------------------------I-B: //background image
	Mat b; medianBlur(image,b,69);
	imwrite("02b.tif",b);

	//----------------------------------image_d (sc) //difference image
	gamma.convertTo(gamma,CV_32F); b.convertTo(b,CV_32F);
	Mat d; 	subtract(gamma, b, d);
	scale01(d); d=d*255; d.convertTo(d,CV_8U);
	imwrite("03sc.tif",d);

	//----------------------------------image_h //homogenized image
	double histogram[256];
	for(int i=0;i<d.rows;i++) for(int j=0;j<d.cols;j++) histogram[d.at<uchar>(i,j)]++;

	int maxhist = histogram[d.at<uchar>(0,0)];
	int greyinmax = d.at<uchar>(0,0);
	for(int i=0; i<d.rows; i++)
	{
		for(int j=0; j<d.cols; j++)
		{
			if(histogram[d.at<uchar>(i,j)]>maxhist)
			{
				maxhist = histogram[d.at<uchar>(i,j)];
				greyinmax = d.at<uchar>(i,j);
			}
		}
	}

	int g=0;
	for(int i=0; i<d.cols; i++)
	{
		for(int j=0; j<d.rows; j++)
		{
			g = d.at<unsigned char>(i,j)+128-greyinmax;
			if(g<0) d.at<unsigned char>(i,j)=0;
			else if(g>255) d.at<unsigned char>(i,j)=255;
			else d.at<unsigned char>(i,j)=g;
		}
	}

	//I_VE: //vessel enhancement image
	dest1 = d;

	Mat d_not, d_not_transf;
	bitwise_not(d,d_not);

    Mat element2 = getStructuringElement( 0, Size(17,17), Point( -1, -1 ));
    morphologyEx(d_not,d_not_transf, 2, element2 );

    d.convertTo(d,CV_32F); d_not.convertTo(d_not,CV_32F);  d_not_transf.convertTo(d_not_transf,CV_32F);

    subtract(d_not,d_not_transf,d);
    scale01(d); d=d*255;
    d.convertTo(d,CV_8U);
    imwrite("05veve.tif",d);

    dest2 = d;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

int main( int argc, char** argv )
{
	//list of strings for recursive preprocessing over DRIVE database

	string listin[] = {"01_tes.tif","02_tes.tif","03_tes.tif","04_tes.tif","05_tes.tif","06_tes.tif","07_tes.tif","08_tes.tif","09_tes.tif","10_tes.tif",
						"11_tes.tif","12_tes.tif","13_tes.tif","14_tes.tif","15_tes.tif","16_tes.tif","17_tes.tif","18_tes.tif","19_tes.tif","20_tes.tif",
						"21_tes.tif","22_tes.tif","23_tes.tif","24_tes.tif","25_tes.tif","26_tes.tif","27_tes.tif","28_tes.tif","29_tes.tif","30_tes.tif",
						"31_tes.tif","32_tes.tif","33_tes.tif","34_tes.tif","35_tes.tif","36_tes.tif","37_tes.tif","38_tes.tif","39_tes.tif","40_tes.tif",};

	string listtifh[] = {"01h.tif","02h.tif","03h.tif","04h.tif","05h.tif","06h.tif","07h.tif","08h.tif","09h.tif","10h.tif",
						 "11h.tif","12h.tif","13h.tif","14h.tif","15h.tif","16h.tif","17h.tif","18h.tif","19h.tif","20h.tif",
						 "21h.tif","22h.tif","23h.tif","24h.tif","25h.tif","26h.tif","27h.tif","28h.tif","29h.tif","30h.tif",
						 "31h.tif","32h.tif","33h.tif","34h.tif","35h.tif","36h.tif","37h.tif","38h.tif","39h.tif","40h.tif",};

	string listtifve[] = {"01ve.tif","02ve.tif","03ve.tif","04ve.tif","05ve.tif","06ve.tif","07ve.tif","08ve.tif","09ve.tif","10ve.tif",
						 "11ve.tif","12ve.tif","13ve.tif","14ve.tif","15ve.tif","16ve.tif","17ve.tif","18ve.tif","19ve.tif","20ve.tif",
						 "21ve.tif","22ve.tif","23ve.tif","24ve.tif","25ve.tif","26ve.tif","27ve.tif","28ve.tif","29ve.tif","30ve.tif",
						 "31ve.tif","32ve.tif","33ve.tif","34ve.tif","35ve.tif","36ve.tif","37ve.tif","38ve.tif","39ve.tif","40ve.tif",};

	//recursive image H and VE generation
	for(int i=0; i<sizeof(listtifh)/4;i++)
	{
		string name = listin[i];
		string saveh = listtifh[i];
		string saveve = listtifve[i];

		Mat image = cv::imread(name),h,ve;
		preprocessing(image,h,ve);

		imwrite(saveh,h);
		imwrite(saveve,ve);
	}

    waitKey(0);
    return 0;
}
