imresize-试图理解双三次插值。

Question

我正在尝试理解这个函数：

function [weights, indices] = contributions(in_length, out_length, ...
                                            scale, kernel, ...
                                            kernel_width, antialiasing)


if (scale < 1) && (antialiasing)
    % Use a modified kernel to simultaneously interpolate and
    % antialias.
    h = @(x) scale * kernel(scale * x);
    kernel_width = kernel_width / scale;
else
    % No antialiasing; use unmodified kernel.
    h = kernel;
end

我不是很明白这行代码的意思

 h = @(x) scale * kernel(scale * x);

我的比例尺是0.5
核函数是三次样条函数。

除此之外，它还有什么含义呢？ 我认为它就像创建一个稍后将被调用的函数一样。

- Gilad

Answer 1

这是对你关于MATLAB中imresize和OpenCV中cv::resize在使用双三次插值时差异的之前 问题的跟进。

我自己也很感兴趣，想要找出差异的原因。这是我的研究结果（如果我有任何错误，请指正）。

---

将图像大小调整视为从大小为M乘以N的输入图像到大小为scaledM乘以scaledN的输出图像的平面变换。

问题在于点不一定适配离散网格，因此为了获得输出图像中像素的强度，我们需要插值某些相邻样本的值（通常按照反向顺序执行，即对于每个输出像素，我们找到其对应的非整数点在输入空间中进行插值）。

这就是插值算法的差异之处，通过选择邻域的大小和给予该邻域中每个点的权重系数来实现。关系可以是一阶或高阶（涉及的变量是从反映射的非整数样本到原始图像网格上的离散点的距离）。通常将更接近的点赋予更高的权重。

查看MATLAB中的imresize，以下是线性和立方卷积核的权重函数：

function f = triangle(x)
    % or simply: 1-abs(x) for x in [-1,1]
    f = (1+x) .* ((-1 <= x) & (x < 0)) + ...
        (1-x) .* ((0 <= x) & (x <= 1));
end

function f = cubic(x)
    absx = abs(x);
    absx2 = absx.^2;
    absx3 = absx.^3;
    f = (1.5*absx3 - 2.5*absx2 + 1) .* (absx <= 1) + ...
        (-0.5*absx3 + 2.5*absx2 - 4*absx + 2) .* ((1 < absx) & (absx <= 2));
end

这些函数基本上返回样本与插值点之间的差距，从而确定样本所占的比重。

这些函数长这个样子:

>> subplot(121), ezplot(@triangle,[-2 2])  % triangle
>> subplot(122), ezplot(@cubic,[-3 3])     % Mexican hat

请注意，线性核（在[-1,0]和[0,1]间隔上的分段线性函数，在其他地方为零）适用于2个相邻点，而立方核（在区间[-2,-1]，[-1,1]和[1,2]上的分段三次函数，在其他地方为零）适用于4个相邻点。

以下是一维情况的示例，显示如何使用立方核从离散点f(x_k)中插值得到值x：

核函数 h(x) 以点 x 为中心，即待插值点的位置。插值值 f(x) 是离散相邻点（左边2个和右边2个）的加权和，乘上这些离散点处插值函数的值。

假设 x 与最近点之间的距离为 d (0 <= d < 1)，则在位置 x 处的插值值为：

f(x) = f(x1)*h(-d-1) + f(x2)*h(-d) + f(x3)*h(-d+1) + f(x4)*h(-d+2)

下面显示了点的顺序（注意 x(k+1)-x(k) = 1）：

x1      x2   x    x3       x4
o--------o---+----o--------o
         \___/
       distance d

现在，由于点是离散的且以均匀间隔采样，并且内核宽度通常很小，因此插值可以简洁地表述为卷积操作：

该概念在二维情况下，首先沿着一个维度进行插值，然后使用上一步骤的结果沿另一个维度进行插值即可实现。以下是双线性插值的示例，在二维中考虑4个相邻点：

在二维中双三次插值使用16个相邻点：

首先，我们使用16个网格样本(粉色)沿着行方向进行插值(红点)。然后，使用上一步中的插值点沿另一个维度进行插值(红线)。在每个步骤中，进行常规的一维插值。对我来说，这些方程太长、太复杂了，无法手工计算！

---

现在如果我们回到MATLAB中的“cubic”函数，它实际上与参考文献中所示的卷积核定义（方程4）匹配。下面是从维基百科中提取的同样内容：

你可以看到在上述定义中，MATLAB选择了a=-0.5的值。

现在MATLAB和OpenCV实现之间的差异是OpenCV选择了a=-0.75的值。

static inline void interpolateCubic( float x, float* coeffs )
{
    const float A = -0.75f;

    coeffs[0] = ((A*(x + 1) - 5*A)*(x + 1) + 8*A)*(x + 1) - 4*A;
    coeffs[1] = ((A + 2)*x - (A + 3))*x*x + 1;
    coeffs[2] = ((A + 2)*(1 - x) - (A + 3))*(1 - x)*(1 - x) + 1;
    coeffs[3] = 1.f - coeffs[0] - coeffs[1] - coeffs[2];
}

