[數字圖像處理]數字圖像的整數倍擴大(數字圖像插值)

1.序言

       在不同的平台下，或者從不同儀器獲得圖像（或者資料），其大小與資料類型，都有着很大的不同。這裡的 大小，就是指的是分辨率。 數字圖像，其實是像素的集中表現形式，像素越多越密集，圖像則可以表現得越精确的。我們将數字圖像的像素數，稱為分辨率。        本文主要介紹了數字圖像的整數倍擴大（也就是分辨率的變化）。事實上，擴大處理可以歸結為，新的像素值的如何決定的插值問題。首先，作為目标，先對理想的插值進行了解。然後，對現有的常用插值方法進行一個讨論，對其性能進行全面的分析。

2.理想的插值處理

       首先，作為目标，先對理想的插值進行說明。為了便于了解，我們假設有如下這樣一個模拟信号

。

左邊是我們假設的模拟型号在時間域内的表現，右邊是其頻譜。先這樣假設着，然後，我們将這個信号進行采樣。

       這裡，

是一個沖擊串函數，所得到的是一個離散的信号

。如下所示，求所得到的離散信号

的離散傅裡葉變換。

看上式，并求其離散信号的振幅譜，所得到的振幅譜變為了一個周期性的，表示為如下形式。

這裡，我們選擇的采樣時間間隔為

。一般的，還會将得到的振幅譜歸一化，将橫軸變為歸一化頻率。所得到的振幅譜的兩個波峰之間差一個2π。

       為了說明理想插值，我們假設如下兩個信号。首先，選擇采樣間隔時間為

，可以得到如下信号①。

然後，将其采樣時間間隔設定為0.5倍（也就是将采樣頻率設定為信号①的2倍），也就是

，可以得到如下的信号②。

觀察上述兩個信号，①與②。我們可以得到這樣兩個資訊。一，兩個波峰之間的歸一化頻率的差都是2π。二，信号②的振幅是信号①的兩倍。對于理想的插值處理，我們進行如下定義。

如果通過某個處理，可以從信号①變換到信号②的話，這樣的插值處理就稱為理想插值。

       下面，從頻域上來觀察一下理想插值處理。如下圖所示，

最左邊是信号①的頻譜，最右邊是信号②的頻譜（也就是目标信号）。首先，在信号①的時間域的相鄰的信号中間，依次插入1個0（擴大倍率是U的時候，插入0的個數是U-1個0。上述例子的U為2，故插入一個0。），所得到的信号的振幅譜就像中間的那樣（這個操作相當于升采樣處理[Upsampling]）。下一步，就是從中間的頻譜得到目标頻譜。可以看到，如果可以将中間信号的頻譜的π處的成分完全衰減，剩餘的部分增幅2倍的話，便可以得到目标信号。将頻域的某個部分衰減，同時增幅某個部分，這正是濾波器的作用。

       是以，為了實作理想的插值，就需要如下濾波器。

根據上圖（濾波器的振幅特性），我們可以知道，為實作理想插值，所需要的濾波器是理想濾波器。之前的數字信号進行中，我們已經知道，理想濾波器是無法實作（因為其機關沖擊響應是無限的）。是以，我們隻能接近理想的濾波器。換句話說，我們可以得到這樣一個結論！某種插值方法的好壞，可以由其相當的濾波器的振幅特性與理想濾波器的接近程度來判斷。越接近理想濾波器，則可以得到更好的插值效果。

3.常用插值方法的分析

       本小節主要對4種常用插值方法，零次保持法、線性插值法、 cubic convolution插值法、B-spline插值法的性能進行分析。

       3.1 零次保持法

       這個方法，也就是最鄰近插值法。舉個栗子來說，就像下圖。

左邊的是輸入信号

，右邊的是輸出信号

。這個方法，相當于如下所示的一維濾波器。

       3.2 線性插值法

       這個方法，就是将相鄰的信号之間以直線相連，将直線上的值作為内插值的方法。可以參考下圖。

左邊的是輸入信号

，右邊的是輸出信号

。這個方法，相當于如下所示的一維濾波器。

這個表示方法可能有點難以了解。當擴大倍數為U時候，t值取-1到1，步進值為1/U，所得到的

就是我們所需要的濾波器的機關沖擊響應。例如，擴大倍數U=2的時候，所需要的h(t) = [0.5 1 0.5]。

       3.3 Cubic Convolution插值法

       這個方法也就是三次卷積插值法，這個方法相當于如下所示的一維濾波器。

