灰度归一化算法比较(SQI,TT等)

gross.h

#ifndef gross_h 
#define gross_h 

void GrossBrajovic(IplImage *src, std::string dstfilename);
void Gross(IplImage *src, std::string filename);
void fill0(double **u, int n);
void interp(double **uf, double **uc, int nf, int nc);
void rstrct(double **uc, double **uf, int nc, int nf);
void sor(double **a, double **b, double **c, double **d, double **e,
	double **f, double **u, int jmax, double rjac);



#endif#pragma once

gross.cpp

#include "stdafx.h"
#include "gross.h"
#define MAXITS 1000
#define EPS 1.0e-5
#define pi 3.1415
using namespace std;

void GrossBrajovic(IplImage *src, std::string dstfilename)
{
	int i;
	int j, l;
	int jmax = src->width; //The image size is width*width.
	double lamta = 2.0;
	double rjac = cos(pi / jmax);
	double **a = (double**)malloc(jmax * sizeof(double *));
	for (i = 0; i < jmax; i++) a[i] = (double *)malloc(jmax * sizeof(double));
	double **b = (double**)malloc(jmax * sizeof(double *));
	for (i = 0; i < jmax; i++) b[i] = (double *)malloc(jmax * sizeof(double));
	double **c = (double**)malloc(jmax * sizeof(double *));
	for (i = 0; i < jmax; i++) c[i] = (double *)malloc(jmax * sizeof(double));
	double **d = (double**)malloc(jmax * sizeof(double *));
	for (i = 0; i < jmax; i++) d[i] = (double *)malloc(jmax * sizeof(double));
	double **e = (double**)malloc(jmax * sizeof(double *));
	for (i = 0; i < jmax; i++) e[i] = (double *)malloc(jmax * sizeof(double));
	double **f = (double**)malloc(jmax * sizeof(double *));
	for (i = 0; i < jmax; i++) f[i] = (double *)malloc(jmax * sizeof(double));
	double **L = (double**)malloc(jmax * sizeof(double *));
	for (i = 0; i < jmax; i++) L[i] = (double *)malloc(jmax * sizeof(double));

	IplImage *src_GB = 0;
	IplImage *src_L = 0;
	IplImage *src_GB_out = 0;
	src_GB = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);
	src_L = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	src_GB_out = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);

	for (j = 0; j<jmax; j++)
	{
		for (l = 0; l<jmax; l++)
		{
			f[j][l] = cvGet2D(src, j, l).val[0];
			L[j][l] = a[j][l] = b[j][l] = c[j][l] = d[j][l] = e[j][l] = 0;
		}
	}

	for (j = 1; j<jmax - 1; j++)
	{
		for (l = 1; l<jmax - 1; l++)
		{
			b[j][l] = -lamta * min(f[j][l], f[j - 1][l]) / (fabs(f[j][l] - f[j - 1][l]) + 1);
			d[j][l] = -lamta * min(f[j][l], f[j][l - 1]) / (fabs(f[j][l] - f[j][l - 1]) + 1);
			a[j][l] = -lamta * min(f[j][l], f[j + 1][l]) / (fabs(f[j][l] - f[j + 1][l]) + 1);
			c[j][l] = -lamta * min(f[j][l], f[j][l + 1]) / (fabs(f[j][l] - f[j][l + 1]) + 1);
			e[j][l] = 1 - a[j][l] - b[j][l] - c[j][l] - d[j][l];
		}
	}
	sor(a, b, c, d, e, f, L, jmax, rjac);

	for (j = 1; j<jmax - 1; j++)
	{
		for (l = 1; l<jmax - 1; l++)
		{
			double q = cvGet2D(src, j, l).val[0] / (L[j][l] + 1);
			cvSet2D(src_GB, j, l, cvScalar(q));
			cvSet2D(src_L, j, l, cvScalar(L[j][l]));
		}
	}
	cvConvertScaleAbs(src_GB, src_GB_out, 190, 0);

	//cvNamedWindow("src_GB", CV_WINDOW_AUTOSIZE);
	//cvShowImage("src_GB", src_GB_out);
	//cvNamedWindow( "src_L", CV_WINDOW_AUTOSIZE );
	//cvShowImage( "src_L", src_L );
	cvSaveImage(dstfilename.c_str(), src_GB_out);

	cvReleaseImage(&src_GB);
	cvReleaseImage(&src_GB_out);
	cvReleaseImage(&src_L);
}

void Gross(IplImage *src, std::string filename)
{
	int i;
	int j, l;
	int n = (src->width); //The image size is width*width. 
	double lamta = 2.0;
	double **If = (double**)malloc(n * sizeof(double *));
	for (i = 0; i < n; i++) If[i] = (double *)malloc(n * sizeof(double));
	double **Lf = (double**)malloc(n * sizeof(double *));
	for (i = 0; i < n; i++) Lf[i] = (double *)malloc(n * sizeof(double));
	int nn = (n + 1) / 2;
	double **Ic = (double**)malloc(nn * sizeof(double *));
	for (i = 0; i < nn; i++) Ic[i] = (double *)malloc(nn * sizeof(double));
	double **Lc = (double**)malloc(nn * sizeof(double *));
	for (i = 0; i < nn; i++) Lc[i] = (double *)malloc(nn * sizeof(double));
	double **a = (double**)malloc(nn * sizeof(double *));
	for (i = 0; i < nn; i++) a[i] = (double *)malloc(nn * sizeof(double));
	double **b = (double**)malloc(nn * sizeof(double *));
	for (i = 0; i < nn; i++) b[i] = (double *)malloc(nn * sizeof(double));
	double **c = (double**)malloc(nn * sizeof(double *));
	for (i = 0; i < nn; i++) c[i] = (double *)malloc(nn * sizeof(double));
	double **d = (double**)malloc(nn * sizeof(double *));
	for (i = 0; i < nn; i++) d[i] = (double *)malloc(nn * sizeof(double));
	double **e = (double**)malloc(nn * sizeof(double *));
	for (i = 0; i < nn; i++) e[i] = (double *)malloc(nn * sizeof(double));
	fill0(Lf, n);  //L fine
	fill0(Lc, nn); //L coarse
	fill0(If, n);  //I fine
	fill0(Ic, nn); //I coarse
	fill0(a, nn);
	fill0(b, nn);
	fill0(c, nn);
	fill0(d, nn);
	fill0(e, nn);

	for (j = 0; j<n; j++)
	{
		for (l = 0; l<n; l++)
		{
			If[j][l] = cvGet2D(src, j, l).val[0];
		}
	}

	rstrct(Ic, If, nn, n);

	//½«ÊäÈë½µµ½´ÖÍø¸ñÉÏ£¬½âͬÑùµÄƫ΢·Ö·½³Ì
	for (j = 1; j<nn - 1; j++)
	{
		for (l = 1; l<nn - 1; l++)
		{
			b[j][l] = -lamta * min(Ic[j][l], Ic[j - 1][l]) / (fabs(Ic[j][l] - Ic[j - 1][l]) + 1);
			d[j][l] = -lamta * min(Ic[j][l], Ic[j][l - 1]) / (fabs(Ic[j][l] - Ic[j][l - 1]) + 1);
			a[j][l] = -lamta * min(Ic[j][l], Ic[j + 1][l]) / (fabs(Ic[j][l] - Ic[j + 1][l]) + 1);
			c[j][l] = -lamta * min(Ic[j][l], Ic[j][l + 1]) / (fabs(Ic[j][l] - Ic[j][l + 1]) + 1);
			e[j][l] = 1 - a[j][l] - b[j][l] - c[j][l] - d[j][l];
		}
	}
	double rjac = cos(pi / nn);
	sor(a, b, c, d, e, Ic, Lc, nn, rjac);

	//½«½âÍêµÄ½á¹ûË«ÏßÐÔ²åÖµµ½Ï¸Íø¸ñÉÏ£¬LfµÄֵΪËùÇóµÄL
	interp(Lf, Lc, n, nn);

	IplImage *src_GB = 0;
	IplImage *src_L = 0;
	IplImage *src_GB_out = 0;
	src_GB = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);
	src_L = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	src_GB_out = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);

	for (j = 2; j<n - 2; j++)
	{
		for (l = 2; l<n - 2; l++)
		{
			double q = cvGet2D(src, j, l).val[0] / (Lf[j][l] + 1);
			cvSet2D(src_GB, j, l, cvScalar(q));
			cvSet2D(src_L, j, l, cvScalar(Lf[j][l]));
		}
	}

	cvConvertScaleAbs(src_GB, src_GB_out, 170, 0);
	//cvNamedWindow( "src_L", CV_WINDOW_AUTOSIZE );
	//cvShowImage( "src_L", src_L );
	cvNamedWindow("src_GB_MG", CV_WINDOW_AUTOSIZE);
	cvShowImage("src_GB_MG", src_GB_out);

	//std::string name="G&B_"+filename;
	//cvSaveImage(name.c_str(), src_GB_out);

}

