Frankot and Chellappa算法

Frankot and Chellappa算法C++代码实现。重建面型的梯度场与测得梯度场的总体误差最小化，通过傅里叶变换将一系列不可积的梯度场映射到频域中一系列可积的基本函数的组合

//傅里叶变换
void Fouriertransform（cv::Mat& src cv::Mat& complexImg）
{
	int M = cv::getOptimalDFTSize（src.rows）;
	int N = cv::getOptimalDFTSize（src.cols）;

	cv::Mat padded;
	cv::copyMakeBorder（src padded 0 M - src.rows 0 N - src.cols cv::BORDER_CONSTANT cv::Scalar::all（0））;

	cv::Mat planes[] = { cv::Mat_（padded） cv::Mat::zeros（padded.size（） CV_32F） };

	cv::merge（planes 2 complexImg）;

	dft（complexImg complexImg）;

}

void dftshift（cv::Mat& ds）
{
	int cx = ds.cols / 2;//图像的中心点x 坐标
	int cy = ds.rows / 2;//图像的中心点y 坐标
	cv::Mat q0 = ds（cv::Rect（0 0 cx cy））;//左上
	cv::Mat q1 = ds（cv::Rect（cx 0 cx cy））;//右上
	cv::Mat q2 = ds（cv::Rect（0 cy cx cy））;//左下
	cv::Mat q3 = ds（cv::Rect（cx cy cx cy））;//右下
	cv::Mat tmp;
	q0.copyTo（tmp）;
	q3.copyTo（q0）;
	tmp.copyTo（q3）;
	q1.copyTo（tmp）;
	q2.copyTo（q1）;
	tmp.copyTo（q2）;
}

//傅里叶逆变换
void InvFouriertransform（cv::Mat& src cv::Mat& dst）
{
	cv::idft（src dst/* cv::DFT_SCALE | cv::DFT_REAL_OUTPUT*/）;
	dst.convertTo（dst CV_32F）;

}

//保存数据点
void SavePointData（cv::Mat& srcstring& filename）
{
	printf（“保存数据中!\n“）;
	int rows = src.rows;
	int cols = src.cols;

	string filepath = “D:\\specular_surface\\calib\\Test\\“;
	filepath = filepath + filename;
	double time1 = cv::getTickCount（）;
	std::ofstream outfile;
	outfile.open（filepath）;
	for （int i = 0; i < rows; i++）
	{
		float * src_ptr = src.ptr（i）;

		for （int j = 0; j < cols; j++）
		{
			outfile << i << ““ << j << ““ << src_ptr[j] << std::endl;
		}
	}
	outfile.close（）;

	printf（“保存文件时间:%f\n“ （cv::getTickCount（） - time1） / cv::getTickFrequency（））;
}


//求取整幅图像绝对值得最大和最小值
void MinMaxAbValue（cv::Mat& src double & maxv double& minv）
{
	if （src.type（） != CV_32FC1）
	{
		src.convertTo（src CV_32FC1）;
	}
	int rows = src.rows;
	int cols = src.cols;

	maxv = 0;
	minv = cv::NORM_INF;

	for （int i = 0; i < rows; i++）
	{
		float * ptr = src.ptr（i）;
		for （int j = 0; j < cols; j++）
		{
			if （abs（ptr[j]） > maxv）
			{
				maxv = abs（ptr[j]）;
			}
			if （abs（ptr[j]） < minv）
			{
				minv = abs（ptr[j]）;
			}
		}
	}

}


//乘以虚数j
void MultiplyJ（cv::Mat& src cv::Mat& dst）
{
	cv::Mat planes[2];

	cv::split（src planes）;

	cv::Mat real image;
	real = planes[0];
	image = planes[1];

	cv::Mat dd = cv::Mat::zeros（src.size（） CV_32FC1）;

	cv::subtract（dd real real）;

	cv::Mat temp = image;
	planes[1] = real;
	planes[0] = image;

	cv::merge（planes 2 dst）;
}


//计算梯度
void CalcGradient（const cv::Mat1f& matPhasetan cv::Mat& matGradient_X cv::Mat& matGradient_Y）
{
	matGradient_X = cv::Mat1f（matPhasetan.rows matPhasetan.cols）;
	matGradient_Y = cv::Mat1f（matPhasetan.rows matPhasetan.cols）;
	for （int i = 0; i < matPhasetan.rows; i++）
	{
		for （int j = 0; j < matPhasetan.cols; j++）