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;
}
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:
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;
}
}