Правильная реализация интерполяции кубическим сплайном

Я использовал один из предложенных алгоритмов, но результаты очень плохие.

Я реализовал алгоритм вики

на Java (код ниже). x(0) - это points.get(0), y(0) - это values[points.get(0)], α - это alfa, а μ - это mi. В остальном то же, что и в псевдокоде вики.

    public void createSpline(double[] values, ArrayList<Integer> points){
    a = new double[points.size()+1];

    for (int i=0; i <points.size();i++)
    {
        a[i] = values[points.get(i)];

    }

    b = new double[points.size()];
    d = new double[points.size()];
    h = new double[points.size()];

    for (int i=0; i<points.size()-1; i++){
        h[i] = points.get(i+1) - points.get(i);
    }

    alfa = new double[points.size()];

    for (int i=1; i <points.size()-1; i++){
        alfa[i] = (double)3 / h[i] * (a[i+1] - a[i]) - (double)3 / h[i-1] *(a[i+1] - a[i]);
    }

    c = new double[points.size()+1];
    l = new double[points.size()+1];
    mi = new double[points.size()+1];
    z = new double[points.size()+1];

    l[0] = 1; mi[0] = z[0] = 0;

    for (int i =1; i<points.size()-1;i++){
        l[i] = 2 * (points.get(i+1) - points.get(i-1)) - h[i-1]*mi[i-1];
        mi[i] = h[i]/l[i];
        z[i] = (alfa[i] - h[i-1]*z[i-1])/l[i];
    }

    l[points.size()] = 1;
    z[points.size()] = c[points.size()] = 0;

    for (int j=points.size()-1; j >0; j--)
    {
        c[j] = z[j] - mi[j]*c[j-1];
        b[j] = (a[j+1]-a[j]) - (h[j] * (c[j+1] + 2*c[j])/(double)3) ;
        d[j] = (c[j+1]-c[j])/((double)3*h[j]);
    }


    for (int i=0; i<points.size()-1;i++){
        for (int j = points.get(i); j<points.get(i+1);j++){
            //                fk[j] = values[points.get(i)];
            functionResult[j] = a[i] + b[i] * (j - points.get(i)) 
                                + c[i] * Math.pow((j - points.get(i)),2)
                                + d[i] * Math.pow((j - points.get(i)),3);
        }
    }

}

Результат, который я получаю, следующий:

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

но это должно быть похоже на это:

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


Я также пытаюсь реализовать алгоритм по-другому в соответствии с: http://www.geos.ed.ac.uk/~yliu23/docs/lect_spline.pdf

Сначала они показывают, как делать линейные сплайны, и это довольно просто. Я создаю функции, которые вычисляют коэффициенты A и B. Затем они расширяют линейный сплайн, добавляя вторую производную. Коэффициенты C и D тоже легко вычислить.

Но проблемы начинаются, когда я пытаюсь вычислить вторую производную. Я не понимаю, как они их рассчитывают.

Поэтому я реализовал только линейную интерполяцию. Результат:

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

Кто-нибудь знает, как исправить первый алгоритм или объяснить мне, как вычислить вторую производную во втором алгоритме?


person Mr Jedi    schedule 30.11.2013    source источник
comment
Попробуйте здесь: codereview.stackexchange.com   -  person ceving    schedule 30.11.2013
comment
@ceving они говорят мне, что проверяют качество кода, а не результаты кода.   -  person Mr Jedi    schedule 30.11.2013
comment
@ceving Поскольку результаты здесь не являются желаемыми, этот вопрос не подходит для проверки кода.   -  person Simon Forsberg    schedule 30.11.2013
comment
будет приятно увидеть контрольные точки на ваших графиках, чтобы увидеть, что на самом деле не так, а что нет   -  person Spektre    schedule 11.12.2013
comment
Возможно, не стоит использовать устаревшие статьи в Википедии, помеченные как запутывающие. Недавняя статья о сплайн-интерполяции, к удивлению, находится на странице сплайн-интерполяции.   -  person Lutz Lehmann    schedule 18.12.2013


