數字影像處理，經典對比度增強演算法

EbowTang發表於2014-07-28
演算法

關於影像增強必須清楚的基本概念

1.影像增強的目的：

1）改善影像的視覺效果，
2）轉換為更適合於人或機器分析處理的形式
3）突出對人或機器分析有意義的資訊
4）抑制無用資訊，提高影像的使用價值
5）增強後的影像並不一定保真


2，影像增強的方法分類：

1）從處理物件分類：灰度影像，（偽）彩色影像
2）從處理策略分類：全域性處理，區域性處理（ROI ROI，Region of Interest Interest）
3）從處理方法分類：空間域（點域運算，即灰度變換；鄰域方法，即空域濾波），頻域方法
4）從處理目的分類：影像銳化，平滑去噪，灰度調整（對比度增強）


3，影像增強的方法之對比度增強

1）灰度變換法

線性變換（已實現）
對數變換（已實現）
指數變換（已實現）

2）直方圖調整法
直方圖均衡化（已實現）
直方圖匹配（未實現）



一，直方圖均衡化 

直方圖均衡化的英文名稱是Histogram Equalization. 

　　影像對比度增強的方法可以分成兩類:一類是直接對比度增強方法;另一類是間接對比度增強方法。直方圖拉伸和直方圖均衡化是兩種最常見的間接對比度增強方法。直方圖拉伸是通過對比度拉伸對直方圖進行調整，從而“擴大”前景和背景灰度的差別，以達到增強對比度的目的，這種方法可以利用線性或非線性的方法來實現;直方圖均衡化則通過使用累積函式對灰度值進行“調整”以實現對比度的增強。
　　直方圖均衡化處理的“中心思想”是把原始影像的灰度直方圖從比較集中的某個灰度區間變成在全部灰度範圍內的均勻分佈。直方圖均衡化就是對影像進行非線性拉伸，重新分配影像畫素值，使一定灰度範圍內的畫素數量大致相同。直方圖均衡化就是把給定影像的直方圖分佈改變成“均勻”分佈直方圖分佈。
　　缺點： 
　　1）變換後影像的灰度級減少，某些細節消失； 
　　2）某些影像，如直方圖有高峰，經處理後對比度不自然的過分增強。 
　　直方圖均衡化是影像處理領域中利用影像直方圖對對比度進行調整的方法。 
　　這種方法通常用來增加許多影像的區域性對比度，尤其是當影像的有用資料的對比度相當接近的時候。通過這種方法，亮度可以更好地在直方圖上分佈。這樣就可以用於增強區域性的對比度而不影響整體的對比度，直方圖均衡化通過有效地擴充套件常用的亮度來實現這種功能。
　　這種方法對於背景和前景都太亮或者太暗的影像非常有用，這種方法尤其是可以帶來X光影像中更好的骨骼結構顯示以及曝光過度或者曝光不足照片中更好的細節。這種方法的一個主要優勢是它是一個相當直觀的技術並且是可逆操作，如果已知均衡化函式，那麼就可以恢復原始的直方圖，並且計算量也不大。這種方法的一個缺點是它對處理的資料不加選擇，它可能會增加背景雜訊的對比度並且降低有用訊號的對比度。
　　

關於程式設計實現，同樣是不呼叫matlab庫函式，自己程式設計實現。這樣可以更深刻地理解直方圖均衡化技術，提高程式設計能力。

實現程式碼（matlab）：


clc;
close all;
clear all;
 
src_img = imread('flyman_gray.bmp');  

figure (1) 
subplot(321),imshow(src_img),title('原影像');%顯示原始影像  
subplot(322),imhist(src_img),title('原影像直方圖');%顯示原始影像直方圖  

matlab_eq=histeq(src_img);         %利用matlab的函式直方圖均衡化
subplot(323),imshow(matlab_eq),title('matlab直方圖均衡化原影像');%顯示原始影像  
subplot(324),imhist(matlab_eq),title('matlab均衡化後的直方圖');%顯示原始影像直方圖 

dst_img=myHE(src_img);             %利用自己寫的函式直方圖均衡化
subplot(325),imshow(dst_img),title('手寫均衡化效果');%顯示原始影像  
subplot(326),imhist(dst_img),title('手寫均衡化直方圖');%顯示原始影像直方圖

直方圖均衡化函式的實現：


function dst_img=myHE(src_img)  

[height,width] = size(src_img);
dst_img=uint8(zeros(height,width));
%進行畫素灰度統計;    
NumPixel = zeros(1,256);%統計各灰度數目，共256個灰度級    
for i = 1:height    
    for j = 1: width    
        NumPixel(src_img(i,j) + 1) = NumPixel(src_img(i,j) + 1) + 1;%對應灰度值畫素點數量增加一    
    end    
