#include "stdafx.h"
#include "histProject.h"
#include "omp.h" // 多核编程

// 功能：直方图投影，提高图像对比度
// 输入：    
//	        src： 输入图像
// 输出：  
//          dst:  输出图像
// 返回值：运算成功返回true，否则返回false
bool HP(cv::Mat& dst, const cv::Mat &src)
{
	try
	{
		// 数据有效性验证
		if (src.empty() == true || src.type() != 0)
		{
			return false;
		}

		// 计算直方图
		int channels[1];
		channels[0] = 0;

		int histSize[1];
		histSize[0] = 256;

		float hranges[2];
		hranges[0] = 0.0f;
		hranges[1] = 255.0f;

		const float * ranges[1];
		ranges[0] = hranges;

		cv::MatND hist;
		cv::calcHist(&src, 1, channels, cv::Mat(), hist, 1, histSize, ranges);

		// 直方图统计
		cv::Mat ZeroValuesStat(1, 256, CV_8UC1);
		ZeroValuesStat = 0;

		double srcMaxVal = 0;
		uchar zeroCount = 0;
		for (int i = 0; i < 256; i++)
		{
			ZeroValuesStat.at<uchar>(0, i) = zeroCount;

			if (hist.at<float>(i) == 0)
			{
				zeroCount++;
			}
			else
			{
				srcMaxVal = i;
			}
		}

		if (zeroCount == 0)
		{
			src.copyTo(dst);
			return true;
		}

		dst.release();
		dst = cv::Mat(src.rows, src.cols, CV_8UC1);

		// 剔除像素个数为0的像素值：向低阶压缩
		#pragma omp parallel for
		for (int r = 0; r < src.rows; r++)
		{
			for (int c = 0; c < src.cols; c++)
			{
				dst.at<uchar>(r, c) = src.at<uchar>(r, c) - ZeroValuesStat.at<uchar>(0, src.at<uchar>(r, c));
			}
		}

		// 灰度阶拉伸
		double maxVal = 255 - ZeroValuesStat.at<uchar>(0, 255);

		double ratio = 1.0;
		if (maxVal > 0)
		{
			ratio = 255.0 / maxVal;
		}

		dst.convertTo(dst, CV_8U, ratio, 0);

		return true;
	}
	catch(cv::Exception &e)
	{
		e.msg;
		return false;
	}
	catch(...)
	{
		return false;
	}
}

// 自适应阈值直方图投影，提高图像对比度
bool SAHP(cv::Mat &dst,const cv::Mat& src)
{
	// 数据有效性验证
	if (src.empty() == true || src.type() != CV_8UC1)
	{
		return false;
	}

	// 计算直方图
	int channels[1];
	channels[0] = 0;

	int histSize[1];
	histSize[0] = 256;

	float hranges[2];
	hranges[0] = 0.0f;
	hranges[1] = 255.0f;

	const float * ranges[1];
	ranges[0] = hranges;

	// 统计直方图
	cv::MatND hist;
	cv::calcHist(&src, 1, channels, cv::Mat(), hist, 1, histSize, ranges);

	// 直方图分析
	cv::Mat ZeroValuesStat(1, 256, CV_8UC1);
	ZeroValuesStat = 0;

	uchar zeroCount = 0;
	int threshold = (src.rows * src.cols) / (256 * 33); // 均值的3%倍
	for (int i = 0; i < 256; i++)
	{  
		ZeroValuesStat.at<uchar>(0, i) = zeroCount;

		if (hist.at<float>(i) < threshold)
		{
			zeroCount++;
		}
	}

 	if (zeroCount == 0)
	{
		src.copyTo(dst);
		return true;
	}

	dst.release();
	dst = cv::Mat(src.rows, src.cols, CV_8UC1);

	// 剔除像素个数为0的像素值：向低阶压缩
	for (int r = 0; r < src.rows; r++)
	{
		for (int c = 0; c < src.cols; c++)
		{
			dst.at<uchar>(r, c) = src.at<uchar>(r, c) - ZeroValuesStat.at<uchar>(0, src.at<uchar>(r, c));
		}
	}

	// 灰度阶拉伸
	double maxVal = 255 - ZeroValuesStat.at<uchar>(0, 255);

	double ratio = 1.0;
	if (maxVal > 0)
	{
		ratio = 255.0 / maxVal;
	}

	dst.convertTo(dst, CV_8UC1, ratio);

	return true;
}


