www.gusucode.com > 《图像配准技术及其Matlab编程实现》--源码程序 > 《图像配准技术及其Matlab编程实现》/match/CubicConvolutionInterpolation.m

    function [newImage]=CubicConvolutionInterpolation(oldImage,w,z)
%参数：oldImage：待插值图像
%      w：反射变换参数x0（与oldImage同维数组）
%      z：反射变换参数y0（与oldImage同维数组）
%      newImage：插值后图像
%功能：利用立方卷积插值法根据几何变换参数进行灰度级插值

[Width,Height]=size(oldImage);
newImage=uint8(zeros(Width,Height));
for i=1:Width
    for j=1:Height
        source_x=w(i,j);
        source_y=z(i,j);
        if(source_x>=Width-3 || source_y>=Height-3 || double(uint16(source_x))<=2 || double(uint16(source_y))<=2)
            newImage(i,j)=0;
        else
            if(source_x/double(uint16(source_x))==1.0) & (source_y/double(uint16(source_y))==1.0)
                newImage(i,j)=oldImage(int16(source_x),int16(source_y));
            else
                %立方体卷积插值法
                a=fix(source_x);
                b=fix(source_y);
                u=(source_x -a);
                v=(source_y -b);
                a11=S(u+1);
                a12=S(u);
                a13=S(u-1);
                a14=S(u-2);
                A=[a11 a12 a13 a14];
                
                b11=double(oldImage(a-1,b-1));
                b12=double(oldImage(a-1,b));
                b13=double(oldImage(a-1,b+1));
                b14=double(oldImage(a-1,b+2));
                b21=double(oldImage(a,b-1));
                b22=double(oldImage(a,b));
                b23=double(oldImage(a,b+1));
                b24=double(oldImage(a,b+2));
                b31=double(oldImage(a+1,b-1));
                b32=double(oldImage(a+1,b));
                b33=double(oldImage(a+1,b+1));
                b34=double(oldImage(a+1,b+2));
                b41=double(oldImage(a+2,b-1));
                b42=double(oldImage(a+2,b));
                b43=double(oldImage(a+2,b+1));
                b44=double(oldImage(a+2,b+2));
                B=[b11 b12 b13 b14;b21 b22 b23 b24;b31 b32 b33 b34;b41 b42 b43 b44];
                c11=S(v+1);
                c21=S(v);
                c31=S(v-1);
                c41=S(v-2);
                C=[c11;c21;c31;c41];
                newImage(i,j)=A*B*C;
            end
        end
    end
end

figure;imshow(newImage);
%插值该函数
function s=S(w)
    w=abs(w);
    if w<1
        s=1-2*w^2+w^3;
    elseif w>=2
        s=0;
    else
        s=4-8*w+5*w^2-w^3;
    end
                
%RunInterpolation('lena.bmp',3);
%RunInterpolation('ct.bmp',3);