[数字图像处理]数字图像的整数倍扩大(数字图像插值)

【修订版】本文已整理到我的个人博客订正了原文的几处错误采样间隔写反了线性插值核应为 1−|t|三次卷积核与 B 样条核的系数也有误。CSDN 这一版保留原样之后的勘误与更新只在博客进行Image Resize : 数字图像的整数倍扩大 | Finns blog1.序言在不同的平台下或者从不同仪器获得图像或者数据其大小与数据类型都有着很大的不同。这里的大小就是指的是分辨率。数字图像其实是像素的集中表现形式像素越多越密集图像则可以表现得越精确的。我们将数字图像的像素数称为分辨率。本文主要介绍了数字图像的整数倍扩大也就是分辨率的变化。事实上扩大处理可以归结为新的像素值的如何决定的插值问题。首先作为目标先对理想的插值进行理解。然后对现有的常用插值方法进行一个讨论对其性能进行全面的分析。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所得到的就是我们所需要的滤波器的单位冲击响应。例如扩大倍数U2的时候所需要的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(x1,y1); 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(U1,:); g_bilin((M*U)U1,:) g_bilin((M*U)1,:); g_bilin(:,1) g_bilin(:,U1); g_bilin(:,(M*U)U1) 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(x3,y3); 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 (a2)*(abs(t).^(3))-(a3)*(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(U1,:); g_cubic((M*U)U1,:) g_cubic((M*U)1,:); g_cubic(:,1) g_cubic(:,U1); g_cubic(:,(M*U)U1) g_cubic(:,(M*U)1); for x 1:1:(U*M16) g_cubic(x,:) filter(H_cubic,1,g_cubic(x,:)); end for y 1:1:(U*N16) 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(x6,y6); 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(U1,:); g_B_sp((M*U)U1,:) g_B_sp((M*U)1,:); g_B_sp(:,1) g_B_sp(:,U1); g_B_sp(:,(M*U)U1) 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(x6,y6); 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原文发于博客 zhoufan900428-CSDN博客