end    
%計算灰度分佈密度    
ProbPixel = zeros(1,256);    
for i = 1:256    
    ProbPixel(i) = NumPixel(i) / (height * width * 1.0);    
end    
%計算累計直方圖分佈    
CumuPixel = zeros(1,256);    
for i = 1:256    
    if i == 1    
        CumuPixel(i) = ProbPixel(i);    
    else    
        CumuPixel(i) = CumuPixel(i - 1) + ProbPixel(i);    
    end    
end    
  
% 指定範圍進行均衡化  
% pixel_max=max(max(I));  
% pixel_min=min(min(I));  
pixel_max=255;  
pixel_min=0;  
%對灰度值進行對映（均衡化）    
for i = 1:height    
    for j = 1: width    
        dst_img(i,j) = CumuPixel(src_img(i,j)+1)*(pixel_max-pixel_min)+pixel_min;    
    end    
end    
return;



為什們和matlab的直方圖不一樣呢？？？



二，指數變換

指數變換（Power-Law ）的公式：S=c*R^r,通過合理的選擇c和r可以壓縮灰度範圍，演算法以c=1.0/255.0, r=2實現。
要做該影像增強變換需要先做歸一化，再指數變換，最後反歸一化
增強效果展示：可以看見，改增強演算法並不能很好的將畫素儘可能的碾平。
指數增強參考程式為：
clc;
close all;
clear all; 
   
% -------------Gamma Transformations-----------------  
%f = imread('Fig0316(4)(bottom_left).tif');   
f = imread('seed.tif');   
Gamma = 0.4;  
g2 = myExpEnhance(f,Gamma);  

figure();  
subplot(221);  imshow(f);  xlabel('a).Original Image');  
subplot(222),imhist(f),title('原影像直方圖');%顯示原始影像直方圖  
subplot(223);  imshow(g2);  xlabel('b).Gamma Transformations \gamma = 0.4');  
subplot(224),imhist(g2),title('增強影像直方圖');%顯示原始影像直方圖
指數增強核心函式為：
function dst_img=myExpEnhance(src_img,Gamma)  
src_img = mat2gray(src_img,[0 255]);%將影像矩陣A中介於amin和amax的資料歸一化處理， 其餘小於amin的元素都變為0， 大於amax的元素都變為1。  
C = 1;  
g2 = C*(src_img.^Gamma); 
%反歸一化
max=255;
min=0;
dst_img=uint8(g2*(max-min)+min);




三，對數變換

       對數變換主要用於將影像的低灰度值部分擴充套件，將其高灰度值部分壓縮，以達到強調影像低灰度部分的目的。變換方法由下式給出。

這裡的對數變換，底數為（v+1），實際計算的時候，需要用換底公式。其輸入範圍為歸一化的【0-1】，其輸出也為【0-1】。對於不同的底數，其對應的變換曲線如下圖所示。

底數越大，對低灰度部分的強調就越強，對高灰度部分的壓縮也就越強。相反的，如果想強調高灰度部分，則用反對數函式就可以了。看下面的實驗就可以很直觀的理解，下圖是某影像的二維傅立葉變換影像，其為了使其灰度部分較為明顯，一般都會使用灰度變換處理一下。

效果圖：


參考程式碼：
clc;
close all;
clear all; 

%-------------Log Transformations-----------------
f = imread('seed.tif');

g_1 = myLogEnhance(f,10);
g_2 = myLogEnhance(f,100);
g_3 = myLogEnhance(f,200);

figure();
subplot(2,2,1);
imshow(f);xlabel('a).Original Image');

subplot(2,2,2);
imshow(g_1);xlabel('b).Log Transformations v=10');

subplot(2,2,3);
imshow(g_2);xlabel('c).Log Transformations v=100');

subplot(2,2,4);
imshow(g_3);
xlabel('d).Log Transformations v=200');

對數變換核心函式
function dst_img=myLogEnhance(src_img,v) 
c=1.0;
src_img = mat2gray(src_img,[0 255]);
g =c*log2(1 + v*src_img)/log2(v+1);
%反歸一化
max=255;
min=0;
dst_img=uint8(g*(max-min)+min);





四，灰度拉伸

灰度拉伸也用於強調影像的某個部分，與伽馬變換與對數變換不同的是，灰度拉昇可以改善影像的動態範圍。可以將原來低對比度的影像拉伸為高對比度影像。實現灰度拉昇的方法很多，其中最簡單的一種就是線性拉伸。而這裡介紹的方法稍微複雜一些。灰度拉伸所用數學式如下所示。

