双三次插值法(Bicubic interpolation)

• 双三次插值法(Bicubic interpolation)相对前两种算法计算过程更为复杂，考虑了待求像素坐标反变换后得到的浮点坐标周围的16个邻近像素。目的坐标M(x,y)和浮点坐标m(i+u,j+v)，需要选取插值基函数来拟合数据，最常用的插值基函数是对sin[(πx)/x]的逼近。
• 计算公式
  

function re_im = dthird(im, p, q)
%双三次插值法实现图像缩放，输入目标图像和行缩放、纵缩放倍数
%ziheng zhou 2016.3.30
im = im2double(im);
[m,n] = size(im);
im_R = im(:,:,1);
im_G = im(:,:,2);
im_B = im(:,:,3);
l = round(m*p);
h = round(n*q)/3;
re_R = zeros(l,h);
re_G = zeros(l,h);
re_B = zeros(l,h);

for dstx = 1:l
    for dsty = 1:h
        x = floor(dstx/p); u = dstx/p - x;
        y = floor(dsty/q); v = dsty/q - y;
        %逼近函数，A = s[]
        sa = 4 - 8*abs(u+1) + 5*abs(u+1)^2 - abs(u+1)^3;
        sb = 1 - 2*abs(u)^2 + abs(u)^3;
        sc = 1 - 2*abs(u-1)^2 + abs(u-1)^3;
        sd = 4 - 8*abs(u-2) + 5*abs(u-2)^2 - abs(u-2)^3;
        A = [sa, sb, sc, sd];
        %逼近函数，C = s[]
        sa = 4 - 8*abs(v+1) + 5*abs(v+1)^2 - abs(v+1)^3;
        sb = 1 - 2*abs(v)^2 + abs(v)^3;
        sc = 1 - 2*abs(v-1)^2 + abs(v-1)^3;
        sd = 4 - 8*abs(v-2) + 5*abs(v-2)^2 - abs(v-2)^3;
        C = [sa; sb; sc; sd];
        %考虑到边界像素，分类讨论得到插值计算区域矩阵 B[]
        xi = max(1,x-1); xx = max(1,x); xa = min(m,x+1); xaa = min(m,x+2);
        yi = max(1,y-1); yy = max(1,y); ya = min(n/3,y+1); yaa = min(n/3,y+2);
        if xi~=x-1||yi~=y-1||xaa~=x+2||yaa~=y+2
            BR = [im_R(xi,yi),im_R(xi,yy),im_R(xi,ya),im_R(xi,yaa);im_R(xx,yi),im_R(xx,yy),im_R(xx,ya),im_R(xx,yaa);
                im_R(xa,yi),im_R(xa,yy),im_R(xa,ya),im_R(xa,yaa);im_R(xaa,yi),im_R(xaa,yy),im_R(xaa,ya),im_R(xaa,yaa)];
            BG = [im_G(xi,yi),im_G(xi,yy),im_G(xi,ya),im_G(xi,yaa);im_G(xx,yi),im_G(xx,yy),im_G(xx,ya),im_G(xx,yaa);
                im_G(xa,yi),im_G(xa,yy),im_G(xa,ya),im_G(xa,yaa);im_G(xaa,yi),im_G(xaa,yy),im_G(xaa,ya),im_G(xaa,yaa)];
            BB = [im_B(xi,yi),im_B(xi,yy),im_B(xi,ya),im_B(xi,yaa);im_B(xx,yi),im_B(xx,yy),im_B(xx,ya),im_B(xx,yaa);
                im_B(xa,yi),im_B(xa,yy),im_B(xa,ya),im_B(xa,yaa);im_B(xaa,yi),im_B(xaa,yy),im_B(xaa,ya),im_B(xaa,yaa)];
        else
             BR = im_R(xi:xaa,yi:yaa);
             BG = im_G(xi:xaa,yi:yaa);
             BB = im_B(xi:xaa,yi:yaa);
        end
        %插值主运算 A*B*C
        re_R(dstx,dsty) = (A*BR)*C;
        re_G(dstx,dsty) = (A*BG)*C;
        re_B(dstx,dsty) = (A*BB*C);
    end
end
rm = cat(3,re_R,re_G,re_B);
re_im = im2uint8(rm);
figure,imshow(re_im);