//percent_i means the weight that the ith derivative takes in the final result
double percent_0 = 0.5;
double percent_1 = 0.5;

//function declaration

double calc_first_derivative(CFile* file1, CFile* file2);
double calc_first_derivative_per_frame(CFile* file1,CFile* file2);
double first_derivative_average(vector<double> &value);

//function definition
double calc_first_derivative(CFile* file1, CFile* file2)
{

	// first_derivative_per_frame saves first derivative of each frame
	vector<double> first_derivative_frame;

	while ((file1 -> GetPosition() != file1 -> GetLength()) && (file2 -> GetPosition() != file2 -> GetLength()))
	{
		first_derivative_frame.push_back(calc_first_derivative_per_frame(file1,file2));
	}

	//calculate total first_derivative
	return first_derivative_average(first_derivative_frame);
}

double calc_first_derivative_per_frame(CFile* file1,CFile* file2)
{
	//set buffer and malloc memory
	unsigned char *y_data1 = new unsigned char[y_size],*y_data2 = new unsigned char[y_size];
	unsigned char *cb_data1 = new unsigned char[cb_size],*cb_data2 = new unsigned char[cb_size];
	unsigned char *cr_data1 = new unsigned char[cr_size],*cr_data2 = new unsigned char[cr_size];

	//read file
	file1 -> Read(y_data1,y_size);
	file1 -> Read(cb_data1,cb_size);
	file1 -> Read(cr_data1,cr_size);
	file2 -> Read(y_data2,y_size);
	file2 -> Read(cb_data2,cb_size);
	file2 -> Read(cr_data2,cr_size);

	//set 2-D buffer
	unsigned char (*y_data1_2d)[width] = new unsigned char[height][width],(*y_data2_2d)[width] = new unsigned char[height][width];
	//unsigned char (*cb_data1_2d)[width] = new unsigned char[height][width],(*cb_data2_2d)[width] = new unsigned char[height][width];
	//unsigned char (*cr_data1_2d)[width] = new unsigned char[height][width],(*cr_data2_2d)[width] = new unsigned char[height][width];

	//horizontal[][] saves the relationship of horizontal items
	//vertical[][] saves the relationship of vertical items
	int (*y_horizontal1)[width - 1] = new int[height][width - 1],(*y_vertical1)[width] = new int[height - 1][width];
	//int (*cb_horizontal1)[width - 1] = new int[height][width - 1],(*cb_vertical1)[width] = new int[height - 1][width];
	//int (*cr_horizontal1)[width - 1] = new int[height][width - 1],(*cr_vertical1)[width] = new int[height - 1][width];
	int (*y_horizontal2)[width - 1] = new int[height][width - 1],(*y_vertical2)[width] = new int[height - 1][width];
	//int (*cb_horizontal2)[width - 1] = new int[height][width - 1],(*cb_vertical2)[width] = new int[height - 1][width];
	//int (*cr_horizontal2)[width - 1] = new int[height][width - 1],(*cr_vertical2)[width] = new int[height - 1][width];

	//convert 1-D to 2-D 
	//(cb & cr must be treated differently from y!)
	for (int h = 0;h < height;++h)
	{
		for (int w = 0;w < width;++w)
		{
			y_data1_2d[h][w] = y_data1[h * width + w];
			//cb_data1_2d[h][w] = cb_data1[h * width + w];
			//cr_data1_2d[h][w] = cr_data1[h * width + w];
			y_data2_2d[h][w] = y_data2[h * width + w];
			//cb_data2_2d[h][w] = cb_data2[h * width + w];
			//cr_data2_2d[h][w] = cr_data2[h * width + w];
		}
	}

	//calculate horizontal and vertical relationship
	////for data1
	for (int h = 0;h < height;++h)
	{
		for (int w = 0;w < width - 1;++w)
		{
			y_horizontal1[h][w] = y_data1_2d[h][w + 1] - y_data1_2d[h][w];
			//cb_horizontal1[h][w] = cb_data1_2d[h][w + 1] - cb_data1_2d[h][w];
			//cr_horizontal1[h][w] = cr_data1_2d[h][w + 1] - cr_data1_2d[h][w];
		}
	}
	for (int h = 0;h < height - 1;++h)
	{
		for (int w = 0;w < width;++w)
		{
			y_vertical1[h][w] = y_data1_2d[h + 1][w] - y_data1_2d[h][w];
			//cb_vertical1[h][w] = cb_data1_2d[h + 1][w] - cb_data1_2d[h][w];
			//cr_vertical1[h][w] = cr_data1_2d[h + 1][w] - cr_data1_2d[h][w];
		}
	}
	////for data2
	for (int h = 0;h < height;++h)
	{
		for (int w = 0;w < width - 1;++w)
		{
			y_horizontal2[h][w] = y_data2_2d[h][w + 1] - y_data2_2d[h][w];
			//cb_horizontal2[h][w] = cb_data2_2d[h][w + 1] - cb_data2_2d[h][w];
			//cr_horizontal2[h][w] = cr_data2_2d[h][w + 1] - cr_data2_2d[h][w];
		}
	}
	for (int h = 0;h < height - 1;++h)
	{
		for (int w = 0;w < width;++w)
		{
			y_vertical2[h][w] = y_data2_2d[h + 1][w] - y_data2_2d[h][w];
			//cb_vertical2[h][w] = cb_data2_2d[h + 1][w] - cb_data2_2d[h][w];
			//cr_vertical2[h][w] = cr_data2_2d[h + 1][w] - cr_data2_2d[h][w];
		}
	}

	//sum
	unsigned __int64 sum_horizontal = 0;
	for (int h = 0;h < height;++h)
	{
		for (int w = 0;w < width - 1;++w)
		{
			sum_horizontal += (y_horizontal1[h][w] - y_horizontal2[h][w]) * (y_horizontal1[h][w] - y_horizontal2[h][w]);
		}
	}

	unsigned __int64 sum_vertical = 0;
	for (int h = 0; h < height - 1;++h)
	{
		for (int w = 0;w < width;++w)
		{
			sum_vertical += (y_vertical1[h][w] - y_vertical2[h][w]) * (y_vertical1[h][w] - y_vertical2[h][w]);
		}
	}

	unsigned __int64 sum = sum_horizontal + sum_vertical;

	//free memory
	//free 1_D arrays
	delete [] y_data1;
	delete [] cb_data1;
	delete [] cr_data1;
	delete [] y_data2;
	delete [] cb_data2;
	delete [] cr_data2;

	//free 2_D arrays
	//delete y_data1_2d
	delete [] y_data1_2d;
	
	//delete [] cb_data1_2d;
	//delete [] cr_data1_2d;

	//delete y_data2_2d
	delete [] y_data2_2d;

	//delete [] cb_data2_2d;
	//delete [] cr_data2_2d;

	//delete y_horizontal1
	delete [] y_horizontal1;

	//delete y_vertical1
	delete [] y_vertical1;

	//delete [] cb_horizontal1;
	//delete [] cb_vertical1;
	//delete [] cr_horizontal1;
	//delete [] cr_vertical1;

	//delete y_horizontal2
	delete [] y_horizontal2;

	//delete y_vertical2
	delete [] y_vertical2;

	//delete [] cb_horizontal2;
	//delete [] cb_vertical2;
	//delete [] cr_horizontal2;
	//delete [] cr_vertical2;

	//return first_derivative_per_frame
	//first_derivative is the result when considered just by first_derivative
	double first_derivative = 0;
	if (sum == 0)
		first_derivative = 100;
	else
	{
		//scheme 1
		//the worst case happens when one value is 255 and the other is -255,so the worst MSE is (255-(-255))^2 = 510*510
		double first_derivative_horizontal = 10 * log10((510 * 510) / double(sum_horizontal / (height * (width - 1))));
		double first_derivative_vertical = 10 * log10((510 * 510) / double(sum_vertical / ((height - 1) * width)));
		first_derivative = (first_derivative_horizontal + first_derivative_vertical) / 2;

		//scheme 2
		//return 10 * log10(510 * 510 / double(sum / (height * (width - 1) + (height - 1) * width)));
	}

	double psnr = calc_psnr(file1,file2);

	//final result
	if (psnr < 30 && first_derivative > 30)
	{
		percent_0 = 0.7;
		percent_1 = 0.3;
	}
	else if (psnr > 30 && first_derivative < 30)
	{
		percent_0 = 0.3;
		percent_1 = 0.7;
	}
	else
	{
		percent_0 = 0.5;
		percent_1 = 0.5;
	}
	
	double result = percent_0 * psnr + percent_1 * first_derivative;
	return result;
}

double first_derivative_average(vector<double> &value)
{
	double sum = 0;
	int i = 0;
	for (vector<double>::iterator iter = value.begin();iter != value.end();++iter)
	{
		++i;
		sum += (*iter);
	}
	return sum / i;
}