與之前的線性插值法一樣，當擴大倍數為U時候，t值取-1到1，步進值為1/U，所得到的

就是我們所需要的濾波器的機關沖擊響應。        這裡有一個參數a，調整參數a可以調整此插值法的性能，如下圖所示。

        3.4 B-Spline插值法

        這個方法相當于如下所示的一維濾波器。

這個方法沒有參數。

        3.5 上述四種方法的分析

        将上述四種方法相當于的濾波器和理想濾波器的振幅特性繪制出來，如下圖所示。

根據上圖，我們可以得到如下這樣幾個結論。        ①：插值效果最差的插值法是零次保持法（最鄰近插值法），這個濾波器跟理想濾波器差得最遠。        ②：比起零次保持法，線性插值法更加接近理想濾波器，是以線性插值法的效果比零次保持法更好。        ③：cubic convolution插值法則更加接近理想濾波器，是以，性能比線性插值法要好。但是，其機關沖擊響應有負值，會産生“振鈴”現象，是以需要小心的調整參數a。         ④： B-spline 插值法，高頻特性最接近理想濾波器，但是低頻特性是距理想濾波器最遠的。

4.常用插值方法的實作（Matlab）與實際插值效果分析

        4.1 評價用圖檔的制作方法

       首先，說明一下實驗用的圖檔。在第2節中，我們通過用不同的采樣頻率，制作出兩個信号①和②。我們将信号②作為目标信号，定義了理想的插值。換句話說，也就是信号②是信号①理想插值的結果。        用相同的思想，我們這裡用Matlab制作出兩枚圖檔。圖檔②是圖檔①理想插值的結果。這兩枚圖檔的制作方法如下圖所示。    

使用上述方法，我們可以得到一下兩枚實驗用圖檔。其分辨率為：256×256  512×512兩枚。

         4.2 實驗順序

         本次實驗的順序如下所示。

本次實驗主要實作的是圖像的分辨率的2倍擴大。由于需要評價插值法的效果，将2倍擴大後的圖檔與目标圖檔做差，得到差分圖檔。本次實驗，我們使用SAD(Sum of Absolute Difference)作為評價參數，也就是差分圖像的絕對值的和。SAD的值越大，插值效果越差。

        4.3 實驗結果與Matlab代碼

close all;
clear all;
clc;

f = imread('./cheer/maru_256(cheer).tif');
f_Goal = imread('./cheer/maru_512(cheer).tif');
f = mat2gray(f,[0 255]);
f_Goal = mat2gray(f_Goal,[0 255]);

U = 2;  %拡大率　

[M,N] = size(f);
g_hold = zeros((M*U),(N*U));
H_hold = ones(1,U);

for x = 0:1:M-1     
   for y = 0:1:N-1
       g_hold((U*x)+1,(U*y)+1) = f(x+1,y+1); 
   end
end

for x = 1:1:U*M 
    g_hold(x,:) = filter(H_hold,1,g_hold(x,:));
end
for y = 1:1:U*N 
    g_hold(:,y) = filter(H_hold,1,g_hold(:,y));
end

g_hold_diff = zeros((M*U),(N*U));
for x = 1:1:(M*U)     
   for y = 1:1:(N*U)
       g_hold_diff(x,y) = abs(f_Goal(x,y) - g_hold(x,y)); 
   end
end

SAD_hold = sum(sum(g_hold_diff))
%SSD_hold = sum(sum(g_hold_diff .^2))
%MSE_hold = sum(sum(g_hold_diff .^2))/((M*U)*(M*U));
%PSNR_hold = 10*log10((1*1)/MSE_hold)

figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image (256x256)');

figure();
subplot(1,2,1);
imshow(g_hold,[0 1]);
xlabel('a).Ruselt of hold (512x512)');

subplot(1,2,2);
imshow(g_hold_diff,[0 1]);
xlabel('b).Difference image (512x512)');

figure();
w = 0:0.01:pi;
Xejw_1 = freqz(H_hold,1,w);
plot(w,abs(Xejw_1));
axis([0,pi,0,2]);grid;
xlabel('\omega [rad]');
ylabel('|H(e^{j\omega_1})|');
f = imread('./cheer/maru_256(cheer).tif');
f_Goal = imread('./cheer/maru_512(cheer).tif');
f = mat2gray(f,[0 255]);
f_Goal = mat2gray(f_Goal,[0 255]);

