�禺画像]
バイキュービックの計算には、16画素のデータそれぞれに対する係数を求める必要がある。

求める画素位置と、実際の画素の位置から、画素間の距離を求め、その距離に応じた係数をそれぞれの画素の輝度値に掛け、その総和を求める画素の輝度値とする。
ここで、距離と係数の関係は sinc 関数を３次近似したものが使われる。
上の画像のピンクの線が sinc 関数で、下の緑と青が近似した値である。見やすくするために sinc 関数は上に0.5だけ移動してある。
ちなみに、sinc 関数とは sin(πx)/(πx) のこと。

この部分を関数にするとこんな感じ。

function fsinc(d)
        return (d < 1) ? d^3 - 2d^2 + 1 : (d < 2) ? -d^3 + 5d^2 - 8d + 4 : 0
end

実は julia には sinc 関数があるので、別の名前にした。

x方向、y方向、それぞれ独立に係数を求め、掛け合わせることで16個の係数を求める。（試しに xy両方を含めた距離で計算してみたが、ノイズっぽい画像になってしまった）

x方向だけ考えると、求める画素の位置は、４つの画素の真ん中にくる。
x値を小さい順に x1, x2, x3, x4 と置くと、求める画素の位置(x)は x2≦x＜x3 となり、x2= floor(x) である。
よって、x と x2 の距離は x - floor(x) であり、x と x1 の距離は x と x2 の距離に１を足したものである。
また、x と x3 の距離は、１から(x と x2 の距離)を引いたものとなり、x と x4 の距離はそれに１を足したものである。
よって、x と x2 の距離（すなわち、x - floor(x) ）を xd とすれば、それぞれ xd+1, xd, 1-xd, 2-xd となる。
y に関しても同様。

これを関数で書くと以下のようになる。

function getCOEFFS (xin, yin)
        # coefficient (a= -1)
        function fsinc(d)
                return (d < 1) ? d^3 - 2d^2 + 1 : (d < 2) ? -d^3 + 5d^2 - 8d + 4 : 0
        end
        # calc coefficients
        xd= xin - floor(xin)
        xc= map(x -> fsinc(x), [xd+1; xd; 1-xd; 2-xd])
        #
        yd= yin - floor(yin)
        yc= map(x -> fsinc(x), [yd+1 yd 1-yd 2-yd])
        #
        return xc * yc
end

関数の中に関数が書けるのは便利。
map を使って、行列の距離値を sinc 関数を使った係数に変換している。

xc は、4行１列の行列。yc は、１行４列の行列。
これらを掛け合わせると、４行４列の行列が得られる。これを戻り値にしている。

完成した関数は以下のとおり。

# ibuff の (xin,yin) から画素値を取得する
function imgpic (ibuff, xin, yin)

        # get 4x4 pixels
        function getPIX4x4 (ibuff, xin, yin)
                xstart= ifloor(xin) - 1
                ystart= ifloor(yin) - 1
                try
                        return ibuff[xstart:(xstart+3), ystart:(ystart+3)]
                catch   # out of range
                        obuff= ibuff[1:4,1:4]
                        (xmax,ymax)= size(ibuff)
                        for x= 1:4, y= 1:4
                                xin= xstart + (x - 1)
                                yin= ystart + (y - 1)
                                xin= (xin < 1) ? 1 : (xin > xmax) ? xmax : xin
                                yin= (yin < 1) ? 1 : (yin > ymax) ? ymax : yin
                                obuff[x,y]= ibuff[xin,yin]
                        end
                        return obuff
                end
        end

        # get 4x4 coefficients
        function getCOEFFS (xin, yin)
                # coefficient (a= -1)
                function fsinc(d)
                        return (d < 1) ? d^3 - 2d^2 + 1 : (d < 2) ? -d^3 + 5d^2 - 8d + 4 : 0
                end
                # calc coefficients
                xd= xin - floor(xin)
                xc= map(x -> fsinc(x), [xd+1; xd; 1-xd; 2-xd])
                #
                yd= yin - floor(yin)
                yc= map(x -> fsinc(x), [yd+1 yd 1-yd 2-yd])
                #
                return xc * yc
        end

        pixels= getPIX4x4(ibuff, xin, yin)
        coeffs= getCOEFFS(xin, yin)

        return sum(pixels .* coeffs) / sum(coeffs)

end