//Successive overrelaxation method for PDE problem
void sor(double **a, double **b, double **c, double **d, double **e,
	double **f, double **u, int jmax, double rjac)
{
	int ipass, j, jsw, l, lsw, n;
	double anorm, anormf = 0.0, omega = 1.0, resid;

	for (j = 1; j<jmax - 1; j++)
		for (l = 1; l<jmax - 1; l++)
			anormf += fabs(f[j][l]);
	for (n = 1; n <= MAXITS; n++)
	{
		anorm = 0.0;
		jsw = 1;
		for (ipass = 1; ipass <= 2; ipass++)
		{
			lsw = jsw;
			for (j = 1; j<jmax - 1; j++)
			{
				for (l = lsw + 1; l<jmax - 1; l += 2)
				{
					resid = a[j][l] * u[j + 1][l]
						+ b[j][l] * u[j - 1][l]
						+ c[j][l] * u[j][l + 1]
						+ d[j][l] * u[j][l - 1]
						+ e[j][l] * u[j][l]
						- f[j][l];
					anorm += fabs(resid);
					u[j][l] -= omega*resid / e[j][l];
				}
				lsw = 3 - lsw;
			}
			jsw = 3 - jsw;
			omega = (n == 1 && ipass == 1 ? 1.0 / (1.0 - 0.5*rjac*rjac) :
				1.0 / (1.0 - 0.25*rjac*rjac*omega));
		}
		if (anorm < EPS*anormf) return;
	}
	cout << "MAXITS exceeded" << endl;
}

//  ²åÖµ
//      | 0.25  0.5  0.25 |
//	P = | 0.5   1    0.5  |
//      | 0.25  0.5  0.25 |
void interp(double **uf, double **uc, int nf, int nc)
{
	int ic, iif, jc, jf;
	for (jc = 1, jf = 2; jc<nc - 1; jc++, jf += 2)
		for (ic = 1; ic<nc - 1; ic++) uf[2 * ic][jf] = uc[ic][jc];
	for (jf = 2; jf<nf - 2; jf += 2)
		for (iif = 3; iif<nf - 2; iif += 2)
			uf[iif][jf] = 0.5*(uf[iif + 1][jf] + uf[iif - 1][jf]);
	for (jf = 3; jf<nf - 2; jf += 2)
		//for (iif = 2; iif<nf - 2; iif += 2)
			uf[iif][jf] = 0.5*(uf[iif][jf + 1] + uf[iif][jf - 1]);
	for (jf = 3; jf<nf - 3; jf += 2)
		for (iif = 3; iif<nf - 3; iif += 2)
			uf[iif][jf] = 0.25*(uf[iif + 1][jf + 1] + uf[iif + 1][jf - 1] + uf[iif - 1][jf - 1] + uf[iif - 1][jf + 1]);
}

void fill0(double **u, int n)
{
	int i, j;
	for (j = 0; j<n; j++)
		for (i = 0; i<n; i++)
			u[i][j] = 0.0;
}

//  ½«Êý×齵ά
//     | 0.0625  0.125  0.0625 |
//  R= | 0.125   0.25   0.125  |
//	   | 0.0625  0.125  0.0625 |
void rstrct(double **uc, double **uf, int nc, int nf)
{
	int ic, iif, jc, jf;

	for (jf = 2, jc = 1; jc<nc - 1; jc++, jf += 2) {
		for (iif = 2, ic = 1; ic<nc - 1; ic++, iif += 2) {
			uc[ic][jc] = 0.25*uf[iif][jf] + 0.125*(uf[iif + 1][jf] + uf[iif - 1][jf]
				+ uf[iif][jf + 1] + uf[iif][jf - 1]) + 0.0625*(uf[iif + 1][jf + 1] +
					uf[iif - 1][jf - 1] + uf[iif - 1][jf + 1] + uf[iif + 1][jf - 1]);
		}
	}
}

test.cpp

#include "stdafx.h"
#include <cv.h>
#include <highgui.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <afx.h>
#include <time.h>
#include "gross.h"
#include <string>
#include <opencv2\opencv.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\core\core.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
#define MAXITS 1000
#define EPS 1.0e-5
#define pi 3.1415
#define _AFXDLL
/* Function prototype */
void Histogram_equalizationg(IplImage *src);
void Histogram_specification(IplImage *src);
void Sobel(IplImage *src);
IplImage* LogarithmTransform(IplImage *src);
void LogEdge(IplImage *src);
void LogDCT(IplImage *src);
void mat2gray(CvArr* src);
void MQI(IplImage *src);
void DMQI(IplImage *src);
void SQI(IplImage *src);
void TT(IplImage *src);
void MultiScaleRetinex(IplImage *img);
void modification(IplImage *src);
void SSR(IplImage *img);
IplImage* DMQIL(IplImage *src);
void sort(double vector[], int size);
IplImage* MQIDenoise(IplImage *src);
IplImage** CreateMQIImage(int masksize, IplImage *src_denoise);
double** CreateGaussionFilter(int size, double sigma);

//std::string filename = "C:\\Users\\˼³¬\\Desktop\\yaleB07_P00A-020E+10_.jpg";
//std::string filename = "C:\\Users\\˼³¬\\Desktop\\1.bmp";
//CString dir = "C:\\Users\\˼³¬\\Desktop\\paper\\database\\new-CAS-PEAL\\Lighting";
//CString dir = "C:\\Users\\˼³¬\\Desktop\\paper\\database\\new-CAS-PEAL\\Lighting";
CString dir = "C:\\Users\\Stefan\\Desktop\\test";
CString dstdir = "C:\\Users\\Stefan\\Desktop\\res\\";
//CString dir = "C:\\Users\\˼³¬\\Desktop\\н¨Îļþ¼Ð (3)";
//CString dstdir = "C:\\Users\\˼³¬\\Desktop\\н¨Îļþ¼Ð (4)\\";
char* dstfilename;

