码迷，mamicode.com
首页 > 其他好文 > 详细

OpenCV2马拉松第9圈——再谈对比度(对比度拉伸，直方图均衡化)

时间：2014-05-11 13:40:50      阅读：450      评论：0      收藏：0      [点我收藏+]

标签：opencv   对比度拉伸   直方图均衡化   自使用直方图均衡化   lookup table   

收入囊中
  • lookup table
  • 对比度拉伸
  • 直方图均衡化
葵花宝典
lookup table是什么东西呢？
举个例子，假设你想把图像颠倒一下，f[i] = 255-f[i],你会怎么做？
for( int i = 0; i < I.rows; ++i)
  for( int j = 0; j < I.cols; ++j )
    I.at<uchar>(i,j) = 255 － I.at<uchar>(i,j);
大部分人应该都会这么做.或者：
for( i = 0; i < nRows; ++i){
        p = I.ptr<uchar>(i);
        for ( j = 0; j < nCols; ++j){
            p[j] = 255 － p[j];
        }
}
或者使用迭代器
MatIterator_<uchar> it, end;
            for( it = I.begin<uchar>(), end = I.end<uchar>(); it != end; ++it)
                *it = 255 － *it;


OpenCV提供了一个更快的方法，如下代码
LUT函数接收src,table和output
table是一个1*256的mat,将对应关系已经map好了
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
using namespace cv;

Mat applyLookUp(const cv::Mat& image,const cv::Mat& lookup) { 
	Mat result;
	cv::LUT(image,lookup,result);
	return result;
}

int main( int, char** argv )
{
  	Mat image,gray;
  	image = imread( argv[1], 1 );
  	if( !image.data )
  		return -1;
  	cvtColor(image, gray, CV_BGR2GRAY);
    
	Mat lut(1,256,CV_8U);
	for (int i=0; i<256; i++) {
		lut.at<uchar>(i)= 255-i;
	}

	Mat out = applyLookUp(gray,lut);
	namedWindow("sample");
	imshow("sample",out);

  	waitKey(0);
  	return 0;
}

第一种：On-The-Fly RA93.7878 milliseconds
第二种：Efficient Way79.4717 milliseconds
第三种：Iterator83.7201 milliseconds
第四种：LUT function32.5759 milliseconds


对比度拉伸又是什么？先来直观地看张图片。
bubuko.com,布布扣

右边是原始图，可以发现，低灰度没有像素，但是左边就比较好，低灰度也有，这就是对比度的拉伸。

公式非常简单：f[i] = 255.0*(i-imin)/(imax-imin)+0.5);
当灰度 < imin , f[i] = 0;
当灰度 > imax, f[i] = 255;
imin < f[i] < imax,就线性映射.

有人就会问了，imin和imax要怎么确定呢？imin取10还是20还是30呢？
我们可以确定一个阀值minvalue，当灰度的个数>这个阀值minvalue时，就确定下来了。看下面这个确定的过程：
histSize[0]就是256
Mat hist= getHistogram(image);
int imin= 0;
for( ; imin < histSize[0]; imin++ )
    if (hist.at<float>(imin) > minValue)
        break;

int imax= histSize[0]-1;
for( ; imax >= 0; imax-- )
    if (hist.at<float>(imax) > minValue)
        break;

一旦imin,imax确定下来，我们就可以填充lookup table了
Mat lookup(1, 256, CV_8U);
for (int i=0; i<256; i++) {
    if (i < imin) lookup.at<uchar>(i)= 0;
    else if (i > imax) lookup.at<uchar>(i)= 255;
    else lookup.at<uchar>(i)= static_cast<uchar>(255.0*(i-imin)/(imax-imin)+0.5);
}


直方图均衡化
这也是一个老生长谈的问题了
我再赘述一下
直方图均衡化的思想就是这样的。假设我有灰度级255的图像，但是都是属于［100，110］的灰度，图像对比度就很低，我应该尽可能拉到整个［0，255］
下面是直方图均衡化的代码，有个累积函数的概念，其实很简单。
我先计算出每个灰度级g(0),g(1)......g(255)点的个数,sum为图像width*height
那么累计函数c(0) = g(0)/sum
c(1) = (g(0)+g(1))/sum
......
c(255) = 1
下面是JAVA代码，改C++非常容易的
public int[][] Histogram_Equalization(int[][] oldmat)  
    {  
        int[][] new_mat = new int[height][width];  
        int[] tmp = new int[256];  
        for(int i = 0;i < width;i++){  
            for(int j = 0;j < height;j++){  
                //System.out.println(oldmat[j][i]);  
                int index = oldmat[j][i];  
                tmp[index]++;  
            }  
        }  
          
        float[] C = new float[256];  
        int total = width*height;  
        //计算累积函数  
        for(int i = 0;i < 256 ; i++){  
            if(i == 0)  
                C[i] = 1.0f * tmp[i] / total;  
            else  
                C[i] = C[i-1] + 1.0f * tmp[i] / total;  
        }  
          
        for(int i = 0;i < width;i++){  
            for(int j = 0;j < height;j++){  
                new_mat[j][i] = (int)(C[oldmat[j][i]] * 255);  
                new_mat[j][i] = new_mat[j][i] + (new_mat[j][i] << 8) + (new_mat[j][i] << 16);  
                //System.out.println(new_mat[j][i]);  
            }  
        }  
        return new_mat;  
    }