U = 2;   %拡大率　

[M,N] = size(f);
g_bilin = zeros((M*U)+6,(N*U)+6);

H_bilin = [1 2 1]/2;

for x = 1:1:M    
   for y = 1:1:N
       g_bilin((U*x)+1,(U*y)+1) = f(x,y); 
   end
end

g_bilin(1,:) = g_bilin(U+1,:); 
g_bilin((M*U)+U+1,:) = g_bilin((M*U)+1,:); 
g_bilin(:,1) = g_bilin(:,U+1); 
g_bilin(:,(M*U)+U+1) = g_bilin(:,(M*U)+1); 

for x = 1:1:(U*M)+6 
    g_bilin(x,:) = filter(H_bilin,1,g_bilin(x,:));
end
for y = 1:1:(U*N)+6 
    g_bilin(:,y) = filter(H_bilin,1,g_bilin(:,y));
end

g_b = zeros((M*U),(N*U));
for x = 1:1:(M*U)
   for y = 1:1:(M*U)
       g_b(x,y) = g_bilin(x+3,y+3); 
   end
end

g_b_diff = zeros((M*U),(N*U));
for x = 1:1:(M*U)    
   for y = 1:1:(N*U)
       g_b_diff(x,y) = abs(f_Goal(x,y) - g_b(x,y)); 
   end
end

SAD_b = sum(sum(g_b_diff))
%SSD_b = sum(sum(g_b_diff .^2))
%MSE_b = sum(sum(g_b_diff .^2))/((M*U)*(M*U));
%PSNR_b = 10*log10((1*1)/MSE_b)

figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image (256x256)');

figure();
subplot(1,2,1);
imshow(g_b,[0 1]);
xlabel('a).Ruselt of bilinear interpolation (512x512)');

subplot(1,2,2);
imshow(g_b_diff,[0 1]);
xlabel('b).Difference image (512x512)');

figure();
w = 0:0.01:pi;
Xejw_2 = freqz(H_bilin,1,w);
plot(w,abs(Xejw_2));
axis([0,pi,0,2]);grid;
xlabel('\omega [rad]');
ylabel('|H(e^{j\omega_1}|');
f = imread('./cheer/maru_256(cheer).tif');
f_Goal = imread('./cheer/maru_512(cheer).tif');
f = mat2gray(f,[0 255]);
f_Goal = mat2gray(f_Goal,[0 255]);

U = 2;   %拡大率　

[M,N] = size(f);
g_cubic = zeros((M*U)+16,(N*U)+16);

a = -0.0001;
t = -1+(1/U):(1/U):1-(1/U);
H_cubic = (a+2)*(abs(t).^(3))-(a+3)*(abs(t).^(2)) + 1;  
t = -2+(1/U):(1/U):-1;
H_cubic = [0 (a)*(abs(t).^(3))-(a*5)*(abs(t).^(2))+(a*8)*abs(t)-4*a H_cubic];  
t = 1:(1/U):2-(1/U);
H_cubic = [H_cubic (a)*(abs(t).^(3))-(a*5)*(abs(t).^(2))+(a*8)*abs(t)-4*a 0];

for x = 1:1:M 
   for y = 1:1:N
       g_cubic((U*x)+1,(U*y)+1) = f(x,y); 
   end
end

g_cubic(1,:) = g_cubic(U+1,:); 
g_cubic((M*U)+U+1,:) = g_cubic((M*U)+1,:); 
g_cubic(:,1) = g_cubic(:,U+1); 
g_cubic(:,(M*U)+U+1) = g_cubic(:,(M*U)+1); 

for x = 1:1:(U*M+16) 
    g_cubic(x,:) = filter(H_cubic,1,g_cubic(x,:));
end
for y = 1:1:(U*N+16) 
    g_cubic(:,y) = filter(H_cubic,1,g_cubic(:,y));
end

g_c = zeros((M*U),(N*U));
for x = 1:1:(M*U)
   for y = 1:1:(M*U)
       g_c(x,y) = g_cubic(x+6,y+6); 
   end
end

g_c_diff = zeros((M*U),(N*U));
for x = 1:1:(M*U) 
   for y = 1:1:(N*U)
       g_c_diff(x,y) = abs(f_Goal(x,y) - g_c(x,y)); 
   end
end