/* main program */
int main(int argc, char* argv[])
{
	
	CFileFind finder;
	BOOL find = finder.FindFile(dir+"\\*.jpg");
	int i,len;
	CString path;	
	char* filepath;
	CString name;
	CString dstfile;

	int cnt=0;
	clock_t start, end;
	double dif;
	int ans;
	char ns;

	while (find)
	{
		find = finder.FindNextFile();
			
		path = finder.GetFilePath();
			name = finder.GetFileName();
			dstfile = dstdir + name;

			/*len = path.GetLength();
			filepath = new char[len];
			for (i = 0; i < len; i++)
				filepath[i] = (char)path.GetAt(i);//°ÑCStringת³Échar*
			*/
			len = WideCharToMultiByte(CP_ACP, 0, path, path.GetLength(), NULL, 0, NULL, NULL);
			filepath = new char[len + 1];
			WideCharToMultiByte(CP_ACP, 0, path, path.GetLength() + 1, filepath, len + 1, NULL, NULL);
			filepath[len] = '\0';
			

			/*len = dstfile.GetLength();
			dstfilename = new char[len];
			for (i = 0; i < len; i++)
				dstfilename[i] = (char)dstfile.GetAt(i);//°ÑCStringת³Échar*
			*/
			len = WideCharToMultiByte(CP_ACP, 0, dstfile, dstfile.GetLength(), NULL, 0, NULL, NULL);
			dstfilename = new char[len + 1];
			WideCharToMultiByte(CP_ACP, 0, dstfile, dstfile.GetLength() + 1, dstfilename, len + 1, NULL, NULL);
			dstfilename[len] = '\0';
			

			IplImage *src = 0; //¶¨ÒåԴͼÏñÖ¸Õë	
			IplImage *rst = 0;

			src = cvLoadImage(filepath, 0);
			/*int flag = 1;
			while (flag)
			{
				int ans;
				char ns;
				cout << " 1. Histogram equaliztion" << endl;
				cout << " 2. Histogram specification" << endl;
				cout << " 3. Gradient(Sobel)" << endl;
				cout << " 4. Logarithm transform" << endl;
				cout << " 5. Log edge" << endl;
				cout << " 6. LogDCT" << endl;
				cout << " 7. MQI" << endl;
				cout << " 8. DMQI" << endl;
				cout << " 9. SQI" << endl;
				cout << "10. Gross & Brajovic(sor)" << endl;
				cout << "11. Gross & Brajovic(Multigrid)" << endl;
				cout << "12. DMQI in log domain" << endl;


				cout << "Please choose algorithm:(1-12)";
				cin >> ans;

				start = clock();
				switch (ans)
				{
				case 1: Histogram_equalizationg(src); break;
				case 2: Histogram_specification(src); break;
				case 3: Sobel(src); break;
				case 4: LogarithmTransform(src); break;
				case 5: LogEdge(src); break;
				case 6: LogDCT(src); break;
				case 7: MQI(src); break;
				case 8: DMQI(src); break;
				case 9: SQI(src); break;
				case 10: GrossBrajovic(src); break;
				case 11: Gross(src, filename); break;
				case 12: DMQIL(src);
				default: break;
				}
				end = clock();
				dif = (double)(end - start) / CLOCKS_PER_SEC;
				printf("It took %f seconds.\n", dif);

				cvNamedWindow("src", CV_WINDOW_AUTOSIZE);
				cvShowImage("src", src);

				cvWaitKey(0);
				cout << "Continue?";
				cin >> ns;
				cout << endl;
				if (ns == 'N' || ns == 'n') flag = 0;

			}
			cvReleaseImage(&src);*/


			
			ans = 15;
			if(cnt==0)
			start = clock();
			switch (ans)
			{
			case 1: Histogram_equalizationg(src); break;
			case 2: Histogram_specification(src); break;
			case 3: Sobel(src); break;
			case 4: LogarithmTransform(src); break;
			case 5: LogEdge(src); break;

			case 6: LogDCT(src); break;
			case 7: MQI(src); break;
			case 8: DMQI(src); break;
			case 9: SQI(src); break;
			case 10: GrossBrajovic(src,dstfilename); break;
			case 11: Gross(src, filepath); break;
			case 12: DMQIL(src); break;
			case 13: TT(src); break;
			case 14: MultiScaleRetinex(src); break;
			case 15: modification(src); break;    //¸Ä½øµÄSQI
			case 16: SSR(src); break;
			default: break;
			}
			
			//cvNamedWindow("src", CV_WINDOW_AUTOSIZE);
			//cvShowImage("src", src);

			cvReleaseImage(&src);
			cnt++;
	}
	finder.Close();
	end = clock();
	dif = (double)(end - start) / CLOCKS_PER_SEC;
	printf("It took %f seconds.\n", dif);
	cvWaitKey(0);
	system("pause");
	return 0;
}

void modification(IplImage *src)
{
	int i, j, x, y, ii;
	int k = 0;
	int size = 31;
	int siz = (size - 1) / 2;
	double sigma = 8.0;//sigma=1.2
	double sum = 0;
	double tao = 0;
	double N = 0;
	double convlution_sum = 0;
	double** im_kernel = (double **)malloc(size * sizeof(double*));
	for (j = 0; j < size; j++) im_kernel[j] = (double *)malloc(size * sizeof(double));

	IplImage *src_guassion = 0; //¾­¹ý¸ß˹Â˲¨µÄͼÏñ
	IplImage *src_sqi = 0;//×ÔÉÌͼÏñ
	src_guassion = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	src_sqi = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);

	cvConvert(src, src_guassion);
	cvSmooth(src_guassion, src_guassion, CV_GAUSSIAN, 51, 51, 8.0);//¸ß˹º¯Êý¾í»ý²Ù×÷

	//»ñµÃ×ÔÉÌͼÏñ
	for (i = 0; i < src->height; i++)
	{
		for (j = 0; j < src->width; j++)
		{
			double q = atan(cvGet2D(src, i, j).val[0] / (cvGet2D(src_guassion, i, j).val[0] + 1));
			cvSet2D(src_sqi, i, j, cvScalar(q));
		}
	}
/*	for (i = 0; i < src_sqi->height; i++)
	{
		for (j = 0; j < src_sqi->width; j++)
		{
			if (cvGet2D(src_sqi3, i, j).val[0] > alpha * cvGet2D(src_sqi1, i, j).val[0])
			{
				cvSet2D(src_sqi, i, j, cvGet2D(src_sqi3, i, j));
			}
			else if (cvGet2D(src_sqi3, i, j).val[0]
				< alpha * cvGet2D(src_sqi1, i, j).val[0] &&
				cvGet2D(src_sqi3, i, j).val[0]
				> beta * cvGet2D(src_sqi1, i, j).val[0])
			{
				cvSet2D(src_sqi, i, j, cvGet2D(src_sqi2, i, j));
			}
			else if (cvGet2D(src_sqi3, i, j).val[0]
				< beta * cvGet2D(src_sqi1, i, j).val[0])
			{
				cvSet2D(src_sqi, i, j, cvGet2D(src_sqi1, i, j));
			}
		}
	}
	*/
	//cvNamedWindow("src_sqi", CV_WINDOW_AUTOSIZE);
	//cvShowImage("src_sqi", src_sqi);
	
	//»ñÈ¡¹âÕÕ·ÖÁ¿
	cvSmooth(src, src_guassion, CV_GAUSSIAN, 51, 51);
	IplImage *log_guassion, *log_sqi;
	log_guassion = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);
	log_sqi = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);

	//¶Ô¹âÕÕ·ÖÁ¿½øÐÐÙ¤ÂíУÕý
	Mat m_FaceImg = cvarrToMat(src_guassion);
	m_FaceImg.convertTo(m_FaceImg, CV_32FC1, 1.0 / 255);
	pow(m_FaceImg, 0.1, m_FaceImg);
	log_guassion = &IplImage(m_FaceImg);
	
	cvMul(log_guassion, src_sqi, src_sqi);//²¹¹â²Ù×÷
	
	IplImage *src_sqi_8u = cvCreateImage(cvGetSize(src_sqi), IPL_DEPTH_8U, src->nChannels);
	cvConvertScaleAbs(src_sqi, src_sqi_8u, 255, 0);
	cvSaveImage(dstfilename, src_sqi_8u);

	for (j = 0; j < size; j++) free(im_kernel[j]);
	free(im_kernel);
	cvReleaseImage(&src_guassion);
}