这可能不是立即显而易见的，但代码确实计算了三次卷积函数的项（在论文中的25号公式右侧列出）。

我们可以借助符号数学工具箱进行验证：

A = -0.5;
syms x
c0 = ((A*(x + 1) - 5*A)*(x + 1) + 8*A)*(x + 1) - 4*A;
c1 = ((A + 2)*x - (A + 3))*x*x + 1;
c2 = ((A + 2)*(1 - x) - (A + 3))*(1 - x)*(1 - x) + 1;
c3 = 1 - c0 - c1 - c2;

这些表达式可以重写为：

>> expand([c0;c1;c2;c3])
ans =
       - x^3/2 + x^2 - x/2
 (3*x^3)/2 - (5*x^2)/2 + 1
 - (3*x^3)/2 + 2*x^2 + x/2
             x^3/2 - x^2/2

从上面的方程式中可以匹配出相应的项。

显然，MATLAB和OpenCV之间的区别在于使用不同的自由项a值。根据论文作者的说法，0.5是首选值，因为与a的其他任何选择相比，它意味着更好的近似误差性能。

- Amro

---

非常好的回答。感谢您验证了符号数学工具箱中的方程式匹配。值得花费这些努力！ - chappjc

@chappjc：谢谢。很遗憾，参数a的值在两个实现中都是硬编码的。如果我们想让它们匹配，在MATLAB中，您需要修改内置的imresize函数（我从不喜欢这样做），而在OpenCV中，您需要编译整个源代码，只为了翻转一个细微的值！有趣的是，如果OP在OpenCV中将a更改为-0.5并验证我们在两个实现之间获得相同的结果，那将是很有意思的。上一次我尝试过，记得编译OpenCV从头开始要花费10分钟。 - Amro

@Gilad 我刚刚验证了将OpenCV更改为a = -0.5f与MATLAB给出的结果完全相同。我已经更新了我的答案以回答你的另一个问题，但我也会在这里添加一些内容。顺便说一下，使用普通硬盘，我的i7笔记本电脑编译大约需要6分钟。 :) 只需禁用CUDA！ - chappjc

我通常只是下载二进制文件。 :D - Gilad

顺便问一下，当你们提到抗锯齿时，是指哪个维基页面上的数值？http://en.wikipedia.org/wiki/Anti-aliasing 我猜是空间... - Gilad

Answer 2

imresize在缩小图像时实现抗锯齿效果，只需扩大立方卷积核，而不是离散的预处理步骤。

对于kernel_width为4个像素（缩放后为8个像素）的情况，其中contributions函数针对每个像素使用10个相邻点，kernel与h（缩放后的卷积核）如下图所示（未归一化，请忽略x轴）：

这比首先在单独的预处理步骤中执行低通滤波或高斯卷积更容易。

立方卷积核在imresize.m底部定义如下：

function f = cubic(x)
% See Keys, "Cubic Convolution Interpolation for Digital Image
% Processing," IEEE Transactions on Acoustics, Speech, and Signal
% Processing, Vol. ASSP-29, No. 6, December 1981, p. 1155.

absx = abs(x);
absx2 = absx.^2;
absx3 = absx.^3;

f = (1.5*absx3 - 2.5*absx2 + 1) .* (absx <= 1) + ...
                (-0.5*absx3 + 2.5*absx2 - 4*absx + 2) .* ...
                ((1 < absx) & (absx <= 2));

引用论文的PDF文件。

相关部分为方程式(15)：

这是以下方程式中的特定版本，用于a=-0.5的插值方程式：

a通常设置为-0.5或-0.75。请注意，a=-0.5对应于连续且具有连续一阶导数的三次Hermite样条。OpenCV似乎使用-0.75。

但是，如果您编辑[OPENCV_SRC]\modules\imgproc\src\imgwarp.cpp并更改代码：

static inline void interpolateCubic( float x, float* coeffs )
{
    const float A = -0.75f;
    ...

to:

static inline void interpolateCubic( float x, float* coeffs )
{
    const float A = -0.50f;
    ...

重新编译OpenCV（提示：禁用CUDA和gpu模块以缩短编译时间），然后您将获得与此相同的结果。请参阅我的另一个答案中与该问题相关的匹配输出。

- chappjc

---

@Gilad 我记得你在研究MATLAB和OpenCV的三次插值，看起来两者之间的差异是MATLAB的a=-0.5，而OpenCV的a=-0.75。 - chappjc

@chappjc：+1 不错的发现。我正好在写关于那个的答案 :) - Amro

@chappjc，我快要理解如何在C++中实现Matlab版本了。我认为我在信号处理方面缺乏一些基础知识，所以谢谢您。 - Gilad

@chappjc：我终于完成了我的回答。花费的时间比我预期的要长！ - Amro

@Gilad，你需要将OpenCV代码翻转使用-0.5。Amro的答案验证了该方程式与该常数除外的匹配。 - chappjc