Бикубическая интерполяция?

Я просмотрел Интернет и с точки зрения бикубической интерполяции не могу найти для нее простого уравнения. Страница Википедии по этому вопросу не очень полезна, поэтому есть ли простой способ узнать, как работает бикубическая интерполяция и как ее реализовать? Я использую его для генерации шума Перлина, но использование билинейной интерполяции слишком прерывисто для моих нужд (я уже пробовал).

Если кто-нибудь может указать мне правильное направление с помощью хорошего веб-сайта или просто ответа, я был бы очень признателен. (кстати, я использую С#)


person Nicholas Pipitone    schedule 04.01.2014    source источник
comment
Вы проверили github о Perlin Noise? Я надеюсь, что это может быть полезно: github.com/marbsydo/PerlinNoise/blob/master. /Noise.cs   -  person Ahmet Kakıcı    schedule 04.01.2014
comment
@AhmetKakıcı Большое спасибо! Кажется, мне потребовалось около 1/2 часа, чтобы запустить его, потому что их код сильно отличался от моего, но я, наконец, понял это. :) Спасибо за помощь.   -  person Nicholas Pipitone    schedule 04.01.2014
comment
Приятно слышать, что вы это сделали, и не забудьте отметить свой ответ как правильный, чтобы помочь новичкам ;)   -  person Ahmet Kakıcı    schedule 04.01.2014


Ответы (4)


Используя это (спасибо Ахмету Какычи, который нашел это), я понял как добавить бикубическую интерполяцию. Для тех, кто также ищет ответ, вот что я использовал:

private float CubicPolate( float v0, float v1, float v2, float v3, float fracy ) {
    float A = (v3-v2)-(v0-v1);
    float B = (v0-v1)-A;
    float C = v2-v0;
    float D = v1;

    return A*Mathf.Pow(fracy,3)+B*Mathf.Pow(fracy,2)+C*fracy+D;
}

Чтобы получить 2D-интерполяцию, я сначала получил x, а затем интерполировал y. Например.

float x1 = CubicPolate( ndata[0,0], ndata[1,0], ndata[2,0], ndata[3,0], fracx );
float x2 = CubicPolate( ndata[0,1], ndata[1,1], ndata[2,1], ndata[3,1], fracx );
float x3 = CubicPolate( ndata[0,2], ndata[1,2], ndata[2,2], ndata[3,2], fracx );
float x4 = CubicPolate( ndata[0,3], ndata[1,3], ndata[2,3], ndata[3,3], fracx );

float y1 = CubicPolate( x1, x2, x3, x4, fracy );

Где ndata определяется как:

float[,] ndata = new float[4,4];
for( int X = 0; X < 4; X++ )
    for( int Y = 0; Y < 4; Y++ )
        //Smoothing done by averaging the general area around the coords.
        ndata[X,Y] = SmoothedNoise( intx+(X-1), inty+(Y-1) );

(intx и inty — это полные значения запрошенных координат. fracx и fracy — дробные части введенных координат, равные x-intx и y-inty соответственно)

person Nicholas Pipitone    schedule 04.01.2014
comment
Привет, Николай! Не могли бы вы присоединиться ко мне в этом вопросе? - person netaholic; 19.08.2015
comment
A*Mathf.Pow(fracy,3)+B*Mathf.Pow(fracy,2)+C*fracy+D; неэффективен с точки зрения времени и точности. Предлагаю реализовать таким образом: D + fracy * (C + fracy * (B + fracy * A)) - person Angelo Mascaro; 22.02.2016
comment
@Nicholas Pipitone, не могли бы вы объяснить нам, что такое десятичные дроби введенных координат, например, координата (x, y) = (1,2) - person VenushkaT; 22.06.2016
comment
@VenushkaT Хм, я, видимо, сделал правильные правки в сентябре 16 года, но не сообщил вам об этом. Это (очень) поздний пинг, но мой последний комментарий теперь должен быть более читабельным. - person Nicholas Pipitone; 29.09.2018

Принял ответ Эске Рана и сделал один вызов (обратите внимание, в приведенном ниже коде используется соглашение о размерах матрицы (j, i), а не изображение (x, y), но это не должно иметь значения для интерполяции):