Ответы (4)


Кубический b-сплайн недавно был описан в серии статей Unser, Thévenaz et al., см. Среди других

М. Унзер, А. Альдроуби, М. Иден, Быстрые B-сплайновые преобразования для непрерывного представления изображений и интерполяции, IEEE Trans. Pattern Anal. Машинный интеллект., т. 13, п. 3, стр. 277-285, март 1991 г.

М. Унсер, Сплайны, идеально подходят для обработки сигналов и изображений, IEEE Signal Proc. Mag., стр. 22–38, ноябрь 1999 г.

и

П. Тевеназ, Т. Блу, М. Унсер, Возвращение к интерполяции, IEEE Trans. по медицинской визуализации, т. 19, нет. 7, pp. 739-758, июль 2000 г.

Вот несколько рекомендаций.

Что такое шлицы?

Сплайны - это кусочно-полиномы, которые плавно соединены между собой. Для сплайна степени n каждый сегмент является полиномом степени n. Куски соединены так, чтобы сплайн был непрерывным до своей производной степени n-1 в узлах, а именно в точках соединения полиномиальных частей.

Как можно построить сплайны?

Сплайн нулевого порядка имеет вид

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

Все остальные шлицы можно построить как

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

где свертка сделана n-1 раз.

Кубические шлицы

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

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

Проблема интерполяции сплайна

Учитывая функцию f(x), выбранную в дискретных целочисленных точках k, задача сплайн-интерполяции состоит в том, чтобы определить приближение s(x) к f(x), выраженное следующим образом

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

где ck - коэффициенты интерполяции, а s(k) = f(k).

Предварительная фильтрация сплайнов

К сожалению, начиная с n=3 г.,

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

так что ck не являются коэффициентами интерполяции. Их можно определить, решив линейную систему уравнений, полученную принуждением s(k) = f(k), а именно:

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

Такое уравнение можно преобразовать в свертку и решить в преобразованном z-пространстве как

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

где

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

Соответственно,

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

Такой путь всегда предпочтительнее, чем решение линейной системы уравнений, например, путем разложения.

Решение вышеуказанного уравнения можно определить, заметив, что

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

где

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

Первая фракция представляет причинный фильтр, а вторая - антикаузальный фильтр. Оба они показаны на рисунках ниже.

Причинный фильтр

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

Антикаузальный фильтр

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

На последнем рисунке

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

Выходные данные фильтров можно выразить следующими рекурсивными уравнениями

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

Вышеупомянутые уравнения могут быть решены путем определения начальных условий для c- и c+. При условии периодической зеркальной входной последовательности fk такой, что

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

тогда можно показать, что начальное условие для c+ может быть выражено как

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

в то время как начальное условие для c+ может быть выражено как

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

person Vitality    schedule 14.12.2014

