MATLAB數字影象處理（二）影象增強

天靖居士發表於2017-03-06

1 影象增強

1.1 直方圖均衡化

對於灰度影象，可以使用直方圖均衡化的方法使得原影象的灰度直方圖修正為均勻的直方圖。

程式碼如下：

I2=histeq(I1);
figure,imshow(I2);
figure,imhist(I2);

原影象為lena的圖片，經過直方圖均衡化後的圖片為Figure 1。

Figure 1

原圖灰度直方圖請參見：MATLAB數字影象處理（一）基礎操作和傅立葉變換，經過修正後的灰度直方圖為Figure 2：

Figure 2

可以看出，與原圖對比，灰度分佈更加均勻。

1.2 灰度變換

在MATLAB中，可以使用imadjust函式對影象進行線性變換。常用的為線性變換。使用的函式為imadjust(I, [low_in, high_in], [low_out high_out],gamma)，其中，in為輸入的圖片變化的灰度範圍，out為輸出的圖片變化的灰度範圍，均為[0 1]之間。也就是將low_in和high_in之間的值對映到low_out和high_out之間。gamma為線性值，若為1則為線性變換。

程式碼和效果如下：

I3=imadjust(I1,[0.3 0.7],[0.1 0.9],1);
figure,imshow(I3);
figure,imhist(I3);

Figure 3

Figure 4

其中，突出了原圖0.3-0.7之間的灰度元素。

1.3 空間域平滑

影象平滑是為了抑制影象噪聲，改善影象質量而進行的圖片處理。在這裡，筆者將首先介紹MATLAB中自帶的均值濾波和中值濾波實現，之後設計自己的模版運算元，再次分別實現了均值濾波和中值濾波。

MATLAB中自帶的均值濾波可以採用fspecial(‘params’,n)函式實現，其中，params代表的是濾波的種類，常用的有average（均值濾波），還有可以在空間域銳化中用到的高斯、拉普拉斯、prewitt、sobel運算元等。n代表的是模版的在這裡僅僅介紹均值濾波。

MATLAB中，均值濾波程式碼如下：

avgModel=fspecial('average',3);
Iavg=filter2(avgModel,I2)/255;
figure,imshow(Iavg);

由於均值濾波對於高斯噪聲效果較好，故對高斯噪聲的圖（參見MATLAB數字影象處理（一）基礎操作和傅立葉變換）進行了影象平滑。效果如下（Figure 5）：

Figure 5

而筆者根據均值濾波原理，自己寫了一段程式碼，如下：

n=3;
model(1:n,1:n)=1;
Iavg=I2;
for i=2:length(Iavg)-1
    for j=2:length(Iavg(2,:))-1
        x=I2(i-(n-1)/2:i+(n-1)/2,j-(n-1)/2:j+(n-1)/2);
        xuint=uint8(model).*uint8(x);
        x_avg=mean(xuint(:));
        Iavg(i,j)=x_avg;
    end
end
figure,imshow(Iavg);

效果（Figure 6）：

Figure 6

可以發現，除了執行速度稍慢，效果與MATLAB自帶的函式相同。

MATLAB中自帶的中值濾波可以採用medfilt2(I,[n n])函式實現，其中，I代表的是需要處理的圖片，n代表的是模版的大小。

MATLAB中，中值濾波程式碼如下：

Imid=medfilt2(I3,[3,3]);
figure,imshow(Imid);

由於中值濾波對於椒鹽噪聲效果較好，故對椒鹽噪聲的圖（參見MATLAB數字影象處理（一）基礎操作和傅立葉變換）進行了影象平滑。效果如下（Figure 7）：

Figure 7

而筆者根據中值濾波原理，自己寫了一段程式碼，如下：

n=3;
Imid=I3;
for i=2:length(I1)-1
    for j=2:length(I1(2,:))-1
        x=I3(i-(n-1)/2:i+(n-1)/2,j-(n-1)/2:j+(n-1)/2);%3*3
        x_mid=median(x(:));
        Imid(i,j)=x_mid;
    end
end
figure,imshow(Imid);

效果（Figure 8）：