void SSR(IplImage *img)
{
	IplImage *A, *fA, *fB, *fC;
	double sigma = 8.0;
	double filter_size = 51;
	// Initialize temp images
	fA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);
	fB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);
	fC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);

	// Compute log image
	cvConvert(img, fA);
	cvLog(fA, fB);

	// Compute log of blured image
	A = cvCloneImage(img);
	//FastFilter(A, sigma);
	cvSmooth(A, A, CV_GAUSSIAN, filter_size, filter_size, sigma);
	cvConvert(A, fA);
	cvLog(fA, fC);

	// Compute difference
	cvSub(fB, fC, fA);

	// Restore
	cvExp(fA, fA);
	
	
	IplImage *src_guassion = 0;
	src_guassion = cvCreateImage(cvGetSize(img), IPL_DEPTH_8U, img->nChannels);
	cvSmooth(img, src_guassion, CV_GAUSSIAN, 31, 31);
	IplImage *log_guassion, *log_ssr;
	log_guassion = cvCreateImage(cvGetSize(img), IPL_DEPTH_32F, img->nChannels);
	log_ssr = cvCreateImage(cvGetSize(img), IPL_DEPTH_32F, img->nChannels);

	Mat m_FaceImg = cvarrToMat(src_guassion);
	m_FaceImg.convertTo(m_FaceImg, CV_32FC1, 1.0 / 255);
	pow(m_FaceImg, 0.5, m_FaceImg);
	log_guassion = &IplImage(m_FaceImg);

	cvNamedWindow("log_guassion", CV_WINDOW_AUTOSIZE);
	cvShowImage("log_guassion", log_guassion);
	cvMul(log_guassion, fA, fA);

	IplImage *src_ssr_8u = cvCreateImage(cvGetSize(fA), IPL_DEPTH_8U, fA->nChannels);
	cvConvertScaleAbs(fA, src_ssr_8u, 255, 0);
	cvSaveImage(dstfilename, src_ssr_8u);

	// Release temp images
	cvReleaseImage(&A);
	cvReleaseImage(&fA);
	cvReleaseImage(&fB);
	cvReleaseImage(&fC);
}

void MultiScaleRetinex(IplImage *img)
{
	int scales = 3;
	int filter_size = 31;
	double weight[3];
	double sigma[3];
	sigma[0] = 3; sigma[1] = 4; sigma[2] = 5;
	weight[0] = (1 * 1.0) / 3; weight[1] = (1 * 1.0) / 3; weight[2] = (1 * 1.0) / 3;
	//int gain = 128;
	//int offset = 128;

	int i;
	int j;
	//double weight;
	IplImage *A, *fA, *fB, *fC, *fD, *fE;

	// Initialize temp images
	fA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);
	fB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);
	fC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);
	fD = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);
	fE = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);


	// Compute log image
	cvConvert(img, fA);
	cvLog(fA, fB);

	// Normalize according to given weights
	//for (i = 0, weight = 0; i < scales; i++)
	//	weight += weights[i];

	//if (weight != 1.0) cvScale(fB, fB, weight);

	// Filter at each scale
	//ÒòΪûÓгõʼ»¯fD£¬ËùÒ԰ѵÚÒ»¸öÑ­»·µ¥¶ÀÄóöÀ´
	A = cvCloneImage(img);
	//FastFilter(A, sigmas[0]);
	cvSmooth(A, A, CV_GAUSSIAN, filter_size, filter_size, sigma[0]);

	cvConvert(A, fA);
	cvLog(fA, fC);
	//cvReleaseImage(&A);

	// Compute weighted difference
	//cvScale(fC, fC, weights[i]);
	cvSub(fB, fC, fD);
	cvScale(fD, fD, weight[0]);
	cvConvert(fD, fE);

	for (i = 1; i < scales; i++) {
		A = cvCloneImage(img);
		//FastFilter(A, sigmas[i]);
		cvSmooth(A, A, CV_GAUSSIAN, filter_size, filter_size, sigma[i]);

		cvConvert(A, fA);
		cvLog(fA, fC);
		//cvReleaseImage(&A);

		// Compute weighted difference
		//cvScale(fC, fC, weights[i]);
		cvSub(fB, fC, fD);
		cvScale(fD, fD, weight[i]);
		cvAdd(fD, fE, fE);
	}

	// Restore
	cvExp(fE, fE);

	for(i=0;i<fE->height;i++)
		for (j = 0; j < fE->width; j++)
		{
			CvScalar value;
			value.val[0] = atan(cvGet2D(fE, i, j).val[0]);
			cvSet2D(fE, i, j, value);
		}
	
	//cvConvertScale(fE, img, gain, offset);

	IplImage *src_msr_8u = cvCreateImage(cvGetSize(fE), IPL_DEPTH_8U, fE->nChannels);
	cvConvertScaleAbs(fE, src_msr_8u, 255, 0);
	cvSaveImage(dstfilename, src_msr_8u);
	// Release temp images
	cvReleaseImage(&A);
	cvReleaseImage(&fA);
	cvReleaseImage(&fB);
	cvReleaseImage(&fC);
}


void TT(IplImage *src)
{
	    Mat m_FaceImg=cvarrToMat(src);
		Mat m_fFaceImg, m_fFaceImg1, m_fFaceImg2;
		double sigma0 = 0.01, sigma1 = 30, gamma = 0.1, alpha = 0.01, tau = 8;
		m_FaceImg.convertTo(m_fFaceImg, CV_32FC1, 1.0 / 255);//8u->32f

		//gamma correction
		pow(m_fFaceImg, gamma, m_fFaceImg);

		//DoG filter
		Mat dst1, dst2;
		int dia = 51;
		GaussianBlur(m_fFaceImg, dst1, Size(dia, dia), sigma0, sigma0);//sigma0ԽС±£ÁôµÄ¸ßƵÐźÅÔ½¶à
		GaussianBlur(m_fFaceImg, dst2, Size(dia, dia), sigma1, sigma1);//sigma1Ô½´ó±£ÁôµÄµÍƵÐźÅÔ½¶à
		m_fFaceImg = dst1 - dst2;
		
		//contrast equalization
		//img=img/(mean2(abs(img).^alpha)^(1/alpha))
		m_fFaceImg1 = abs(m_fFaceImg);
		pow(m_fFaceImg1, alpha, m_fFaceImg1);
		m_fFaceImg /= pow(mean(m_fFaceImg1).val[0], 1 / alpha);

		//img=img/(mean2(min(tau,abs(img)).^alpha)^(1/alpha));
		m_fFaceImg1 = abs(m_fFaceImg);
		m_fFaceImg1 = (cv::min)(tau, m_fFaceImg1);
		pow(m_fFaceImg1, alpha, m_fFaceImg1);
		m_fFaceImg = m_fFaceImg / pow(mean(m_fFaceImg1).val[0], 1 / alpha) / tau;

		//tanh x=(e^(x)-e^(-x))/(e^x+e^(-x))
		exp(m_fFaceImg, m_fFaceImg1);
		exp((0 - m_fFaceImg), m_fFaceImg2);
		m_fFaceImg = m_fFaceImg1 - m_fFaceImg2;
		m_fFaceImg1 = m_fFaceImg1 + m_fFaceImg2;
		m_fFaceImg = m_fFaceImg / m_fFaceImg1;

		normalize(m_fFaceImg, m_fFaceImg, 0, 1, CV_MINMAX);//¹éÒ»»¯
		convertScaleAbs(m_fFaceImg, m_FaceImg, 255);//32fת³É8u

		//m_FaceImg -= 10;
		//namedWindow("TT", WINDOW_AUTOSIZE);
		//imshow("TT", m_FaceImg);
		IplImage *tt = &IplImage(m_FaceImg);
		cvSaveImage(dstfilename, tt);
}

