//piece of code for features extraction
#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) //return max value of 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) //return min value of 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;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Mat computeFeatures(cv::Mat& h, cv::Mat& ve, int& w, int& w2, int& n_features) //extract the 7 features
{
	//w2 = 17, w = 9;
	h.convertTo(ve, CV_32F);
	ve.convertTo(h, CV_32F);

	Mat temp(w,w,CV_32FC1);
	Mat features(n_features,(h.rows-w2+1)*(h.cols-w2+1),CV_32FC1);

	int counter = 0;
	for(int i=0; i<h.cols-w2+1;i++) //extraction out of I_H, f1 f2 f3 f4 f5
	{
		for(int j=0; j<h.rows-w2+1;j++)
		{
			temp = h(cv::Rect(i+floor(w2/2)-floor(w/2),j+floor(w2/2)-floor(w/2),w,w));
			float max2=cvMax(temp); float min2=cvMin(temp);

			features.at<float>(0,counter)=h.at<float>(i+floor(w2/2),j+floor(w2/2))-min2;
			features.at<float>(1,counter)=max2-h.at<float>(i+floor(w2/2),j+floor(w2/2));

			Scalar mu, var; meanStdDev(temp, mu, var);
			features.at<float>(2,counter)=h.at<float>(i+floor(w2/2),j+floor(w2/2))-mu[0];
			features.at<float>(3,counter)=var[0];
			features.at<float>(4,counter)=h.at<float>(i+floor(w2/2),j+floor(w2/2));

			counter++;
		}
	}
	//extraction out of I_VE, f6 f7
	//gaussian valued 17x17 image

	Mat gauss(w2,w2,CV_32FC1);
	float var2 = (1.7)*(1.7);
	for(int i=0; i<gauss.cols; i++)
	{
		for(int j=0; j<gauss.rows; j++)
		{
			gauss.at<float>(j,i)=exp(-1*((i-floor(w2/2))*(i-floor(w2/2))+(j-floor(w2/2))*(j-floor(w2/2)))/(2*var2));
		}
	}

    //I_VE
	Mat temp2(w2,w2,CV_32FC1),usata(w2,w2,CV_32FC1);
	int counter2 = 0;

	for(int i=0; i<ve.cols-w2+1;i++)
	{
		for(int j=0; j<ve.rows-w2+1;j++)
		{
			temp2 = ve(cv::Rect(i,j,w2,w2));
			usata = temp2.mul(gauss);

			//-----------m0,0
			Scalar somma = sum(usata);
			float m00 = somma[0];
			if (m00 == 0) m00 = somma[0]+1;

			//-----------m1,0
			float m10 = 0;
			for(int a=0; a<w; a++) // a --> rows
			{
				for(int b=0; b<w; b++)
				{
					m10 = m10 + a*usata.at<float>(a,b);
				}
			}
			//-----------m0,1
			float m01 = 0;
			for(int c=0; c<w; c++) //c --> rows
			{
				for(int d=0; d<w; d++)
				{
					m01 = m01 + d*usata.at<float>(c,d);
				}
			}
			//-----------ni20
			float ni20 = 0;
			for(int e=0; e<w; e++) // e --> rows
			{
				for(int f=0; f<w; f++)
				{
					ni20 = ni20 + ((1/m00)*(1/m00))*((e-m10/m00)*(e-m10/m00)*usata.at<float>(e,f));
				}
			}
			//-----------ni02
			float ni02 = 0;
			for(int g=0; g<w; g++) // g --> rows
			{
				for(int h=0; h<w; h++)
				{
					ni02 = ni02 + ((1/m00)*(1/m00))*((h-m01/m00)*(h-m01/m00)*usata.at<float>(g,h));
				}
			}
			//-----------ni11
			float ni11 = 0;
			for(int m=0; m<w; m++) // m --> rows
			{
				for(int n=0; n<w; n++)
				{
					ni11 = ni11 + ((1/m00)*(1/m00))*((m-m10/m00)*(n-m01/m00)*usata.at<float>(m,n));
				}
			}
			//-----------fi1 e fi2
			float fi1, fi2;
			fi1 = ni20 + ni02;
			fi2 = fi1*fi1 + 4*((1/m00)*(1/m00)*(1/m00)*(1/m00))*ni11*ni11;

			features.at<float>(5,counter2)=abs(log(fi1));
			features.at<float>(6,counter2)=abs(log(fi2));
			counter2++;
		}
	}

	//features normalization
	Mat featrigatemp(1,features.cols,CV_32FC1);
	Scalar media, varianza;
	for(int i=0; i<features.rows;i++)
	{
		for(int j=0; j<features.cols;j++) featrigatemp.at<float>(0,j) = features.at<float>(i,j);
		meanStdDev(featrigatemp, media, varianza);
		for(int h=0; h<features.cols; h++) features.at<float>(i,h)=(features.at<float>(i,h)-media[0])/varianza[0];
	}
	transpose(features,features);
	return features;
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//----------------------------------------------------------------------------------main

int main( int argc, char** argv )
{
	//parameters
	int w=9, w2=17, n_features=7;

	//list of strings for recursive features generation over DRIVE database

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

	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 features generation

	for(int i=0; i<sizeof(listyml)/4;i++)
	{
		string nameh = listtifh[i];
		string nameve = listtifve[i];
		string save = listyml[i];

		Mat h = imread(nameh,0);
		Mat ve = imread(nameve,0);
		Mat features = computeFeatures(h,ve,w,w2,n_features);

		cv::FileStorage storage(save, cv::FileStorage::WRITE);
		storage << "features" << features;
		storage.release();
	}
}