Извините, но ваш исходный код для меня действительно нечитаемый, поэтому я придерживаюсь теории. Вот несколько советов:

  1. SPLINE кубики

    СПЛАЙН - это не интерполяция, а приближение, чтобы использовать их, вам не нужно никаких производных. Если у вас десять точек: p0,p1,p2,p3,p4,p5,p6,p7,p8,p9, тогда кубический сплайн начинается / заканчивается тройной точкой. Если вы создадите функцию для «рисования» участка кубической кривой SPLINE, тогда для обеспечения непрерывности последовательность вызовов будет такой:

        spline(p0,p0,p0,p1);
        spline(p0,p0,p1,p2);
        spline(p0,p1,p2,p3);
        spline(p1,p2,p3,p4);
        spline(p2,p3,p4,p5);
        spline(p3,p4,p5,p6);
        spline(p4,p5,p6,p7);
        spline(p5,p6,p7,p8);
        spline(p6,p7,p8,p9);
        spline(p7,p8,p9,p9);
        spline(p8,p9,p9,p9);
    

    не забывайте, что кривую SPLINE для p0,p1,p2,p3 рисуйте только кривую "между" p1,p2 !!!

  2. Кубики БЕЗЬЕ

    4-точечные коэффициенты Безье можно вычислить следующим образом:

        a0=                              (    p0);
        a1=                    (3.0*p1)-(3.0*p0);
        a2=          (3.0*p2)-(6.0*p1)+(3.0*p0);
        a3=(    p3)-(3.0*p2)+(3.0*p1)-(    p0);
    
  3. Интерполяция

    чтобы использовать интерполяцию, вы должны использовать полиномы интерполяции. Их много, но я предпочитаю использовать кубики ... похожие на BEZIER / SPLINE / NURBS ..., поэтому

    • p(t) = a0+a1*t+a2*t*t+a3*t*t*t where t = <0,1>

    Осталось только вычислить a0,a1,a2,a3. У вас есть 2 уравнения (p(t) и его вывод по t) и 4 точки из набора данных. Вы также должны обеспечить непрерывность ... Таким образом, первая производная для точек соединения должна быть одинаковой для соседних кривых (t=0,t=1). Это приводит к системе 4 линейных уравнений (используйте t=0 и t=1). Если вы его выведете, оно создаст простое уравнение, зависящее только от координат входной точки:

        double  d1,d2;
        d1=0.5*(p2-p0);
        d2=0.5*(p3-p1);
        a0=p1;
        a1=d1;
        a2=(3.0*(p2-p1))-(2.0*d1)-d2;
        a3=d1+d2+(2.0*(-p2+p1));
    

    последовательность вызовов такая же, как для SPLINE

[Примечания]

  1. разница между интерполяцией и аппроксимацией заключается в следующем:

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

  2. переменные:

    • a0,... p0,... are vectors (number of dimensions must match the input points)
    • t скалярный
  3. рисовать кубик из коэффициентов a0,..a3

    просто сделайте что-нибудь вроде этого:

        MoveTo(a0.x,a0.y);
        for (t=0.0;t<=1.0;t+=0.01)
         {
         tt=t*t;
         ttt=tt*t;
         p=a0+(a1*t)+(a2*tt)+(a3*ttt);
         LineTo(p.x,p.y);
         }
    
person Spektre    schedule 11.12.2013
comment
Можно выполнять интерполяцию кубических сплайнов, где значения функций, производные и вторые производные совпадают в точках интерполяции. Просто нужно решить большую, но ленточную систему линейных уравнений. - person Lutz Lehmann; 18.12.2013
comment
зачем делать это в обратном направлении, кубическая кривая имеет предел связности c1 (значение и первое совпадение деривации), если требуется более высокая связность, тогда используйте больше точек 5 точек - ›c2, 6 точек -› c3 и так далее ... с 4 точками для достижения связности c2 это больше не СПЛАЙН ... и даже не 4-х точечный куб, потому что вы также решаете близкие точки ... Конечно, я могу пропустить что-то вроде новых уравнений, обнаруженных в недавнем прошлом, но я сильно сомневаюсь в этом - person Spektre; 18.12.2013
comment
Я считаю, что вы пропустили несколько очень старых формул. Да, по значениям и производным вы можете построить кусочно-кубическую функцию. Если у вас нет данных о производных, то есть большая свобода в их выборе. Эта свобода ограничивается требованием минимизированных вторых производных. И действительно, это приводит к дважды дифференцируемой кусочно-кубической с линейными концами, которая однозначно определяется трехдиагональной линейной системой. См. сплайн-интерполяцию - person Lutz Lehmann; 18.12.2013
comment
И обычно термин сплайн зарезервирован для тех кусочно-кубических функций, которые имеют минимальную кривизну (или вторую производную) для данных ограничений. С кубическими сплайнами получается не лучше, чем C2, третья производная почти всегда кусочно-постоянная со скачками. Если скачков нет, то третья производная постоянна, а функция - кубический многочлен. - person Lutz Lehmann; 18.12.2013