// 高斯无约束直方图投影算法： 质量一般、速度快
bool GNRHP(cv::Mat& dst, const cv::Mat &src, bool bDeNoise)
{
	try
	{
		// 数据有效性验证
		if (src.empty() == true || src.type() != CV_8UC1)
		{
			return false;
		}

		// 备份
		cv::Mat srcCopy;
		src.copyTo(srcCopy);

		// 第一步：散粒噪声滤波
		if (bDeNoise == true)
		{
			cv::medianBlur(srcCopy, srcCopy, 3);
		}


		// 第二步：高斯无约束细节增强

		// 计算参数 T
		double minVal = 0, maxVal = 0; 
		cv::minMaxLoc(srcCopy, &minVal, &maxVal);

		double T = 2 * (maxVal - minVal);

		if (T == 0)
		{
			src.copyTo(dst);
			return false;
		}

		// 底图
		cv::Mat BaseImg;
		srcCopy.convertTo(BaseImg, CV_32FC1);

		// 高斯模糊 获取基图 3 * 1.667 = 5
		cv::GaussianBlur(BaseImg, BaseImg, cv::Size(5,5), 1.667);

		// 细节图
		cv::Mat DetailImg;
		cv::Mat srcF;
		srcCopy.convertTo(srcF, CV_32FC1);
		DetailImg = srcF - BaseImg;


		// 计算比例系数 自适应
		cv::minMaxLoc(BaseImg, &minVal, &maxVal);
		double BaseImgSpan = maxVal - minVal;
		cv::minMaxLoc(DetailImg, &minVal, &maxVal);
		double DetailImgSpan = maxVal - minVal;
		double ratio = BaseImgSpan / DetailImgSpan * 2.0;

		// 图像合成 可改进
		cv::Mat dstimg;
		cv::addWeighted(BaseImg, 1.0, DetailImg, ratio, 0, dstimg);

		cv::minMaxLoc(dstimg, &minVal, &maxVal);
		dstimg = dstimg - minVal;
		maxVal -= minVal;
		dstimg.convertTo(dstimg, CV_8U, 255.0 / maxVal, 0);

		// 细节增强图像HP(直方图投影）
		return HP(dst, dstimg);
	}
	catch(cv::Exception &e)
	{
		e.msg;
		return false;
	}
	catch(...)
	{
		return false;
	}
}


// 高斯约束直方图投影算法： 质量好、速度慢
bool GRHP(cv::Mat& dst, const cv::Mat &src, bool bDeNoise)
{
	try
	{
		// 数据有效性验证
		if (src.empty() == true || src.type() != 0)
		{
			return false;
		}

		// 备份
		cv::Mat srcCopy;
		src.copyTo(srcCopy);

		// 第一步：散粒噪声滤波
		if (bDeNoise == true)
		{
			cv::medianBlur(srcCopy, srcCopy, 3); // 用中值滤波代替
		}


		// 高斯滤波核尺寸
		int kernelSize = 7;
		double delta = 2.0;

		// 边缘复制
		cv::copyMakeBorder(srcCopy, srcCopy, kernelSize, kernelSize, kernelSize, kernelSize, BORDER_REPLICATE);

		// 第二步：高斯约束细节增强

		// 计算参数 T
		double minVal = 0, maxVal = 0; 
		cv::minMaxLoc(srcCopy, &minVal, &maxVal);

		double T = 2 * (maxVal - minVal);

		if (T == 0)
		{
			src.copyTo(dst);
			return false;
		}

		// 底图
		cv::Mat BaseImg;
		srcCopy.convertTo(BaseImg, CV_32FC1);

#pragma omp parallel for
		for (int r = kernelSize; r < srcCopy.rows - kernelSize; r++)
		{
			for (int c = kernelSize; c < srcCopy.cols - kernelSize; c++)
			{
				double val = 0;
				double sumS = 0;

				for (int i = -kernelSize; i < kernelSize + 1; i++)
				{
					for (int j = -kernelSize; j < kernelSize + 1; j++)
					{
						// 高斯滤波器系数：GC
						double GC = exp(-(i * i + j * j)/ (2 * delta * delta));

						// 灰度差值约束系数：DR
						double diff = srcCopy.at<uchar>(r + i, c + j) - srcCopy.at<uchar>(r, c);
						double DR = T / (pow(diff, 2) + T);

						// 合成系数
						double CC = GC * DR;

						// 累和
						val += CC * srcCopy.at<uchar>(r + i, c + j);
						sumS += CC;
					}
				}

				val = val / sumS;
				BaseImg.at<float>(r, c) = float(val);
			}
		}

		// 细节图
		cv::Mat DetailImg;
		cv::Mat srcF;
		srcCopy.convertTo(srcF, CV_32FC1);
		DetailImg = srcF - BaseImg;

		// BaseImg  DetailImg 裁边
		BaseImg   = BaseImg(cv::Range(kernelSize, BaseImg.rows - kernelSize),     cv::Range(kernelSize, BaseImg.cols - kernelSize));
		DetailImg = DetailImg(cv::Range(kernelSize, DetailImg.rows - kernelSize), cv::Range(kernelSize, DetailImg.cols - kernelSize));


		// 计算比例系数 自适应
		cv::minMaxLoc(BaseImg, &minVal, &maxVal);
		double BaseImgSpan = maxVal - minVal;
		cv::minMaxLoc(DetailImg, &minVal, &maxVal);
		double DetailImgSpan = maxVal - minVal;
		double ratio = BaseImgSpan / DetailImgSpan * 2.0;

		// 图像合成 可改进
		cv::Mat dstimg;
		cv::addWeighted(BaseImg, 1.0, DetailImg, ratio, 0, dstimg);

		cv::minMaxLoc(dstimg, &minVal, &maxVal);
		dstimg = dstimg - minVal;
		maxVal -= minVal;
		dstimg.convertTo(dstimg, CV_8U, 255.0 / maxVal, 0);

		// 细节增强图像HP(直方图投影）
		return HP(dst, dstimg);
	}
	catch(cv::Exception &e)
	{
		e.msg;
		return false;
	}
	catch(...)
	{
		return false;
	}
}