/// <summary>
/// Holds extension methods.
/// </summary>
public static class Extension
{
    /// <summary>
    /// Performs a bicubic interpolation over the given matrix to produce a
    /// [<paramref name="outHeight"/>, <paramref name="outWidth"/>] matrix.
    /// </summary>
    /// <param name="data">
    /// The matrix to interpolate over.
    /// </param>
    /// <param name="outWidth">
    /// The width of the output matrix.
    /// </param>
    /// <param name="outHeight">
    /// The height of the output matrix.
    /// </param>
    /// <returns>
    /// The interpolated matrix.
    /// </returns>
    /// <remarks>
    /// Note, dimensions of the input and output matrices are in
    /// conventional matrix order, like [matrix_height, matrix_width],
    /// not typical image order, like [image_width, image_height]. This
    /// shouldn't effect the interpolation but you must be aware of it
    /// if you are working with imagery.
    /// </remarks>
    public static float[,] BicubicInterpolation(
        this float[,] data, 
        int outWidth, 
        int outHeight)
    {
        if (outWidth < 1 || outHeight < 1)
        {
            throw new ArgumentException(
                "BicubicInterpolation: Expected output size to be " +
                $"[1, 1] or greater, got [{outHeight}, {outWidth}].");
        }

        // props to https://stackoverflow.com/a/20924576/240845 for getting me started
        float InterpolateCubic(float v0, float v1, float v2, float v3, float fraction)
        {
            float p = (v3 - v2) - (v0 - v1);
            float q = (v0 - v1) - p;
            float r = v2 - v0;

            return (fraction * ((fraction * ((fraction * p) + q)) + r)) + v1;
        }

        // around 6000 gives fastest results on my computer.
        int rowsPerChunk = 6000 / outWidth; 
        if (rowsPerChunk == 0)
        {
            rowsPerChunk = 1;
        }

        int chunkCount = (outHeight / rowsPerChunk) 
                         + (outHeight % rowsPerChunk != 0 ? 1 : 0);

        var width = data.GetLength(1);
        var height = data.GetLength(0);
        var ret = new float[outHeight, outWidth];

        Parallel.For(0, chunkCount, (chunkNumber) =>
        {
            int jStart = chunkNumber * rowsPerChunk;
            int jStop = jStart + rowsPerChunk;
            if (jStop > outHeight)
            {
                jStop = outHeight;
            }

            for (int j = jStart; j < jStop; ++j)
            {
                float jLocationFraction = j / (float)outHeight;
                var jFloatPosition = height * jLocationFraction;
                var j2 = (int)jFloatPosition;
                var jFraction = jFloatPosition - j2;
                var j1 = j2 > 0 ? j2 - 1 : j2;
                var j3 = j2 < height - 1 ? j2 + 1 : j2;
                var j4 = j3 < height - 1 ? j3 + 1 : j3;
                for (int i = 0; i < outWidth; ++i)
                {
                    float iLocationFraction = i / (float)outWidth;
                    var iFloatPosition = width * iLocationFraction;
                    var i2 = (int)iFloatPosition;
                    var iFraction = iFloatPosition - i2;
                    var i1 = i2 > 0 ? i2 - 1 : i2;
                    var i3 = i2 < width - 1 ? i2 + 1 : i2;
                    var i4 = i3 < width - 1 ? i3 + 1 : i3;
                    float jValue1 = InterpolateCubic(
                        data[j1, i1], data[j1, i2], data[j1, i3], data[j1, i4], iFraction);
                    float jValue2 = InterpolateCubic(
                        data[j2, i1], data[j2, i2], data[j2, i3], data[j2, i4], iFraction);
                    float jValue3 = InterpolateCubic(
                        data[j3, i1], data[j3, i2], data[j3, i3], data[j3, i4], iFraction);
                    float jValue4 = InterpolateCubic(
                        data[j4, i1], data[j4, i2], data[j4, i3], data[j4, i4], iFraction);
                    ret[j, i] = InterpolateCubic(
                        jValue1, jValue2, jValue3, jValue4, jFraction);
                }
            }
        });

        return ret;
    }
}
person mheyman    schedule 08.04.2019

Я немного запутался в использовании многочлена третьей степени.