См. сплайн-интерполяцию, хотя они дают только полезный пример 3x3. Для получения дополнительных точек выборки, скажем, N + 1 пронумерованных x[0..N] со значениями y[0..N], нужно решить следующую систему для неизвестного k[0..N]

              2*    c[0]    * k[0] + c[0] * k[1] == 3*   d[0];
c[0] * k[0] + 2*(c[0]+c[1]) * k[1] + c[1] * k[2] == 3*(d[0]+d[1]);
c[1] * k[1] + 2*(c[1]+c[2]) * k[2] + c[2] * k[3] == 3*(d[1]+d[2]);
c[2] * k[2] + 2*(c[2]+c[3]) * k[3] + c[3] * k[4] == 3*(d[2]+d[3]);
...
c[N-2]*k[N-2]+2*(c[N-2]+c[N-1])*k[N-1]+c[N-1]*k[N] == 3*(d[N-2]+d[N-1]);
c[N-2]*k[N-1] + 2*c[N-1]        *  k[N]          == 3*   d[N-1];

где

c[k]=1/(x[k+1]-x[k]); d[k]=(y[k+1]-y[k])*c[k]*c[k];

Это можно решить, используя итерацию Гауса-Зейделя (разве она не была изобретена именно для этой системы?) Или ваш любимый решатель пространства Крылова.

person Lutz Lehmann    schedule 18.12.2013

// Файл: cbspline-test.cpp

// C++ Cubic Spline Interpolation using the Boost library.
// Interpolation is an estimation of a value within two known values in a sequence of values. 

// From a selected test function y(x) = (5.0/(x))*sin(5*x)+(x-6), 21 numbers of (x,y) equal-spaced sequential points were generated. 
// Using the known 21 (x,y) points, 1001 numbers of (x,y) equal-spaced cubic spline interpolated points were calculated.
// Since the test function y(x) is known exactly, the exact y-value for any x-value can be calculated.
// For each point, the interpolation error is calculated as the difference between the exact (x,y) value and the interpolated (x,y) value. 
// This program writes outputs as tab delimited results to both screen and named text file

// COMPILATION: $ g++ -lgmp -lm -std=c++11 -o cbspline-test.xx cbspline-test.cpp
// EXECUTION:   $ ./cbspline-test.xx 

// GNUPLOT 1:   gnuplot> load "cbspline-versus-exact-function-plots.gpl"
// GNUPLOT 2:   gnuplot> load "cbspline-interpolated-point-errors.gpl"

#include <iostream>
#include <fstream>
#include <boost/math/interpolators/cubic_b_spline.hpp>
using namespace std;

