图像处理学习笔记之MATLAB中imhist、imadjust、stretchlim函数实现

vector<int> imhist(Mat &srcImage, unsigned int n = 256)
{
	CV_Assert(srcImage.channels() == 1);
	vector<int> hist(n, 0);
	double a = n / 256.0;
	int index = 0;
	int rows = srcImage.rows;
	int cols = srcImage.cols;
	for (int i = 0; i < rows;i++)
	{
		uchar* pdata = srcImage.ptr<uchar>(i);
		for (int j = 0; j < cols;j++)
		{
			index = a*pdata[j];
			++hist[index];
		}
	}
	return hist;
}

void stretchlim(Mat &srcImage, Mat& lowHigh, vector<double> tol = vector<double> {0.01, 0.99})
{
    int nbins = 256;
    double tol_low, tol_high;
    double low, high;
    lowHigh.create (2, srcImage.channels (), CV_64F);
    switch (tol.size())
    {
    case 0:
        tol_low = 0.01;
        tol_high = 0.99;
        break;
    case 1:
        tol_low = tol[0];
        tol_high = 1 - tol[0];
        break;
    case 2:

        tol_low = tol[0];
        tol_high = tol[1];
        break;
    default:
        break;
    }
    if(tol_low > tol_high)
        return;
    int sum = srcImage.rows * srcImage.cols;
    //通道分离
    vector<Mat> channels;
    split(srcImage, channels);

    for(int i = 0; i < srcImage.channels(); i++)
    {
        //获取灰度直方图
        vector<int> hist = imhist(channels[i], nbins);
        for(int j = 0; j < hist.size (); ++j)
        {
            //计算累积分布
            double cumsum = std::accumulate(hist.begin(), hist.begin() + 1 + j, 0.0);
            //得到 >tol_low的累积分布概率的灰度等级
            if ((cumsum / sum) > tol_low)  // > tol_low
            {
                low = j / 255.0;
                break;
            }
        }
        //得到 >tol_high的分布概率的灰度等级
        for (int k = 0; k < hist.size(); k++)
        {
            double cumsum = std::accumulate(hist.begin(), hist.begin() + 1 + k, 0.0);
            if ((cumsum / sum) >= tol_high) // > tol_high
            {
                high = k / 255.0;
                break;
            }
        }
        if (low == high)
        {
            lowHigh.ptr<double>(i)[0] = 0;
            lowHigh.ptr<double>(i)[1] = 1;
        }
        else
        {
            lowHigh.ptr<double>(i)[0] = low;
            lowHigh.ptr<double>(i)[1] = high;
        }
    }
}

void imadjust(Mat& src, Mat& dst,Mat& lowHighIn, Mat&lowHighOut, double gamma=1)
{
	CV_Assert(src.data != NULL);

	int chl = src.channels();
	int rowNum = src.rows;
	int colNum = src.cols;

	//通道分离
	vector<Mat> channels;
	split(src, channels);
	
	
	//设置默认值
	if (lowHighIn.data==NULL)
	{
		lowHighIn=Mat::zeros(chl,2,CV_64F);
		for (int i = 0; i < chl; i++)
		{
			lowHighIn.at<double>(i, 1) = 1;
		}
	}
	
	if (lowHighOut.data==NULL)
	{
		lowHighOut = Mat::zeros(chl, 2, CV_64F);
		for (int i = 0; i < chl; i++)
		{
			lowHighOut.at<double>(i, 1) = 1;
		}
	}
	for (int m = 0; m < chl;m++)
	{
		//gamma校正查表
		vector<double> lookuptable(256, 0);
		vector<uchar> img(256,0);
		for (int i = 0; i < 256; i++)
		{
			lookuptable[i] = i / 255.0;
			if (lookuptable[i]<=lowHighIn.at<double>(m,0))
			{
				lookuptable[i] = lowHighIn.at<double>(m, 0);
			}
			if (lookuptable[i] >= lowHighIn.at<double>(m, 1))
			{
				lookuptable[i] = lowHighIn.at<double>(m, 1);
			}
			lookuptable[i] = (lookuptable[i] - lowHighIn.at<double>(m, 0)) / (lowHighIn.at<double>(m, 1) - lowHighIn.at<double>(m, 0));
			lookuptable[i] = pow(lookuptable[i], gamma);
			lookuptable[i] = lookuptable[i] * (lowHighOut.at<double>(m, 1) - lowHighOut.at<double>(m,0))+lowHighOut.at<double>(m,0);
			img[i] = lookuptable[i] * 255;
		}
		for (int j = 0; j < rowNum;j++)
		{
			for (int k = 0; k < colNum; k++)
			{
				channels[m].at<uchar>(j, k) = img[channels[m].at<uchar>(j, k)];
			}
		}
	}
	merge(channels, dst);
}

int main()
{
    Mat srcImage = imread("高圆圆.jpg");
	Mat grayImage, dstImage, lowHigh;
	cvtColor(srcImage, grayImage, CV_BGR2GRAY);
	Mat lh1,lh2;
	stretchlim(grayImage, lh1);
	imadjust(grayImage,dstImage, lh1, lh2,1);
    imshow("处理后的图像", dstImage);
    waitKey(0);
    return 0;
}