bubuko.com,布布扣



这是效果图，可以看到原来的图像被拉伸了

自适应直方图均衡化
AHE算法通过计算图像的局部直方图，然后重新分布亮度来来改变图像对比度。因此，该算法更适合于改进图像的局部对比度以及获得更多的图像细节。
想像以下一幅图像，左上角是黑乎乎的一团，但是其他区域很正常，如果只用HE,那么黑乎乎的那团是没法有多大改进的。
于是，你可以把那黑乎乎的一团当作一张图片，对那一部分进行HE,其实这就是AHE了，就是把图片分片处理，8*8是常用的选择。
然后，你就可以写一个循环来操作，算法和HE是一模一样的，当然可以工作，只是速度比较慢。
正如我下面代码所写的，利用双线性插值。
我以前写CLAHE时候看的博客找不到了T_T   http://m.blog.csdn.net/blog/gududeyhc/8997009这里有但是远远没我以前看的那篇讲的清楚，如果你去看Pizer的论文估计要花很多的时间。下面是我用Java写的CLAHE.
CLAHE比AHE多了裁剪补偿的操作
/* 
     * CLAHE 
     * 自适应直方图均衡化 
     */  
    public int[][] AHE(int[][] oldmat,int pblock)  
    {  
        int block = pblock;  
        //将图像均匀分成等矩形大小，8行8列64个块是常用的选择  
        int width_block = width/block;  
        int height_block = height/block;  
          
        //存储各个直方图  
        int[][] tmp = new int[block*block][256];  
        //存储累积函数  
        float[][] C = new float[block*block][256];  
          
        //计算累积函数  
        for(int i = 0 ; i < block ; i ++)  
        {  
            for(int j = 0 ; j < block ; j++)  
            {  
                int start_x = i * width_block;  
                int end_x = start_x + width_block;  
                int start_y = j * height_block;  
                int end_y = start_y + height_block;  
                int num = i+block*j;  
                int total = width_block * height_block;  
                for(int ii = start_x ; ii < end_x ; ii++)  
                {  
                    for(int jj = start_y ; jj < end_y ; jj++)  
                    {  
                        int index = oldmat[jj][ii];  
                        tmp[num][index]++;  
                    }  
                }  
                  
                  
                //裁剪操作  
                int average = width_block * height_block / 255;  
                int LIMIT = 4 * average;  
                int steal = 0;  
                for(int k = 0 ; k < 256 ; k++)  
                {  
                    if(tmp[num][k] > LIMIT){  
                        steal += tmp[num][k] - LIMIT;  
                        tmp[num][k] = LIMIT;  
                    }  
                      
                }  
                  
                int bonus = steal/256;  
                //hand out the steals averagely  
                for(int k = 0 ; k < 256 ; k++)  
                {  
                    tmp[num][k] += bonus;  
                }  
                  
                //计算累积分布直方图  
                for(int k = 0 ; k < 256 ; k++)  
                {  
                    if( k == 0)  
                        C[num][k] = 1.0f * tmp[num][k] / total;  
                    else  
                        C[num][k] = C[num][k-1] + 1.0f * tmp[num][k] / total;  
                }  
                  
            }  
        }  
          
        int[][] new_mat = new int[height][width];  
        //计算变换后的像素值  
        //根据像素点的位置，选择不同的计算方法  
        for(int  i = 0 ; i < width; i++)  
        {  
            for(int j = 0 ; j < height; j++)  
            {  
                //four coners  
                if(i <= width_block/2 && j <= height_block/2)  
                {  
                    int num = 0;  
                    new_mat[j][i] = (int)(C[num][oldmat[j][i]] * 255);  
                }else if(i <= width_block/2 && j >= ((block-1)*height_block + height_block/2)){  
                    int num = block*(block-1);  
                    new_mat[j][i] = (int)(C[num][oldmat[j][i]] * 255);  
                }else if(i >= ((block-1)*width_block+width_block/2) && j <= height_block/2){  
                    int num = block-1;  
                    new_mat[j][i] = (int)(C[num][oldmat[j][i]] * 255);  
                }else if(i >= ((block-1)*width_block+width_block/2) && j >= ((block-1)*height_block + height_block/2)){  
                    int num = block*block-1;  
                    new_mat[j][i] = (int)(C[num][oldmat[j][i]] * 255);  
                }  
                  
                //four edges except coners  
                else if( i <= width_block/2 )  
                {  
                    //线性插值  
                    int num_i = 0;  
                    int num_j = (j - height_block/2)/height_block;  
                    int num1 = num_j*block + num_i;  
                    int num2 = num1 + block;  
                    float p =  (j - (num_j*height_block+height_block/2))/(1.0f*height_block);  
                    float q = 1-p;  
                    new_mat[j][i] = (int)((q*C[num1][oldmat[j][i]]+ p*C[num2][oldmat[j][i]])* 255);  
                }else if( i >= ((block-1)*width_block+width_block/2)){  
                    //线性插值  
                    int num_i = block-1;  
                    int num_j = (j - height_block/2)/height_block;  
                    int num1 = num_j*block + num_i;  
                    int num2 = num1 + block;  
                    float p =  (j - (num_j*height_block+height_block/2))/(1.0f*height_block);  
                    float q = 1-p;  
                    new_mat[j][i] = (int)((q*C[num1][oldmat[j][i]]+ p*C[num2][oldmat[j][i]])* 255);  
                }else if( j <= height_block/2 ){  
                    //线性插值  
                    int num_i = (i - width_block/2)/width_block;  
                    int num_j = 0;  
                    int num1 = num_j*block + num_i;  
                    int num2 = num1 + 1;  
                    float p =  (i - (num_i*width_block+width_block/2))/(1.0f*width_block);  
                    float q = 1-p;  
                    new_mat[j][i] = (int)((q*C[num1][oldmat[j][i]]+ p*C[num2][oldmat[j][i]])* 255);  
                }else if( j >= ((block-1)*height_block + height_block/2) ){  
                    //线性插值  
                    int num_i = (i - width_block/2)/width_block;  
                    int num_j = block-1;  
                    int num1 = num_j*block + num_i;  
                    int num2 = num1 + 1;  
                    float p =  (i - (num_i*width_block+width_block/2))/(1.0f*width_block);  
                    float q = 1-p;  
                    new_mat[j][i] = (int)((q*C[num1][oldmat[j][i]]+ p*C[num2][oldmat[j][i]])* 255);  
                }  
                  
                //inner area  
                else{  
                    int num_i = (i - width_block/2)/width_block;  
                    int num_j = (j - height_block/2)/height_block;  
                    int num1 = num_j*block + num_i;  
                    int num2 = num1 + 1;  
                    int num3 = num1 + block;  
                    int num4 = num2 + block;  
                    float u = (i - (num_i*width_block+width_block/2))/(1.0f*width_block);  
                    float v = (j - (num_j*height_block+height_block/2))/(1.0f*height_block);  
                    new_mat[j][i] = (int)((u*v*C[num4][oldmat[j][i]] +   
                            (1-v)*(1-u)*C[num1][oldmat[j][i]] +  
                            u*(1-v)*C[num2][oldmat[j][i]] +  
                            v*(1-u)*C[num3][oldmat[j][i]]) * 255);  
  
                }  
                  
                new_mat[j][i] = new_mat[j][i] + (new_mat[j][i] << 8) + (new_mat[j][i] << 16);         
            }  
        }  
          
        return new_mat;  
          
    }