// ========================================================
int main(int argc, char* argv[]) {
// ========================================================

    // Vector size 21 = 20 space intervals, size 1001 = 1000 space intervals  
    std::vector<double> x_knot(21), x_exact(1001), x_cbspline(1001);
    std::vector<double> y_knot(21), y_exact(1001), y_cbspline(1001); 
    std::vector<double> y_diff(1001);

    double x;   double x0 = 0.0;    double t0 = 0.0;    
    double xstep1 = 0.5;    double xstep2 = 0.01;       

    ofstream outfile;                                   // Output data file tab delimited values
    outfile.open("cbspline-errors-1000-points.dat");    // y_diff = (y_exact - y_cbspline)

    // Handling zero-point infinity (1/x) when x = 0
    x0 = 1e-18; 
    x_knot[0] = x0; 
    y_knot[0] = (5.0/(x0))*sin(5*x0)+(x0-6);            // Selected test function

    for (int i = 1; i <= 20; i++) { // Note: Loop i begins with 1 not 0, 20 points 
        x = (i*xstep1); 
        x_knot[i] = x;
        y_knot[i] = (5.0/(x))*sin(5*x)+(x-6); 
    }

    // Execution of the cubic spline interpolation from Boost library
    // NOTE: Using xstep1 = 0.5 based on knot points stepsize (for 21 known points) 
    boost::math::cubic_b_spline<double> spline(y_knot.begin(), y_knot.end(), t0, xstep1);

    // Again handling zero-point infinity (1/x) when x = 0
    x_cbspline[0] = x_knot[0];  
    y_cbspline[0] = y_knot[0];

    for (int i = 1; i <= 1000; ++i) {       // 1000 points.
        x = (i*xstep2);  
        x_cbspline[i] = x;
        // NOTE: y = spline(x) is valid for any value of x in xrange [0:10] 
        // meaning, any x within range [x_knot.begin() and x_knot.end()]
        y_cbspline[i] = spline(x);   
    }

    // Write headers for screen display and output file
    std::cout << "# x_exact[i] \t y_exact[i] \t y_cbspline[i] \t (y_diff[i])" << std::endl;  
    outfile   << "# x_exact[i] \t y_exact[i] \t y_cbspline[i] \t (y_diff[i])" << std::endl;  

    // Again handling zero-point infinity (1/x) when x = 0
    x_exact[0] = x_knot[0]; 
    y_exact[0] = y_knot[0];
     y_diff[0] = (y_exact[0] - y_cbspline[0]);
    std::cout << x_exact[0] << "\t" << y_exact[0] << "\t" << y_cbspline[0] << "\t" << y_diff[0] << std::endl;  // Write to screen
    outfile   << x_exact[0] << "\t" << y_exact[0] << "\t" << y_cbspline[0] << "\t" << y_diff[0] << std::endl;  // Write to text file

    for (int i = 1; i <= 1000; ++i) {       // 1000 points
        x = (i*xstep2);  
        x_exact[i] = x;
        y_exact[i] = (5.0/(x))*sin(5*x)+(x-6); 
         y_diff[i] = (y_exact[i] - y_cbspline[i]);
        std::cout << x_exact[i] << "\t" << y_exact[i] << "\t" << y_cbspline[i] << "\t" << y_diff[i] << std::endl;  // Write to screen
        outfile   << x_exact[i] << "\t" << y_exact[i] << "\t" << y_cbspline[i] << "\t" << y_diff[i] << std::endl;  // Write to text file
    }
    outfile.close();

return(0);
}
// ========================================================
/*
# GNUPLOT script 1
# File: cbspline-versus-exact-function-plots.gpl

set term qt persist size 700,500
set xrange [-1:10.0]
set yrange [-12.0:20.0]
# set autoscale         # set xtics automatically
set xtics 0.5
set ytics 2.0
# set xtic auto         # set xtics automatically
# set ytic auto         # set ytics automatically
set grid x
set grid y
set xlabel "x" 
set ylabel "y(x)"
set title "Function points y(x) = (5.0/(x))*sin(5*x)+(x-6)"
set yzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
set xzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
show xzeroaxis
show yzeroaxis
show grid

plot 'cbspline-errors-1000-points.dat' using 1:2 with lines linetype 3 linewidth 1.0 linecolor 'blue' title 'exact function curve', 'cbspline-errors-1000-points.dat' using 1:3 with lines linecolor 'red' linetype 3 linewidth 1.0 title 'cubic spline interpolated curve'

*/
// ========================================================
/*
# GNUPLOT script 2
# File: cbspline-interpolated-point-errors.gpl

set term qt persist size 700,500
set xrange [-1:10.0]
set yrange [-2.50:2.5]
# set autoscale         # set xtics automatically
set xtics 0.5
set ytics 0.5
# set xtic auto         # set xtics automatically
# set ytic auto         # set ytics automatically
set grid x
set grid y
set xlabel "x" 
set ylabel "y"
set title "Function points y = (5.0/(x))*sin(5*x)+(x-6)"
set yzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
set xzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
show xzeroaxis
show yzeroaxis
show grid

plot 'cbspline-errors-1000-points.dat' using 1:4 with lines linetype 3 linewidth 1.0 linecolor 'red' title 'cubic spline interpolated errors'

*/
// ========================================================
person warazdr    schedule 19.11.2019