双三次插值及Matlab实现

双三次插值及Matlab实现

一、简单实例

  采用简单实例进行对双三次插值的介绍，由于双三次插值对于目标图像的某一像素进行估计时，所采用的像素信息为其周围16个像素点信息，因此不同于最近邻插值和双线性插值，此时假设有$5\times 5$大小的灰度图$src$，如下所示：
$$\left[\begin{array}{ccccc}0.1& 0.1& 0.1& 0.1& 0.1\\ 0.2& 0.2& 0.2& 0.2& 0.2\\ 0.3& 0.3& 0.3& 0.3& 0.3\\ 0.4& 0.4& 0.4& 0.4& 0.4\\ 0.5& 0.5& 0.5& 0.5& 0.5\end{array}\right]$$目标图像$dst$期望大小为$7\times 7$，则若要构建原始图像与目标图像之间的对应关系，根据之前的论述（最近邻插值及Matlab实现，后续有关扩展原始图像的内容也参看此篇），其映射公式为：
$$\{\begin{array}{l}srcX=(dstX+0.5)\times (srcWidth/dstWidth)-0.5\\ srcY=(dstY+0.5)\times (srcHeight/dstHeight)-0.5\end{array}$$其中$srcX$、$srcY$分别为目标图像$dstX$、$dstY$处所对应的原图像坐标，那么，假设现需要求解目标图像坐标为$(4,4)$位置处，即为目标图像中心处的像素估计值，则有：
$$\{\begin{array}{l}srcX=(4+0.5)\times (5/7)-0.5=2.714\\ srcY=(4+0.5)\times (5/7)-0.5=2.714\end{array}$$双三次插值所使用的像素信息为所计算的对应原始图像位置处周围的16个像素点，由于此时计算所得的坐标为$(2.714,2.714)$，因此其所对应的16个像素如下图所示：

红色点处即为所求点。双三次插值使用BiCubic基函数，其定义如下：
$$W(x)=\{\begin{array}{ll}(a+2)|x{|}^3-(a+3)|x{|}^2+1& \text{for |x| }⩽1\\ a|x{|}^3-5a|x{|}^2+8a|x|-4a& \text{for 1<|x|<2}\\ 0& \text{otherwise}\end{array}$$其中参数$a$为常系数，常取$a=-0.5$；变量$x$为红色点与16个对应像素之间的横向或纵向距离。双三次插值需要求解$x$、$y$两个方向上的权值，以左上方第一个像素点为例：①横向距离为$x=1.724$，则$W(x)=(-0.5)|1.724{|}^3-5(-0.5)|1.724{|}^2+8(-0.5)|1.724|-4(-0.5)=-0.0276$；②纵向距离为$y=1.724$，则$W(y)=-0.0276$；③第一个像素点对估计点造成的影响，即其权重则为$W(1)=W(x)\times W(y)=0.00076176$。同理，求得其余15个点对应的权重，再与各个点对应的数值相乘求和，所得结果即为估计点的数值大小。由于数据量大，此处便不再列出最终结果，且同理，此处需要用到边界的对称扩展。
  通过上述简单例子可以看出，双三次插值对于目标图像的每一个像素点的估计采用了其对应的16个像素点的加权组合，计算复杂度上升，但插值结果考虑了周围像素的变化率，使得放大效果的锐度得到了提高，如下图所示为利用双三次插值实现图像放大三倍的结果以及比较。

二、Matlab实现

Matlab实现如下：

function R = bicubic(src, scale)
%% 双三次插值
src = double(src) / 255;
% 判断是灰度图还是RGB图像
if ismatrix(src)
    R = zeros(floor(size(src) * scale));
else
    R = zeros([floor(size(src, 1, 2) * scale), 3]);
end
[dstM, dstN, ~] = size(R);
% 扩展原图像
misrc = zeros([size(src, 1, 2) + 2 * floor(scale), size(R, 3)]);
for i = 1 : size(R, 3)
    tmp = padarray(src(:, :, i), [floor(scale), floor(scale)], 'symmetric');
    misrc(:, :, i) = tmp;
end
 
%逐像素点赋值
for dstX = 1 : dstM
    for dstY = 1 : dstN
        srcX = floor((dstX + 0.5) / scale - 0.5);
        srcY = floor((dstY + 0.5) / scale - 0.5);
        u = ((dstX + 0.5) / scale - 0.5) - srcX;
        v = ((dstY + 0.5) / scale - 0.5) - srcY;
        X1 = zeros(4, 4);  
        X2 = zeros(4, 4);  
        W1 = ones(4, 4); 
        W2 = ones(4, 4);   
        % Bicubic基函数
        for i = 1 : 4
            for j = 1 : 4
                X1(i, j) = abs(u - i + 2);
                X2(i, j) = abs(v - j + 2);
                if X1(i, j) <= 1
                    W1(i, j) = 1.5 * (X1(i, j)) ^ 3 - 2.5 * (X1(i, j)) ^ 2 + 1;
                else
                    if X1(i, j) < 2
                        W1(i, j) = (-0.5) * (X1(i, j)) ^ 3 + 2.5 * (X1(i, j)) ^ 2 - 4 * X1(i, j) + 2;
                    else
                        W1(i, j) = 0;
                    end
                end
                if X2(i, j) <= 1
                    W2(i, j) = 1.5 * (X2(i, j)) ^ 3 - 2.5 * (X2(i, j)) ^ 2 + 1;
                else
                    if  X2(i, j) < 2
                        W2(i, j) = (-0.5) * (X2(i, j)) ^ 3 + 2.5 * (X2(i, j)) ^ 2 - 4 * X2(i, j) + 2;
                    else
                        W2(i, j) = 0;
                    end
                end
            end
        end
        W = W1 .* W2;
        Z = ones(4, 4);  %16个源像素点矩阵
        O = ones(4, 4);  %16个加权后的源像素点矩阵
        for dstC = 1 : size(R, 3)
            for i = 1 : 4
                for j = 1 : 4
                    Z(i, j) = misrc(srcX - 2 + i + round(scale), srcY - 2 + j + round(scale), dstC);
                    O(i, j) = W(i, j) .* Z(i,j);
                end
            end
            O1 = sum(sum(O));
            R(dstX, dstY, dstC) = O1;
        end
    end
end
end