难道直方图均衡化的代码要让我们自己写？当然不是，下面就是API


初识API
C++: void equalizeHist(InputArray src, OutputArray dst)
 
  • src – Source 8-bit single channel image.
  • dst – Destination image of the same size and type as src .
TAT,这是目前为止最简单的 API了
内部好像不是用自适应直方图均衡化来做


荷枪实弹
先给出对比度拉伸的源代码
有一个我们上次用过的直方图类，加了一个拉伸的方法
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
using namespace cv;

Mat applyLookUp(const cv::Mat& image,const cv::Mat& lookup) { 
	Mat result;
	cv::LUT(image,lookup,result);
	return result;
}

class Histogram1D {
private:
	int histSize[1]; // number of bins
	float hranges[2]; // min and max pixel value
	const float* ranges[1];
	int channels[1];

public:
	Histogram1D() {
		histSize[0]= 256;
		hranges[0]= 0.0;
		hranges[1]= 255.0;
		ranges[0]= hranges;
		channels[0]= 0; // by default, we look at channel 0
	}
	
	Mat getHistogram(const cv::Mat &image) {
		Mat hist;
		calcHist(&image,1,channels,Mat(),hist,1,histSize,ranges);
		return hist;
	}
	
	Mat getHistogramImage(const cv::Mat &image){
		Mat hist= getHistogram(image);
		double maxVal=0;
		double minVal=0;
		minMaxLoc(hist, &minVal, &maxVal, 0, 0);
		Mat histImg(histSize[0], histSize[0],CV_8U,Scalar(255));
		int hpt = static_cast<int>(0.9*histSize[0]);
		for( int h = 0; h < histSize[0]; h++ ) {
			float binVal = hist.at<float>(h);
			int intensity = static_cast<int>(binVal*hpt/maxVal);
			line(histImg,Point(h,histSize[0]),
			Point(h,histSize[0]-intensity),
			Scalar::all(0));
		}
		return histImg;
	}
	