Да, он дает правильные значения в 0 и 1, но производные соседних ячеек не подходят, насколько я могу вычислить. Если данные сетки линейны, они даже не возвращают строку....

И это не точечная симметрия в x = 0,5

Многочлен, который соответствует 0 и 1 И также имеет те же производные для соседних ячеек и, следовательно, является гладким, (почти) так же легко вычисляется.

(и сводится к линейной форме, если это соответствует данным)

В маркировке p и m - это сокращение от плюса и минуса, например. vm1 это v(-1)

    //Bicubic convolution algorithm, cubic Hermite spline
    static double CubicPolateConv
            (double vm1, double v0, double vp1, double vp2, double frac) {
        //The polynomial of degree 3 where P(x)=f(x) for x in {0,1}
        //and P'(1) in one cell matches P'(0) in the next, gives a continuous smooth curve.
        //And we also wants it to reduce nicely to a line, if that matches the data
        //P(x)=Dx³+Cx²+Bx+A=((Dx+C)x+B)x+A
        //P'(x)=3Dx²+2Cx+B
        //P(0)=A       =v0
        //P(1)=D+C+B+A =Vp1
        //P'(0)=B      =(vp1-vm1)/2
        //P'(1)=3D+2C+B=(vp2-v0 )/2
        //Subtracting expressions for A and B from D+C+B+A  
        //D+C =vp1-B-A = (vp1+vm1)/2 - v0
        //Subtracting that twice and a B from the P'(1)
        //D=(vp2-v0)/2 - 2(D+C) -B =(vp2-v0)/2 - (Vp1+vm1-2v0) - (vp1-vm1)/2
        // = 3(v0-vp1)/2 + (vp2-vm1)/2
        //C=(D+C)-D = (vp1+vm1)/2 - v0 - (3(v0-vp1)/2 + (vp2-vm1)/2)
        // = vm1 + 2vp1 - (5v0+vp2)/2;
        //It is quite easy to calculate P(½)
        //P(½)=D/8+C/4+B/2+A = (9*(v0+vp1)-(vm1+vp2))/16
        //i.e. symmetric in its uses, and a mean of closest adjusted by mean of next ones

        double B = (vp1 - vm1) / 2;
        double DpC =(vp1 -v0) -B; //D+C+B+A - A - B
        double D = (vp2 - v0) / 2 - 2 * DpC - B;
        double C = DpC - D;
        //double C = vm1 + 2 * vp1 - (5 * v0 + vp2) / 2;
        //double D = (3*(v0 - vp1) + (vp2 - vm1)) / 2;

        return ((D * frac + C) * frac + B) * frac + A;
    }

Вдохновленный комментарием ManuTOO ниже, я сделал это также как четвертый порядок, с необязательным параметром, который вы можете установить/вычислить по своему усмотрению, не нарушая плавность кривой. Это в основном то же самое с добавленным членом на всем протяжении расчетов. И если вы установите E на ноль, это будет идентично приведенному выше. (Очевидно, что если данные на самом деле находятся в строке, ваш расчет E должен гарантировать, что E равно нулю, чтобы получить линейный вывод)

        //The polynomial of degree 4 where P(x)=f(x) for x in {0,1}
        //and P'(1) in one cell matches P'(0) in the next, gives a continuous smooth curve.
        //And we also wants the it to reduce nicely to a line, if that matches the data
        //With high order quotient weight of your choice....
        //P(x)=Ex⁴+Dx³+Cx²+Bx+A=((((Ex+D)x+C)x+B)x+A
        //P'(x)=4Ex³+3Dx²+2Cx+B
        //P(0)=A          =v0
        //P(1)=E+D+C+B+A  =Vp1
        //P'(0)=B         =(vp1-vm1)/2
        //P'(1)=4E+3D+2C+B=(vp2-v0 )/2
        //Subtracting Expressions for A, B and E from the E+D+C+B+A 
        //D+C =vp1-B-A -E = (vp1+vm1)/2 - v0 -E
        //Subtracting that twice, a B and 4E from the P'(1)
        //D=(vp2-v0)/2 - 2(D+C) -B -4E =(vp2-v0)/2 - (Vp1+vm1-2v0-2E) - (vp1-vm1)/2 -4E
        // = 3(v0-vp1)/2 + (vp2-vm1)/2 -2E
        //C=(D+C)-(D) = (vp1+vm1)/2 - v0 -E - (3(v0-vp1)/2 + (vp2-vm1)/2 -2E)
        // = vm1 + 2vp1 - (5v0+vp2)/2 +E;

        double E = ????????; //Your feed.... If Zero, cubic, see below
        double B = (vp1 - vm1) / 2;
        double DpC =(vp1 -v0) -B -E; //E+D+C+B+A - A - B -E
        double D = (vp2 - v0) / 2 - 2 * DpC - B - 4 * E;
        double C = DpC - D;

        return (((E * frac + D) * frac + C) * frac + B) * frac + v0;