//Histogram equaliztion
void Histogram_equalizationg(IplImage *src)
{
	IplImage *src_hist = 0;
	src_hist = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	cvEqualizeHist(src, src_hist);
	//cvNamedWindow("src_HE", CV_WINDOW_AUTOSIZE);
	//cvShowImage("src_HE", src_hist);
	//std::string name = "HE_" + filename;
	cvSaveImage(dstfilename, src_hist);
}

//Histogram specification
void Histogram_specification(IplImage *src)
{
	IplImage *src_ref = cvLoadImage("2.tif", 0);
	IplImage *dst = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);

	CvHistogram* hist_src = 0;
	CvHistogram* hist_ref = 0;

	CvMat* lut_src = 0;
	CvMat* lut_ref = 0;

	int i, hist_sz = 256;
	CvSize img_sz;
	float scale;
	float* h;
	int sum = 0;
	hist_src = cvCreateHist(1, &hist_sz, CV_HIST_ARRAY);
	hist_ref = cvCreateHist(1, &hist_sz, CV_HIST_ARRAY);
	lut_src = cvCreateMat(1, 256, CV_8UC1);
	lut_ref = cvCreateMat(1, 256, CV_8UC1);

	cvCalcArrHist((CvArr**)&src, hist_src);
	cvCalcArrHist((CvArr**)&src_ref, hist_ref);

	img_sz = cvGetSize(src);
	scale = 255.f / (img_sz.width*img_sz.height);

	h = (float*)cvPtr1D(hist_src->bins, 0);
	for (i = 0; i < hist_sz; i++)
	{
		sum += cvRound(h[i]);
		lut_src->data.ptr[i] = (uchar)cvRound(sum*scale);
	}
	lut_src->data.ptr[0] = 0;

	h = (float*)cvPtr1D(hist_ref->bins, 0);
	for (i = 0; i < hist_sz; i++)
	{
		sum += cvRound(h[i]);
		lut_ref->data.ptr[i] = (uchar)cvRound(sum*scale);
	}
	lut_ref->data.ptr[0] = 0;

	int i1, i2;
	CvMat* lut = cvCreateMat(1, 256, CV_8UC1);
	for (i1 = 0; i1<256; i1++)
	{
		for (i2 = 0; i2<256; i2++)
		{
			if (lut_src->data.ptr[i1] >= lut_ref->data.ptr[i2])
			{
				lut->data.ptr[i1] = i2;
			}
		}
	}

	cvLUT(src, dst, lut);

	//cvNamedWindow("src_HS", CV_WINDOW_AUTOSIZE);
	//cvShowImage("src_HS", dst);
	//std::string name = "HS_" + filename;
	cvSaveImage(dstfilename, dst);

	cvReleaseHist(&hist_src);
	cvReleaseHist(&hist_ref);
	cvReleaseMat(&lut);
	cvReleaseMat(&lut_src);
	cvReleaseMat(&lut_ref);
}

//Sobel
void Sobel(IplImage *src)
{
	IplImage *src_sobel_x = 0;
	IplImage *src_sobel_y = 0;
	IplImage *src_sobel1 = 0;
	IplImage *src_sobel2 = 0;
	IplImage *src_sobel_8u = 0;
	IplImage *src_sobel_8u_max = 0;
	IplImage *src_sobel_Q = 0;

	//ÓÃÓÚ½ÓÊܾ­¹ýÌݶȼÆËãµÄͼ£¬cvSobelÒªÇóÊä³öͼÖÁÉÙΪ16λ
	src_sobel_x = cvCreateImage(cvGetSize(src), IPL_DEPTH_16S, src->nChannels);
	src_sobel_y = cvCreateImage(cvGetSize(src), IPL_DEPTH_16S, src->nChannels);
	src_sobel1 = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	src_sobel2 = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	//´Ó16λתΪ8λµÄ½á¹û
	src_sobel_8u = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	//ÉÌͼ
	src_sobel_Q = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);

	cvSobel(src, src_sobel_x, 1, 0, 3);
	cvConvertScaleAbs(src_sobel_x, src_sobel1, 1, 0);
	cvSobel(src, src_sobel_y, 0, 1, 3);
	cvConvertScaleAbs(src_sobel_y, src_sobel2, 1, 0);
	//Á½¸öËã×ÓÏà¼Ó
	cvAdd(src_sobel1, src_sobel2, src_sobel_8u);
	//cvNamedWindow( "src_sobel_8u", CV_WINDOW_AUTOSIZE );
	//cvShowImage( "src_sobel_8u", src_sobel_8u );

	src_sobel_8u_max = cvCreateImage(cvSize(src->width, src->height), IPL_DEPTH_8U, src->nChannels);
	cvDilate(src_sobel_8u, src_sobel_8u_max, NULL, 1);

	for (int i = 0; i < src->height; i++)
	{
		for (int j = 0; j < src->width; j++)
		{
			double q = cvGet2D(src_sobel_8u, i, j).val[0] / (cvGet2D(src_sobel_8u_max, i, j).val[0] + 1);
			cvSet2D(src_sobel_Q, i, j, cvScalar(q));
		}
	}
	cvConvertScaleAbs(src_sobel_Q, src_sobel_8u, 255, 0);

	//cvNamedWindow( "src_sobel_8u_max", CV_WINDOW_AUTOSIZE );
	//cvShowImage( "src_sobel_8u_max", src_sobel_8u_max );
	cvNamedWindow("src_sobel_Q", CV_WINDOW_AUTOSIZE);
	cvShowImage("src_sobel_Q", src_sobel_8u);
	//std::string name = "Sobel_" + filename;
	//cvSaveImage(name.c_str(), src_sobel_8u);

	cvReleaseImage(&src_sobel_x);
	cvReleaseImage(&src_sobel_y);
	cvReleaseImage(&src_sobel1);
	cvReleaseImage(&src_sobel2);
	cvReleaseImage(&src_sobel_8u_max);
	cvReleaseImage(&src_sobel_Q);
}

void sort(double vector[], int size)
{
	double t;
	for (int i = 0; i<size; i++)
	{
		for (int j = i + 1; j<size; j++)
		{
			if (vector[i] < vector[j])
			{
				t = vector[i];
				vector[i] = vector[j];
				vector[j] = t;
			}
		}
	}
}

//Logarithm Transform
IplImage* LogarithmTransform(IplImage *src)
{
	IplImage *src_log = 0;
	src_log = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);

	for (int i = 0; i < src_log->height; i++)
	{
		for (int j = 0; j < src_log->width; j++)
		{
			double q = log(cvGet2D(src, i, j).val[0] + 1);
			cvSet2D(src_log, i, j, cvScalar(q));
		}
	}

	mat2gray(src_log);
	//cvNamedWindow("src_log", CV_WINDOW_AUTOSIZE);
	//cvShowImage("src_log", src_log);

	IplImage *src_log_8u = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	cvConvertScaleAbs(src_log, src_log_8u, 255, 0);
	//std::string name = "log_" + filename;
	//cvSaveImage(name.c_str(), src_log_8u);
	cvSaveImage(dstfilename, src_log);

	return src_log;
}