同樣的，其輸入r為【0-1】，其輸出s也為【0-1】。這個式子再熟悉不過了，跟巴特沃斯高通濾波器像極了，其輸入輸出關係也大致能猜到是個什麼形狀的。但是，這裡就出現一個問題了，輸入為0時候，式子無意義了。所以，在用Matlab計算的時候，將其變為如下形式。

這裡的eps，就是Matlab裡面，一個很小數。如此做的話，式子變得有意義了。但是，其輸入範圍為【0-1】的時候，其輸出範圍變為了。輸出範圍大致為【0-1】，為了精確起見，使用mat2gray函式將其歸一化到精確的[0-1]。呼叫格式如下。



五，線性拉伸

為了突出感興趣的目標或者灰度區間，相對抑制那些不感興趣的灰度區域，可採用分段線性法，常用的是三段線性變換




參考程式：

clc;
close all;
clear all; 

I=imread('seed.tif'); 
[m,n,k]=size(I);
figure (1)
imshow('seed.tif');title(' 原影像'); 
mid=mean(mean(I));
%橫軸
fa=20; fb=80;
%縱軸
ga=50; gb=230;

J=myLinearEnhance(I,fa,fb,ga,gb);
figure (2)
imshow(J);title(' 線性拉伸影像'); 

pixel_f=1:256;
pixel_g=zeros(1,256);

%三段斜率,小於1表示該段將會被收縮
k1=double(ga/fa); 
k2=(gb- ga)/(fb- fa);
k3=(256- gb)/(256- fb);
for i=1:256
    if i <= fa
        pixel_g(i)= k1*i;
    elseif fa < i && i <= fb
        pixel_g(i)= k2*( i- fa)+ ga;
    else
        pixel_g(i)= k3*( i - fb)+ gb;
    end
end
figure (3)
plot(pixel_f,pixel_g);


核心函式：

function dst_img=myLinearEnhance(src_img,fa,fb,ga,gb)  

[height,width] = size(src_img);
dst_img=uint8(zeros(height,width));

src_img=double(src_img);

%三段斜率
k1=ga/fa; 
k2=(gb- ga)/(fb- fa);
k3=(255- gb)/(255- fb);
for i=1:height
    for j=1:width
            if src_img(i,j) <= fa
                dst_img(i,j)= k1*src_img(i,j);
            elseif fa < src_img(i,j) && src_img(i,j) <= fb
                dst_img(i,j)= k2*( src_img(i,j)- fa)+ ga;
            else
                dst_img(i,j)= k3*( src_img(i,j)- fb)+ gb;
            end
    end
end
dst_img=uint8(dst_img);




附錄：

附錄網上的另一份講解：
直方圖均衡化演算法分為三個步驟，第一步是統計直方圖每個灰度級出現的次數，第二步是累計歸一化的直方圖，第三步是計算新的畫素值。
第一步：
for(i=0;i<height;i++)
for(j=0;j<width;j++)
n[s[i][j]]++;

for(i=0;i<L;i++)
p[i]=n[i]/(width*height);

這裡，n[i]表示的是灰度級為i的畫素的個數，L表示的是最大灰度級，width和height分別表示的是原始影像的寬度和高度，所以，p[i]表示的就是灰度級為i的畫素在整幅影像中出現的概率(其實就是p[]這個陣列儲存的就是這幅影像的歸一化之後的直方圖)。
第二步：
for(i=0;i<=L;i++)
for(j=0;j<=i;j++)
c[i]+=p[j];

c[]這個陣列儲存的就是累計的歸一化直方圖。
第三步：
max=min=s[0][0];
for(i=0;i<height;i++)
for(j=0;j<width;j++)
if(max<s[i][j]){
max=s[i][j];
}else if(min>s[i][j]){
min=s[i][j];
}

找出畫素的最大值和最小值。
for(i=0;i<height;i++)
for(j=0;j<width;j++)
t[i][j]=c[s[i][j]]*(max-min)+min;

t[][]就是最終直方圖均衡化之後的結果。


收錄優秀程式碼：

這份程式碼寫得不錯，學習了，原部落格地址見參考資源【3】！

#include <stdio.h>
#include <iostream>
#include "fftw3.h"
#include "string"
#include "vector"
#include <windows.h>
#include <opencv2/legacy/legacy.hpp>
#include <opencv2/nonfree/nonfree.hpp>//opencv_nonfree模組：包含一些擁有專利的演算法，如SIFT、SURF函式原始碼。 
#include "opencv2/core/core.hpp"
#include "opencv2/features2d/features2d.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <opencv2/nonfree/features2d.hpp>