SAD_c = sum(sum(g_c_diff))
%SSD_c = sum(sum(g_c_diff .^2))
%MSE_c = sum(sum(g_c_diff .^2))/((M*U)*(M*U));
%PSNR_c = 10*log10((1*1)/MSE_c)

figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image (256x256)');

figure();
subplot(1,2,1);
imshow(g_c,[0 1]);
xlabel('a).Ruselt of Cubic convolution (a=-0.5)(512x512)');

subplot(1,2,2);
imshow(g_c_diff,[0 1]);
xlabel('b).Difference image (512x512)');

figure();
w = 0:0.01:pi;
Xejw_3 = freqz(H_cubic,1,w);
plot(w,abs(Xejw_3));
axis([0,pi,0,5]);grid;
xlabel('\omega [rad] (a = -0.41)');
ylabel('|H(e^{j\omega_1}|');
f = imread('./cheer/maru_256(cheer).tif');
f_Goal = imread('./cheer/maru_512(cheer).tif');
f = mat2gray(f,[0 255]);
f_Goal = mat2gray(f_Goal,[0 255]);

U = 2;   %拡大率　

[M,N] = size(f);
g_B_sp = zeros((M*U)+16,(N*U)+16);

t = -1+(1/U):(1/U):1-(1/U);
H_B_sp = (1/2)*(abs(t).^(3))-(abs(t).^(2)) + (2/3);  
t = -2+(1/U):(1/U):-1;
H_B_sp = [0 (-1/6)*(abs(t).^(3))+(abs(t).^(2))-(2)*abs(t)+(4/3) H_B_sp];  
t = 1:(1/U):2-(1/U);
H_B_sp = [ H_B_sp (-1/6)*(abs(t).^(3))+(abs(t).^(2))-(2)*abs(t)+(4/3) 0];  

for x = 1:1:M  
   for y = 1:1:N
       g_B_sp((U*x)+1,(U*y)+1) = f(x,y); 
   end
end

g_B_sp(1,:) = g_B_sp(U+1,:); 
g_B_sp((M*U)+U+1,:) = g_B_sp((M*U)+1,:); 
g_B_sp(:,1) = g_B_sp(:,U+1); 
g_B_sp(:,(M*U)+U+1) = g_B_sp(:,(M*U)+1); 

for x = 1:1:(M*U)+16 
    g_B_sp(x,:) = filter(H_B_sp,1,g_B_sp(x,:));
end
for y = 1:1:(N*U)+16
    g_B_sp(:,y) = filter(H_B_sp,1,g_B_sp(:,y));
end

g_B = zeros((M*U),(N*U));
for x = 1:1:(M*U)
   for y = 1:1:(M*U)
       g_B(x,y) = g_B_sp(x+6,y+6); 
   end
end

g_B_diff = zeros((M*U),(N*U));
for x = 1:1:(M*U) 
   for y = 1:1:(N*U)
       g_B_diff(x,y) = abs(f_Goal(x,y) - g_B(x,y)); 
   end
end

SAD_B = sum(sum(g_B_diff))
%SSD_B = sum(sum(g_B_diff .^2))
%MSE_B = sum(sum(g_B_diff .^2))/((M*U)*(M*U));
%PSNR_B = 10*log10((1*1)/MSE_B)


figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image (512x512)');

figure();
subplot(1,2,1);
imshow(g_B,[0 1]);
xlabel('b).Ruselt of B-spline (512x512)');

subplot(1,2,2);
imshow(g_B_diff,[0 1]);
xlabel('b).Difference image (512x512)');

figure();
w = 0:0.01:pi;
Xejw_4 = freqz(H_B_sp,1,w);
plot(w,abs(Xejw_4));
axis([0,pi,0,5]);grid;
xlabel('\omega [rad] (a = -1)');
ylabel('|H(e^{j\omega_1}|');
           

參考文獻 貴家 仁志，“ディジタル畫像の解像度変換 — 任意サイズへの畫像の変換，” インターフェース，CQ出版，1998年6月

貴家 仁志，“ディジタル畫像の解像度変換 — 畫像の拡大，” インターフェース，CQ出版，1998年4月

貴家 仁志，“ディジタル畫像の表現と階調変換，色間引きの原理，” インターフェース，CQ出版，1998年2月

《Digital Image Processing》 Rafael C. Gonzalez / Richard E. Woods

原文發于 部落格： http://blog.csdn.net/thnh169/