Figure 8

可以發現，除了執行速度稍慢，效果與MATLAB自帶的函式相同。

1.4 空間域銳化

影象銳化有多種方法，最常用的銳化運算元包括梯度運算元、Roberts運算元、Prewitt運算元、Sobel運算元以計算梯度的梯度銳化法，Laplacian運算元等。接下來，筆者將一一介紹相關的運算元。

筆者根據梯度運算元的計算原理，編寫了如下程式碼：

Igrad=I1;
for i=1:length(Igrad)-1
    for j=1:length(Igrad(2,:))-1
        x=I1(i,j+1)-I1(i,j);
        y=I1(i+1,j)-I1(i,j);
        grad=max(abs(x),abs(y));
        Igrad(i,j)=grad;%¸³Öµ
    end
end
figure,imshow(Igrad);

效果（Figure 9）：

Figure 9

MATLAB中自帶了Prewitt運算元和Sobel運算元的函式，程式碼及效果如下：

model=fspecial('prewitt');
Iprewitt=filter2(model,I1);
figure,imshow(Iprewitt);

Figure 10

model=fspecial('sobel');
Isobel=filter2(model,I1);
figure,imshow(Isobel);

Figure 11

事實上，可以發現這幾種運算元之間不同的地方只是在於模版，因此筆者編寫了一個可以自定義各種模版的程式碼，如下：

n=2
modelx=[-1 0;0 1];
modely=[0 -1;1 0];
Iend=I1;
Idouble=double(I1);
for i=1:length(Idouble)-5+n
    for j=1:length(Idouble(2,:))-5+n
        area=Idouble(i:i+n-1,j:j+n-1);
        x=area.*modelx;
        y=area.*modely;
        grad=max(abs(sum(x(:))),abs((sum(y(:)))));
        Iend(i,j)=grad;%¸³Öµ
    end
end
figure,imshow(Iend);

可以顯示出Roberts運算元的效果：

Figure 12

與Prewitt運算元、Sobel運算元經過對比後發現效果相同，說明有效。

對於Laplacian運算元，同樣可以使用MATLAB中自帶的函式，也可以進行自定義。自帶函式使用方法如下：

model=fspecial('laplacian');
Ilaplacian=filter2(model,I1);
figure,imshow(Ilaplacian,[]);

效果：

Figure 13

套用前面自定義模版後的程式碼如下：

n=3;
model=[0 -1 0;-1 5 -1;0 -1 0];
Iend=I1;
Idouble=double(I1);
for i=1:length(Idouble)-5+n
    for j=1:length(Idouble(2,:))-5+n
        area=Idouble(i:i+n-1,j:j+n-1);
        x=area.*model;
        grad=sum(x(:));
        Iend(i,j)=grad;
    end
end
figure,imshow(Iend);

效果：

Figure 14

與教科書上的效果相同。

1.5 頻率域影象增強

頻率域增強包括高通濾波器和低通濾波器，主要是使用傅立葉變換和逆變換對影象進行處理。在MATLAB中沒有直接的濾波器函式，但是可以結合原理與MATLAB中的數學公式，構造各種濾波器。在本文中，筆者分別構造了理想、Butterworth、指數的高通和低通濾波器，並一一進行了實驗。

首先需要提前設定截止頻率D0和階數n。之後就可以套用具體的公式進行計算了。部分程式碼如下：

n=6;%濾波器的階數
D0=80;%濾波器的截止頻率
[u,v]=meshgrid(-M/2:(M/2-1),-N/2:(N/2-1));%產生離散資料
D=sqrt(u.^2+v.^2);
%各種不同的濾波器
H1=double(D<=D0);%理想低通濾波器
H2=1./(1+(D/D0).^(2*n));%Butterworth低通濾波器
H3=exp(-(D/D0).^n);%指數低通濾波器
H4=double(D>=D0);%理想高通濾波器
H5=1./(1+(D0./D).^(2*n));%Butterworth高通濾波器
H6=exp(-(D0./D).^n);%指數高通濾波器

之後，可以使用各種不同的濾波器對圖片進行處理。效果如下：