//Log edge
void LogEdge(IplImage *src)
{
	IplImage *src_log = 0;
	IplImage *src_log_edge = 0;
	double min = 0;
	int ii, jj;

	src_log = LogarithmTransform(src);
	src_log_edge = cvCreateImage(cvSize(src->width - 4, src->height - 4), src_log->depth, src->nChannels);

	for (int i = 2; i < src->height - 2; i++)
	{
		for (int j = 2; j < src->width - 2; j++)
		{

			min = cvGet2D(src, i - 2, j - 2).val[0];
			for (ii = -2; ii <= 2; ii++)
			{
				for (jj = -2; jj <= 2; jj++)
				{
					if (cvGet2D(src_log, i + ii, j + jj).val[0]<min)
						min = cvGet2D(src_log, i + ii, j + jj).val[0];
				}
			}
			double change = cvGet2D(src_log, i, j).val[0] - min;
			cvSet2D(src_log_edge, i - 2, j - 2, cvScalar(change));
		}
	}
	cvNamedWindow("src_log_edge", CV_WINDOW_AUTOSIZE);
	cvShowImage("src_log_edge", src_log_edge);

	IplImage *src_log_edge_8u = cvCreateImage(cvGetSize(src_log_edge), IPL_DEPTH_8U, src->nChannels);
	cvConvertScaleAbs(src_log_edge, src_log_edge_8u, 255, 0);
	//std::string name = "logedge_" + filename;
	//cvSaveImage(name.c_str(), src_log_edge_8u);

	cvReleaseImage(&src_log);
}

//LogDCT
void LogDCT(IplImage *src)
{
	int i, j;
	IplImage *src_log = 0;
	IplImage *src_log_double = 0;
	IplImage *src_logDCT_origin = 0;
	IplImage *src_logDCT = 0;

	src_log = LogarithmTransform(src);
	src_logDCT_origin = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);
	src_logDCT = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);

	cvDCT(src_log, src_logDCT_origin, CV_DXT_FORWARD);

	int Ddis = 20;
	double u = 100;
	CvScalar s;

	for (i = 0; i < Ddis + 1; i++)
	{
		for (j = 0; j < Ddis + 1; j++)
		{
			s = cvGet2D(src_logDCT_origin, i, j);
			s.val[0] = s.val[0]/50;
			if ((i + j) <= (Ddis + 1))	cvSet2D(src_logDCT_origin, i, j, s);
		}
	}
	s= cvGet2D(src_logDCT_origin, 0, 0);
	s.val[0] *= 1.1;
	cvSet2D(src_logDCT_origin, 0, 0, s);

	double avg = log(u)* sqrt(src->height * src->width);
	cvSet2D(src_logDCT_origin, 0, 0, cvScalar(avg));

	cvDCT(src_logDCT_origin, src_logDCT, CV_DXT_INVERSE);
	mat2gray(src_logDCT);

	cvNamedWindow("src_logDCT", CV_WINDOW_AUTOSIZE);
	cvShowImage("src_logDCT", src_logDCT);

	IplImage *src_logDCT_8u = cvCreateImage(cvGetSize(src_logDCT), IPL_DEPTH_8U, src->nChannels);
	cvConvertScaleAbs(src_logDCT, src_logDCT_8u, 255, 0);
	//std::string name = "logDCT_" + filename;
	//cvSaveImage(name.c_str(), src_logDCT_8u);
}

void mat2gray(CvArr* src)
{
	int i, j;
	int height = cvGetSize(src).height;
	int width = cvGetSize(src).width;
	double max, min;

	max = min = cvGet2D(src, 0, 0).val[0];
	for (i = 0; i < height; i++)
	{
		for (j = 0; j < width; j++)
		{
			if (cvGet2D(src, i, j).val[0]>max) max = cvGet2D(src, i, j).val[0];
			else if (cvGet2D(src, i, j).val[0]<min) min = cvGet2D(src, i, j).val[0];
		}
	}
	double dif = max - min;
	for (i = 0; i < height; i++)
	{
		for (j = 0; j < width; j++)
		{
			double q = (cvGet2D(src, i, j).val[0] - min) / dif;
			cvSet2D(src, i, j, cvScalar(q));
		}
	}
}

//MQI
void MQI(IplImage *src)
{
	IplImage *src_denoise = 0;
	IplImage **src_MQI = 0;
	//int masksize = 11;
	int masksize = 31;

	src_denoise = src_denoise = MQIDenoise(src);

	src_MQI = CreateMQIImage(masksize, src_denoise);

	IplImage *src_MQI_8u = cvCreateImage(cvGetSize(src_MQI[0]), IPL_DEPTH_8U, src->nChannels);
	cvConvertScaleAbs(src_MQI[0], src_MQI_8u, 255, 0);

	//cvNamedWindow("src_MQI", CV_WINDOW_AUTOSIZE);
	//cvShowImage("src_MQI", src_MQI_8u);

	//std::string name = "MQI_" + filename;
	cvSaveImage(dstfilename, src_MQI_8u);

	cvReleaseImage(&src_denoise);
	cvReleaseImage(&src_MQI[0]);
	cvReleaseImage(&src_MQI[1]);
}

//DMQI
void DMQI(IplImage *src)
{
	IplImage *src_denoise = 0;
	IplImage **src_MQI[3] = { 0,0,0 };
	IplImage *src_DMQI = 0;
	int i, j;
	int mqi = 0, close = 1;
	int masksize = 7;
	int step = 2;
	int masknumber = 3;
	double alpha = 1.8, beta = 1.35;

	src_denoise = MQIDenoise(src);
	src_DMQI = cvCreateImage(cvGetSize(src_denoise), IPL_DEPTH_32F, src_denoise->nChannels);;

	for (i = 0; i<masknumber; i++)
	{
		src_MQI[i] = CreateMQIImage(masksize, src_denoise);
		masksize += step;
	}

	//cvNamedWindow( "src_MQI_5", CV_WINDOW_AUTOSIZE );
	//cvShowImage( "src_MQI_5", src_MQI[0][mqi] );	
	//cvNamedWindow( "src_MQI_7", CV_WINDOW_AUTOSIZE );
	//cvShowImage( "src_MQI_7", src_MQI[1][mqi] );	
	//cvNamedWindow( "src_MQI_9", CV_WINDOW_AUTOSIZE );
	//cvShowImage( "src_MQI_9", src_MQI[2][mqi] );

	for (i = 0; i < src_denoise->height; i++)
	{
		for (j = 0; j < src_denoise->width; j++)
		{
			if (cvGet2D(src_MQI[2][close], i, j).val[0] > alpha * cvGet2D(src_MQI[0][close], i, j).val[0])
			{
				cvSet2D(src_DMQI, i, j, cvGet2D(src_MQI[2][mqi], i, j));
			}
			else if (cvGet2D(src_MQI[2][close], i, j).val[0]
				< alpha * cvGet2D(src_MQI[0][close], i, j).val[0] &&
				cvGet2D(src_MQI[2][close], i, j).val[0]
				> beta * cvGet2D(src_MQI[0][close], i, j).val[0])
			{
				cvSet2D(src_DMQI, i, j, cvGet2D(src_MQI[1][mqi], i, j));
			}
			else if (cvGet2D(src_MQI[2][close], i, j).val[0]
				< beta * cvGet2D(src_MQI[0][close], i, j).val[0])
			{
				cvSet2D(src_DMQI, i, j, cvGet2D(src_MQI[0][mqi], i, j));
			}
		}
	}

	//IplImage *src_DMQI_out = 0;
	//src_DMQI_out=cvCreateImage( cvGetSize(src_denoise), IPL_DEPTH_8U, src->nChannels);
	//cvConvertScaleAbs( src_DMQI, src_DMQI_out, 180, 0 );

	cvNamedWindow("src_DMQI", CV_WINDOW_AUTOSIZE);
	cvShowImage("src_DMQI", src_DMQI);

	IplImage *src_DMQI_8u = cvCreateImage(cvGetSize(src_DMQI), IPL_DEPTH_8U, src->nChannels);
	cvConvertScaleAbs(src_DMQI, src_DMQI_8u, 255, 0);
	//std::string name = "DMQI_" + filename;
	//cvSaveImage(name.c_str(), src_DMQI_8u);

	cvReleaseImage(&src_denoise);
	for (i = 0; i<3; i++) {
		for (j = 0; j<2; j++) {
			cvReleaseImage(&src_MQI[i][j]);
		}
	}
	cvReleaseImage(&src_DMQI);
}