using namespace cv;
using namespace std;

class hisEqt
{
public:
	hisEqt::hisEqt();
	hisEqt::~hisEqt();
public:
	int w;
	int h;
	int nlen;

	int *pHis;
	float *pdf;

	//=====求畫素分佈概率密度====  
	void  getPdf();

	//======統計畫素個數=======  
	void getHis(unsigned char*imgdata);

	//==========畫統計分佈直方圖===============  
	void drawHistogram(const float*pdf,Mat &hist1);  

	//===========直方圖均衡化==========  
	void hisBal();

	//====直方圖均衡化後的影像===  
	void imgBal(unsigned char* img);
};


hisEqt::hisEqt() :nlen(0){
	pHis = new int[256 * sizeof(int)];
	memset(pHis, 0, 256 * sizeof(int));
	pdf = new float[255 * sizeof(float)];
	memset(pdf, 0, 255 * sizeof(float));
}

hisEqt::~hisEqt(){
	delete[]pHis;
	delete[]pdf;
}


//======統計畫素個數=======  
void hisEqt::getHis(unsigned char*imgdata){
	for (int i = 0; i<nlen; i++)
	{
		pHis[imgdata[i]]++;
	}
}


//=====求畫素分佈概率密度====  
void hisEqt::getPdf(){
	for (int k = 0; k<256; k++)
	{
		pdf[k] = pHis[k] / float(nlen);
	}
}

//===========直方圖均衡化==========  
void hisEqt::hisBal(){
	for (int k = 1; k<256; k++)
	{
		pdf[k] += pdf[k - 1];
	}
	for (int k = 0; k<256; k++)
	{
		pHis[k] = 255 * pdf[k];
	}
}

//====直方圖均衡化  
void hisEqt::imgBal(unsigned char* img){
	for (int i = 0; i<nlen; i++)
	{
		img[i] = pHis[img[i]];
	}
}


void hisEqt::drawHistogram(const float *pdf, Mat& hist1){
	for (int k = 0; k<256; k++)
	{
		if (k % 2 == 0)
		{
			Point a(k, 255), b(k, 255 - pdf[k] * 2550);
			line(hist1,
				a,
				b,
				Scalar(0, 0, 255),
				1);
		}
		else
		{
			Point a(k, 255), b(k, 255 - pdf[k] * 2550);
			line(hist1,
				a,
				b,
				Scalar(0, 255, 0),
				1);
		}
	}
}


int main()
{
	Mat image = imread("Fig0651(a)(flower_no_compression).tif");
	if (!image.data)
		return -1;

	Mat hist2(256, 256, CV_8UC3, Scalar(0, 0, 0));
	Mat hist1(256, 256, CV_8UC3, Scalar(0, 0, 0));

	Mat imgOut = Mat(image.rows, image.cols, CV_8UC3, Scalar(0, 0, 0));
	vector<Mat> planes;
	int chn = image.channels();
	if (chn == 3)
	{
		split(image, planes);
	}
	while (chn)
	{
		chn--;
		unsigned char* imageData = new unsigned char[sizeof(unsigned char)*(image.cols*image.rows)];
		memcpy(imageData, planes[chn].data, planes[chn].cols*planes[chn].rows);
		hisEqt his;//自定義的類
		his.nlen = image.rows*image.cols;
		his.getHis(imageData);
		his.getPdf();

		//  //======畫原圖直方圖並儲存============  
		his.drawHistogram(his.pdf, hist1);
		string pic_name = "hisline";
		pic_name = pic_name + to_string(chn);
		pic_name=pic_name+	".jpg";
		imwrite(pic_name, hist1);

		his.hisBal();
		his.getPdf();
		//  //======畫均衡化後直方圖並儲存============  
		his.drawHistogram(his.pdf, hist2);
		string pic_name0 = "his_balanceline";
		pic_name0 = pic_name0 + to_string(chn);
		pic_name0 = pic_name0 + ".jpg";
		imwrite(pic_name0, hist2);

		//  //=====影像均衡化===  
		his.imgBal(imageData);
		memcpy(planes[chn].data, imageData, planes[chn].cols*planes[chn].rows);
		delete[] imageData;
		imageData = NULL;
	}
	merge(planes, imgOut);//單通道合併
	imwrite("result.jpg", imgOut);
	return 0;
}



參考資源

【1】http://blog.csdn.net/xiajun07061225/article/details/6910129

【2】數字影像處理，岡薩雷斯著

【3】http://blog.csdn.net/bettyshasha/article/details/46940805

【4】http://blog.csdn.net/terryzero/article/details/6043821

【5】http://www.myexception.cn/image/1450848.html

相關文章