Massive artifacting in bicubic interpolation; how to fix?

Question

I'm trying to implement a bicubic interpolation algorithm to reconstruct higher-resolution data from a heightmap. After some falstarts and sets of instructions containing nearly-incomprehensible math (it's been some years since I had calculus and I no longer remember much beyond the basics) , I've found the article "Bicubic Interpolation for Image Scaling" by Paul Bourke containing what appeared to be a fairly simple and easy to implement algorithm. http://paulbourke.net/texture_colour/imageprocess/

However, instead of producing an interpolation result even remotely resembling the one on Wikipedia, instead I'm getting this (from the same input data):

What is causing the error, and more importantly - how to fix it?

PHP code below (yes, this probably should - and will - be reimplemented in C; when it's working)

class BicubicInterpolator
{
    private $data;
    public function Set_data($d)
    {
        $this->data=$this->denull($d);
    }
    public function Interpolate($dx,$dy)
    {   
        $r=0;
        for ($m=-1; $m<2; $m++)
            for ($n=-1; $n<2; $n++)
                $r+=$this->data[$m+1][$n+1] * $this->R($m-$dx) * $this->R($dy-$n);
        return $r;
    }
    private function denull($d)
    {
        //Substituting null values with nearest known values as per "A Review of Some Image Pixel Interpolation Algorithms" by Don Lancaster (supposed to produce same output as example image)
        if ($d[0][1]===null) for ($i=0; $i<4; $i++) $d[0][$i]=$d[1][$i];
        if ($d[1][0]===null) for ($i=0; $i<4; $i++) $d[$i][0]=$d[$i][1];
        if ($d[3][1]===null) for ($i=0; $i<4; $i++) $d[3][$i]=$d[2][$i];
        if ($d[1][3]===null) for ($i=0; $i<4; $i++) $d[$i][3]=$d[$i][2];
        return $d;
    }
    function R($x)
    {
        return (      $this->P($x+2)
            - 4 * $this->P($x+1)
            + 6 * $this->P($x)
            - 4 * $this->P($x-1) )/6;
    }
    function P($x)
    {
        if ($x>0) return $x*$x*$x;
        return 0;
    }

Answer 1

In the end I have switched to a different algorithm, based on a combination of the one outlined by Don Lancaster with derivative and cross-derivative formulae from Chapter 3.6 of "Numerical Recipes in C" (2nd Edition), p136.

This is combined with two minor tweaks:

1. The function calculation caches an intermediary set of values for last y-coordinate (a successive Interpolate(x,y) call with the same y argument will require four multiplications and three additions, cutting down on processing time)
2. One set of derivatives is accessible as public, allowing the last two to be passed along as first two for a grid cell in coordinates (x, y+1) and halving the amount of calculations needed to get each of those derivative sets for each cell.

This is the implementation, working with no apparent glitches:

class BicubicInterpolator
{
    private $last_y;
    private $last_y_a;
    private $a;
    public $x;
    public function __construct()
    {
        for ($i=0;$i<4;$i++)
            $this->x[$i]=false;
    }
    public function Set_data($d)
    {
        $d=$this->denull($d);
        $x=$this->x;
        for ($j=1; $j<3; $j++)
            for ($k=1; $k<3; $k++)
            {
                $r=($j-1)*2+($k-1);
                $w[$r]=$d[$j][$k];
                //Derivatives and cross derivatives calculated as per Numerical Recipes in C, 2nd edition.
                if (!$x[$r]) $x[$r]=( $d[$j][$k+1] - $d[$j][$k-1] ) / 2;
                $y[$r]=( $d[$j+1][$k] - $d[$j-1][$k] ) / 2;
                $z[$r]=( $d[$j+1][$k+1]-$d[$j+1][$k-1]-$d[$j-1][$k+1]+$d[$j-1][$k-1] )/4;
            }
        $this->x=$x;
        /* Coefficient calculation as per "A Review of Some Image Pixel Interpolation Algorithms" by Don Lancaster, 
        + addressing changed to (x,y) instead of (y,x)
        + reformulated to minimize the number of multiplications required */
        $this->a[0][0] = $w[0];
        $this->a[1][0] = $y[0];
        $this->a[2][0] = 3*($w[2]-$w[0])-2*$y[0]-$y[2];
        $this->a[3][0] = 2*($w[0]-$w[2])+$y[0]+$y[2];
        $this->a[0][1] = $x[0];
        $this->a[1][1] = $z[0];
        $this->a[2][1] = 3*($x[2]-$x[0])-2*$z[0]-$z[2];
        $this->a[3][1] = 2*($x[0]-$x[2])+$z[0]+$z[2];
        $this->a[0][2] = 3*($w[1]-$w[0])-2*$x[0]-$x[1];
        $this->a[1][2] = 3*($y[1]-$y[0])-2*$z[0]-$z[1];
        $this->a[2][2] = 9*($w[0]-$w[1]-$w[2]+$w[3])+6*($x[0]-$x[2]+$y[0]-$y[1])+3*($x[1]-$x[3]+$y[2]-$y[3])+4*$z[0]+2*($z[1]+$z[2])+$z[3];
        $this->a[3][2] = 6*($w[1]+$w[2]-$w[3]-$w[0])+4*($x[2]-$x[0])+3*($y[1]-$y[0]-$y[2]+$y[3])+2*($x[3]-$z[0]-$z[2]-$x[1])-$z[1]-$z[3];
        $this->a[0][3] = 2*($w[0]-$w[1])+$x[0]+$x[1];
        $this->a[1][3] = 2*($y[0]-$y[1])+$z[0]+$z[1];
        $this->a[2][3] = 6*($w[1]+$w[2]-$w[0]-$w[3])+3*(-$x[0]-$x[1]+$x[2]+$x[3])+4*($y[1]-$y[0])+2*($y[3]-$y[2]-$z[0]-$z[1])-$z[2]-$z[3];
        $this->a[3][3] = 4*($w[0]-$w[1]-$w[2]+$w[3])+2*($x[0]+$x[1]-$x[2]-$x[3]+$y[0]-$y[1]+$y[2]-$y[3])+$z[0]+$z[1]+$z[2]+$z[3];

        $this->last_y=false;
    }
    public function Interpolate($x,$y)
    {
        if ($y!==$this->last_y)
        {
            for ($i=0; $i<4; $i++)
            {
                $this->last_y_a[$i]=0;
                for ($j=0; $j<4; $j++)
                    $this->last_y_a[$i]+=$this->a[$j][$i]*pow($y,$j);
            }
            $this->last_y=$y;
        }
        $r=0;
        for ($i=0; $i<4; $i++)
            $r+=$this->last_y_a[$i]*pow($x,$i);
        return $r;
    }
    private function denull($d)
    {
        //Substituting null values with nearest known values 
        //as per "A Review of Some Image Pixel Interpolation Algorithms" by Don Lancaster
        if ($d[0][1]===null)    for ($i=0; $i<4; $i++)  $d[0][$i]=$d[1][$i];
        if ($d[1][0]===null)    for ($i=0; $i<4; $i++)  $d[$i][0]=$d[$i][1];
        if ($d[3][1]===null)    for ($i=0; $i<4; $i++)  $d[3][$i]=$d[2][$i];
        if ($d[1][3]===null)    for ($i=0; $i<4; $i++)  $d[$i][3]=$d[$i][2];
        return $d;
    }
}