Массивные артефакты при бикубической интерполяции; как исправить?

Я пытаюсь реализовать алгоритм бикубической интерполяции для восстановления данных с более высоким разрешением из карты высот. После некоторых фальстартов и наборов инструкций, содержащих почти непонятную математику (прошло несколько лет с тех пор, как я занимался исчислением, и я больше не помню ничего, кроме основ), я нашел статью «Бикубическая интерполяция для масштабирования изображения» Пола Бурка, содержащую то, что оказался довольно простым и легко реализуемым алгоритмом. http://paulbourke.net/texture_colour/imageprocess/

Однако вместо получения результата интерполяции, даже отдаленно напоминающего тот, что в Википедии, вместо этого я получаю это (из тех же входных данных):

введите здесь описание изображения

Что вызывает ошибку, и что более важно - как ее исправить?

PHP-код ниже (да, это, вероятно, должно быть и будет реализовано на C, когда оно заработает)

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

person Michał Gawlas    schedule 17.09.2012    source источник
comment
Вот пример из Википедии, о котором я говорю: upload.wikimedia.org/wikipedia/commons/thumb/d/d5/   -  person Michał Gawlas    schedule 17.09.2012


Ответы (1)


В конце концов я переключился на другой алгоритм, основанный на комбинации алгоритма, описанного Доном Ланкастером, с формулами производных и перекрестных производных из главы 3.6 «Численных рецептов в C» (2-е издание), стр. 136.

Это сочетается с двумя небольшими изменениями:

  1. Вычисление функции кэширует промежуточный набор значений для последней координаты y (последовательный вызов Interpolate(x,y) с тем же аргументом y потребует четырех умножений и трех сложений, что сокращает время обработки)
  2. Один набор производных доступен как общедоступный, что позволяет передавать последние два как первые два для ячейки сетки с координатами (x, y+1) и вдвое сокращает объем вычислений, необходимых для получения каждого из этих наборов производных для каждой ячейки. .

Это реализация, работающая без явных глюков:

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;
    }
}
person Michał Gawlas    schedule 18.09.2012
comment
что именно эта функция принимает в качестве параметра? Как я должен представить изображение? - person Sami; 18.03.2013
comment
Он принимает массив значений 4x4, представляющий значения сетки в 16 точках, четыре внутренних из которых отмечают углы области, значения которой вы будете получать — имейте в виду, что адресация, используемая Interpolate(), находится в диапазоне 0 к 1. Интерполяция сеток размером более 4x4 может быть выполнена путем их разделения (действительно, это то, что я сделал тогда). Он был разработан для интерполяции числовых данных, которые затем были преобразованы в изображения в искусственных цветах. Можно выполнить интерполяцию реальных изображений, но вам придется разделить цветовые каналы на отдельные массивы. - person Michał Gawlas; 07.08.2013