混合STR图谱识别全流程：GBDT分类、比例估计与去噪复现
2026/11/4 3:11:30
【修订版】本文已整理到我的个人博客，订正了原文的几处错误：采样间隔写反了，线性插值核应为 1−|t|，三次卷积核与 B 样条核的系数也有误。CSDN 这一版保留原样，之后的勘误与更新只在博客进行：Image Resize : 数字图像的整数倍扩大 | Finn's blog
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 评价用图片的制作方法
用相同的思想，我们这里用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月
貴家 仁志，“ディジタル画像の解像度変換 — 画像の拡大，” インターフェース，CQ出版，1998年4月
貴家 仁志，“ディジタル画像の表現と階調変換，色間引きの原理，” インターフェース，CQ出版，1998年2月
《Digital Image Processing》 Rafael C. Gonzalez / Richard E. Woods
原文发于博客： zhoufan900428-CSDN博客