IplImage* MQIDenoise(IplImage *src)
{
	int i, j;
	IplImage *src_denoise;
	src_denoise = cvCreateImage(cvSize(src->width - 2, src->height - 2), src->depth, src->nChannels);

	for (i = 1; i < src->height - 1; i++) //denoise
	{
		for (j = 1; j < src->width - 1; j++)
		{
			double vector[8] = { 0 };
			vector[0] = cvGet2D(src, i - 1, j - 1).val[0];
			vector[1] = cvGet2D(src, i - 1, j).val[0];
			vector[2] = cvGet2D(src, i - 1, j + 1).val[0];
			vector[3] = cvGet2D(src, i, j - 1).val[0];
			vector[4] = cvGet2D(src, i, j + 1).val[0];
			vector[5] = cvGet2D(src, i + 1, j - 1).val[0];
			vector[6] = cvGet2D(src, i + 1, j).val[0];
			vector[7] = cvGet2D(src, i + 1, j + 1).val[0];
			sort(vector, 8);
			if (cvGet2D(src, i, j).val[0] > vector[0])
			{
				cvSet2D(src_denoise, i - 1, j - 1, cvScalar(vector[0]));
			}
			else if (cvGet2D(src, i, j).val[0] < vector[7])
			{
				cvSet2D(src_denoise, i - 1, j - 1, cvScalar(vector[7]));
			}
			else
			{
				cvSet2D(src_denoise, i - 1, j - 1, cvScalar(cvGet2D(src, i, j).val[0]));
			}
		}
	}
	return src_denoise;
}

/* Funtion return two images. The first one is MQI image and
* the second one is the result of close operation.
*/
IplImage** CreateMQIImage(int masksize, IplImage *src_denoise)
{
	IplImage *temp = 0;
	IplImage *src_denoise_q = 0;
	IplImage *src_denoise_c = 0;
	IplImage **rn = (IplImage**)malloc(2 * sizeof(IplImage*));

	src_denoise_c = cvCreateImage(cvGetSize(src_denoise), src_denoise->depth, src_denoise->nChannels);
	src_denoise_q = cvCreateImage(cvGetSize(src_denoise), IPL_DEPTH_32F, src_denoise->nChannels);

	IplConvKernel *Kernel;
	int nCols = masksize;
	int nRows = masksize;
	int anchorX = (masksize - 1) / 2;
	int anchorY = (masksize - 1) / 2;
	Kernel = cvCreateStructuringElementEx(nCols, nRows, anchorX, anchorY, CV_SHAPE_RECT, NULL);
	cvMorphologyEx(src_denoise, src_denoise_c, temp, Kernel, CV_MOP_CLOSE, 1);

	for (int i = 0; i < src_denoise->height; i++)
	{
		for (int j = 0; j < src_denoise->width; j++)
		{
			double q = cvGet2D(src_denoise, i, j).val[0] / (cvGet2D(src_denoise_c, i, j).val[0] + 1);
			cvSet2D(src_denoise_q, i, j, cvScalar(q));
		}
	}

	cvReleaseStructuringElement(&Kernel);

	rn[0] = src_denoise_q;
	rn[1] = src_denoise_c;
	return rn;
}

//SQI
void SQI(IplImage *src)
{

	/*Mat m_src = cvarrToMat(src);
	m_src.convertTo(m_src, CV_32FC1, 1.0 / 255);
	pow(m_src, 0.8, m_src);
	m_src.convertTo(m_src, CV_8U, 255);
	src = &IplImage(m_src);*/

	int i, j, x, y, ii;
	int k = 0;
	int size = 31;
	int siz = (size - 1) / 2;
	double sigma = 8.0;//sigma=1.2
	double sum = 0;
	double tao = 0;
	double N = 0;
	double convlution_sum = 0;
	double** im_kernel = (double **)malloc(size * sizeof(double*));
	for (j = 0; j < size; j++) im_kernel[j] = (double *)malloc(size * sizeof(double));

	IplImage *src_guassion = 0; //¾­¹ý¸ß˹Â˲¨µÄͼÏñ
	IplImage *src_sqi = 0;//×ÔÉÌͼÏñ
	src_guassion = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	src_sqi = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);

	for (i = 0; i<src->height; i++)
	{
		for (j = 0; j<src->width; j++)
		{
			double** kernel = CreateGaussionFilter(size, sigma);//Éú³É¸ß˹ºËº¯Êý
			//»ñÈ¡ÔÚ¾í»ý´°¿ÚÀïÃæµÄËùÓеãµÄ»Ò¶ÈÖµ£¬±£´æÔÚim_kernelÖÐ
			for (x = i - siz; x <= i + siz; x++)
			{
				for (y = j - siz; y <= j + siz; y++)
				{
					if (x<0 || x >= src->height || y<0 || y >= src->width)
					{
						im_kernel[x - i + siz][y - j + siz] = 0;
					}
					else
					{
						im_kernel[x - i + siz][y - j + siz] = cvGet2D(src, x, y).val[0];
					}
				}
			}

			//¼ÆËãtaoµÄÖµ£¬taoÊǾí»ý´°¿ÚÖеÄÏñËصã¾ùÖµ
			sum = 0;
			tao = 0;
			for (x = 0; x<size; x++)
			{
				for (y = 0; y<size; y++)
				{
					sum = sum + im_kernel[x][y];
				}
			}
			tao = sum / (size*size);

			//¼ÆËãȨÖØ£¬È¨ÖØÓÉÏñËصãµÄ»Ò¶ÈÖµºÍtaoµÄ´óС¹ØϵÀ´¾ö¶¨
			for (x = 0; x<size; x++)
			{
				for (y = 0; y<size; y++)
				{
					if (im_kernel[x][y] < tao) kernel[x][y] = 0;
				}
			}

			//¼ÆËã¹éÒ»»¯Òò×ÓN
			N = 0;
			for (x = 0; x<size; x++)
			{
				for (y = 0; y<size; y++)
				{
					N = N + kernel[x][y];
				}
			}

			//¾í»ý²Ù×÷
			convlution_sum = 0;
			for (x = 0; x<size; x++)
			{
				for (y = 0; y<size; y++)
				{
					//convlution
					//convlution_sum = convlution_sum + kernel[x][y] * im_kernel[size-x-1][size-y-1] / N;

					//correlation
					convlution_sum += kernel[x][y] * im_kernel[x][y] / N;//Ö´Ðоí»ý²Ù×÷
				}
			}

			cvSet2D(src_guassion, i, j, cvScalar(convlution_sum));
			for (ii = 0; ii < size; ii++) free(kernel[ii]);
			free(kernel);
		}
	}

	//cvNamedWindow( "src_guassion", CV_WINDOW_AUTOSIZE );
	//cvShowImage( "src_guassion", src_guassion );

	for (i = 0; i < src->height; i++)
	{
		for (j = 0; j < src->width; j++)
		{
			//»ñÈ¡·´Éä·ÖÁ¿²¢ÇÒ½«·´Éä·ÖÁ¿Í¨¹ý·´ÕýÇк¯Êý½øÐд¦Àí
			double q = atan(cvGet2D(src, i, j).val[0] / (cvGet2D(src_guassion, i, j).val[0] + 1));
			cvSet2D(src_sqi, i, j, cvScalar(q));
		}
	}
	

	//cvNamedWindow("src_sqi", CV_WINDOW_AUTOSIZE);
	//cvShowImage("src_sqi", src_sqi);
	
	/*cvSmooth(src, src_guassion, CV_GAUSSIAN, 15, 15);
	IplImage *log_guassion, *log_sqi;
	log_guassion = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);
	log_sqi = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);
	
	   /*for (i = 0; i<src_guassion->height; i++)
		   for (j = 0; j < src_guassion->width; j++)
		   {
			   cvSet2D(src_guassion, i, j, 60);
		   }*/

	/*Mat m_FaceImg = cvarrToMat(src_guassion);
	m_FaceImg.convertTo(m_FaceImg, CV_32FC1, 1.0 / 255);
	pow(m_FaceImg, 0.2, m_FaceImg);
	log_guassion = &IplImage(m_FaceImg);
	

	   /*for(i=0;i<log_guassion->height;i++)
		   for (j = 0; j < log_guassion->width; j++)
		   {
			   cvSet2D(log_guassion, i, j, cvGet2D(log_guassion, i, j).val[0]*(4.0/5));
		   }*/
	//cvNamedWindow("log_guassion", CV_WINDOW_AUTOSIZE);
	//cvShowImage("log_guassion", log_guassion);

	/*cvMul(log_guassion, src_sqi, src_sqi);
	
	   /*
	   cvLog(log_guassion, log_guassion);
	   cvLog(src_sqi, log_sqi);
	//cvAddWeighted(log_guassion, 0.7, log_sqi, 0.3, 0.0, log_sqi);
	//cvAdd(log_sqi, log_guassion,log_sqi);
	for (i = 0; i<log_guassion->height; i++)
		for (j = 0; j < log_guassion->width; j++)
		{
			cvSet2D(log_sqi, i, j, cvGet2D(log_guassion, i, j).val[0] + cvGet2D(log_sqi, i, j).val[0]);
		}
	cvExp(log_sqi, src_sqi);*/
	/*
	Mat m_FaceImg1 = cvarrToMat(src_sqi);
	pow(m_FaceImg1, 0.9, m_FaceImg1);
	src_sqi = &IplImage(m_FaceImg1);
	*/

	

	IplImage *src_sqi_8u = cvCreateImage(cvGetSize(src_sqi), IPL_DEPTH_8U, src->nChannels);
	cvConvertScaleAbs(src_sqi, src_sqi_8u, 255, 0);
	cvSaveImage(dstfilename, src_sqi_8u);

	for (j = 0; j < size; j++) free(im_kernel[j]);
	free(im_kernel);
	cvReleaseImage(&src_guassion);