ДОБАВИТЬ: Цитата из предложения @ManTOO ниже:

         double E = (v0 - vm1 + vp1 - vp2) * m_BicubicSharpness;

При значении m_BicubicSharpness 1,5 это очень близко к Bicubic Sharper в Photoshop; лично я установил его на 1,75 для небольшого дополнительного удара.

Обратите внимание, что если данные находятся в строке, это предложение сводится к нулю.

person Eske Rahn    schedule 12.02.2017
comment
Интерполяция - это плохо поставленная проблема, поэтому я думаю, что проблема заключается в том, что работает лучше. Если бы это было сглаживание, мы могли бы сказать, что это работает лучше или это. Идея состоит в том, чтобы подогнать функцию к точкам, но вы не знаете исходный сигнал, поэтому, так сказать, нет наилучшего соответствия. - person Minato; 13.02.2017
comment
Спасибо за отзыв. Я отредактировал предложение, чтобы получить непрерывную и гладкую кривую, которая также сводится к линии (если данные соответствуют линии). - и включил доказательство. - person Eske Rahn; 13.02.2017
comment
@Eske Rahn, спасибо за улучшенный метод! Вы случайно не знаете, можно ли его настроить, чтобы сделать результат немного четче? (аналогично бикубической резкости в Photoshop) - person ManuTOO; 16.06.2020
comment
@ManuTOO Извините. Я представил только простое математическое гладкое решение низкого порядка. Если вы предпочитаете, чтобы он был более резким, вы можете использовать полином более высокого порядка, чтобы получить более крутые края. - person Eske Rahn; 17.06.2020
comment
@ Эске Ран, добавление порядка полинома даст мне новую переменную, которую я мог бы настроить, чтобы сделать все более плавным или четким, верно? (Мне не нравятся английские термины в математике, поэтому я просто хочу убедиться, что правильно понял :-)) - person ManuTOO; 18.06.2020
comment
@ManuTOO, да, но вы, скорее всего, добавили бы два, чтобы получить пятый порядок, чтобы он оставался нечетным. Это /может быть/ не нужно, но моя интуиция подсказывает мне, что это должно быть странно... - person Eske Rahn; 19.06.2020
comment
@ Эске Ран, я только что сделал это с 4-м заказом, и он работает хорошо. В очередной раз благодарим за помощь ! :-) - person ManuTOO; 20.06.2020
comment
@ManuTOO, изменил ответ, включая четвертый порядок, - person Eske Rahn; 20.06.2020
comment
Для резкости используйте это: double E = (v0 - vm1 + vp1 - vp2) * m_BicubicSharpness; При значении m_BicubicSharpness 1,5 это очень близко к Bicubic Sharper в Photoshop; лично я установил его на 1,75 для дополнительного удара. Я хотел бы предположить, что он становится более гладким, когда он отрицательный, но у меня есть некоторые сомнения по этому поводу. :-п - person ManuTOO; 21.06.2020

Примечание. Бикубическая формула, данная Николаем выше (в качестве ответа), содержит ошибку. интерполируя синусоиду, я смог увидеть, что она неверна.

Правильная формула:

float A = 0.5f * (v3 - v0) + 1.5f * (v1 - v2);
float B = 0.5f * (v0 + v2) - v1 - A;
float C = 0.5f * (v2 - v0);
float D = v1;

Для получения информации см. https://www.paulinternet.nl/?page=bicubic.

person BobO    schedule 27.08.2020