	Mat stretch(const cv::Mat &image, int minValue=0) {
		Mat hist= getHistogram(image);
		int imin= 0;
		for( ; imin < histSize[0]; imin++ )
			if (hist.at<float>(imin) > minValue)
				break;

		int imax= histSize[0]-1;
		for( ; imax >= 0; imax-- )
			if (hist.at<float>(imax) > minValue)
				break;

		Mat lookup(1, 256, CV_8U);
		for (int i=0; i<256; i++) {
			if (i < imin) lookup.at<uchar>(i)= 0;
			else if (i > imax) lookup.at<uchar>(i)= 255;
			else lookup.at<uchar>(i)= static_cast<uchar>(255.0*(i-imin)/(imax-imin)+0.5);
		}
		Mat result;
		result= applyLookUp(image,lookup);
		return result;
	}
	
};

int main( int, char** argv )
{
  	Mat image,gray;
  	image = imread( argv[1], 1 );
  	if( !image.data )
  		return -1;
  	cvtColor(image, gray, CV_BGR2GRAY);
    namedWindow("original");
	imshow("original",gray);
	
    	Histogram1D h;
	Mat streteched = h.stretch(gray,100);
	namedWindow("sample");
	imshow("sample",streteched);
	
	namedWindow("histogram1");
	imshow("histogram1",h.getHistogramImage(gray));
	namedWindow("histogram2");
	imshow("histogram2",h.getHistogramImage(streteched));

  	waitKey(0);
  	return 0;
}

再来看直方图均衡化代码
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
#include <stdio.h>

using namespace cv;
using namespace std;

int main( int, char** argv )
{
  Mat src, dst;

  const char* source_window = "Source image";
  const char* equalized_window = "Equalized Image";

  /// Load image
  src = imread( argv[1], 1 );

  if( !src.data )
    { cout<<"Usage: ./Histogram_Demo <path_to_image>"<<endl;
      return -1;
    }

  /// Convert to grayscale
  cvtColor( src, src, CV_BGR2GRAY );

  /// Apply Histogram Equalization
  equalizeHist( src, dst );

  /// Display results
  namedWindow( source_window, CV_WINDOW_AUTOSIZE );
  namedWindow( equalized_window, CV_WINDOW_AUTOSIZE );

  imshow( source_window, src );
  imshow( equalized_window, dst );

  /// Wait until user exits the program
  waitKey(0);

  return 0;


举一反三
我在上面给出了CLAHE的JAVA代码
这是一个很好的学习材料，双线性插值加速，附带剪裁补偿
如果你有时间，应该认真去看看，这是当初花了一天的时间写的TAT



计算机视觉讨论群：162501053
转载请注明：http://blog.csdn.net/abcd1992719g


OpenCV2马拉松第9圈——再谈对比度(对比度拉伸，直方图均衡化),布布扣,bubuko.com

OpenCV2马拉松第9圈——再谈对比度(对比度拉伸，直方图均衡化)

标签：opencv   对比度拉伸   直方图均衡化   自使用直方图均衡化   lookup table   

原文地址：http://blog.csdn.net/abcd1992719g/article/details/25483395
踩
(0)
赞
(0)
   
举报
评论 一句话评论（0）
登录后才能评论！
分享档案
更多>
2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)
周排行
mamicode.com排行更多图片
更多
友情链接
兰亭集智   国之画   百度统计   站长统计   阿里云   chrome插件   新版天听网
关于我们 - 联系我们 - 留言反馈
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码！