//	cvReleaseImage(&src_sqi);


}

double** CreateGaussionFilter(int size, double sigma)
{
	int i, j;
	int siz = (size - 1) / 2;
	double sum = 0;

	double** kernel = (double**)malloc(size * sizeof(double *));
	for (i = 0; i < size; i++) kernel[i] = (double *)malloc(size * sizeof(double));

	for (i = -siz; i <= siz; i++)	// create the mask for gaussian filter
	{
		for (j = -siz; j <= siz; j++)
		{
			kernel[i + siz][j + siz] = exp((-(pow(i, 2) + pow(j, 2)) / (2 * sigma*sigma)));
		}
	}

	for (i = 0; i<size; i++)
	{
		for (j = 0; j<size; j++)
		{
			sum += kernel[i][j];
		}
	}

	for (i = 0; i<size; i++)
	{
		for (j = 0; j<size; j++)
		{
			kernel[i][j] = kernel[i][j] / sum;
		}
	}
	return kernel;
}

//Morphological Quotient Image in Logarithm Domain
IplImage* DMQIL(IplImage *src)
{
	int i, j;

	// Step 2: take the logarithm transformation.
	IplImage *src_log = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	cvConvertScaleAbs(LogarithmTransform(src), src_log, 255, 0);

	// Step 3: preprocess src_log to reduce the generated noise.
	IplImage *src_GRADIENT = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	IplImage *src_dilate = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	IplImage *src_denoise = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);

	IplImage *temp = cvCloneImage(src_log);
	int b1_masksize = 9;
	IplConvKernel *Kernel;
	int nCols = b1_masksize;
	int nRows = b1_masksize;
	int anchorX = (b1_masksize - 1) / 2;
	int anchorY = (b1_masksize - 1) / 2;
	Kernel = cvCreateStructuringElementEx(nCols, nRows, anchorX, anchorY, CV_SHAPE_RECT, NULL);
	cvMorphologyEx(src_log, src_GRADIENT, temp, Kernel, CV_MOP_GRADIENT, 1);
	cvDilate(src_log, src_dilate, Kernel, 1);

	int sigma = 10;
	for (i = 0; i<src->height; i++)
	{
		for (j = 0; j<src->width; j++)
		{
			if (cvGet2D(src_GRADIENT, i, j).val[0] < sigma)
			{
				cvSet2D(src_denoise, i, j, cvGet2D(src_dilate, i, j));
			}
			else
			{
				cvSet2D(src_denoise, i, j, cvGet2D(src_log, i, j));
			}
		}
	}

	// Step 4: conduct closing operation on src_denoise, and moreover,
	//         the structure elements can be chosen dynamically.
	IplImage *src_mgratio = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);
	temp = cvCloneImage(src_denoise);
	//¸üÐÂÐÎ̬ѧÌݶÈ
	cvMorphologyEx(src_denoise, src_GRADIENT, temp, Kernel, CV_MOP_GRADIENT, 1);

	//µÃµ½mgratio
	for (i = 0; i<src->height; i++)
	{
		for (j = 0; j<src->width; j++)
		{
			double ratio = cvGet2D(src_GRADIENT, i, j).val[0] / cvGet2D(src_denoise, i, j).val[0];
			cvSet2D(src_mgratio, i, j, cvScalar(ratio));
		}
	}

	IplImage* src_close[3] = { 0,0,0 };
	int s = 0, m = 1, l = 2;  //¶ÔÓ¦´óÖÐСÈý¸ö½á¹¹ÔªËØ
	int masksize = 15;
	int step = 2;
	int masknumber = 3;
	double alpha = 0.8, beta = 2.1;

	for (i = 0; i<masknumber; i++)
	{
		src_close[i] = (CreateMQIImage(masksize, src_denoise))[1];
		masksize += step;
	}

	IplImage* src_Dclose = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	for (i = 0; i<src->height; i++)
	{
		for (j = 0; j<src->width; j++)
		{
			if (cvGet2D(src_mgratio, i, j).val[0] <= alpha)
			{
				cvSet2D(src_Dclose, i, j, cvGet2D(src_close[s], i, j));
			}
			else if (cvGet2D(src_mgratio, i, j).val[0] >= beta)
			{
				cvSet2D(src_Dclose, i, j, cvGet2D(src_close[l], i, j));
			}
			else
			{
				cvSet2D(src_Dclose, i, j, cvGet2D(src_close[m], i, j));
			}
		}
	}

	// Step 5: divide I by L to get quotient image R
	IplImage* src_DMQIL = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, src->nChannels);
	for (i = 0; i < src->height; i++)
	{
		for (j = 0; j < src->width; j++)
		{
			double q = cvGet2D(src_denoise, i, j).val[0] / (cvGet2D(src_Dclose, i, j).val[0] + 1);
			cvSet2D(src_DMQIL, i, j, cvScalar(q));
		}
	}
	mat2gray(src_DMQIL);

	IplImage* src_DMQIL_8U = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, src->nChannels);
	cvConvertScaleAbs(src_DMQIL, src_DMQIL_8U, 230, 0);


	cvNamedWindow("src_DMQIL", CV_WINDOW_AUTOSIZE);
	cvShowImage("src_DMQIL", src_DMQIL_8U);

	return src_DMQIL_8U;
}