hulp gevraagd bij een cubic spline berekening

[michelw]

Beste wiskunde guru's,

Ik ben nieuw op dit forum en heb een vraag hoe je een cubic spline berekening opzet van begin tot eind.
De bedoeling is om de functie op te zetten in C++ voor mijn cadcam programma.
Heb al verschillende video's op youtube bekeken en voorbeelden gezien, maar ik kom er niet uit. Mischien weet iemand
een begrijpelijk voorbeeld.

[siep]

Stap 1 en 2 in jouw schema leveren de constanten a, b, c en d voor de derdegraadsfunctie
y=ax^3+bx^2+cx+d
waarvan de grafiek precies door de 4 gegeven punten P1 t/m P4 loopt (zie plaatje hierboven).
Hier is y een functie van x, dus elke x levert precies één waarde voor y.
Dit is niet wat je zoekt.

De 3 splines in jouw plaatje worden ieder gevormd door 2 derdegraadsfuncties per spline:
x=d_xt^3+c_xt^2+b_xt+a_x
y=d_yt^3+c_yt^2+b_yt+a_y
dus x als functie van t en y als functie van diezelfde t per spline, waarbij steeds t \in [0, 1]
(constanten a t/m d zijn hier veranderd in d t/m a voor de aansluiting met onderstaande web-pagina's).

De berekening van de spline constanten (a t/m d) vind je hier:
http://mathworld.wolfram.com/CubicSpline.html
Vergelijking (18) op die pagina (met de tridiagonale matrix) kan je oplossen via dit algoritme:
https://en.wikipedia.org/wiki/Tridiagon ... _algorithm.

Als ik dit omzet in C, is dit de kern die je kan gebruiken:
(de robuustheid, een flexibele input/output, etc zal je waarschijnlijk allemaal nog willen veranderen)

Code: Selecteer alles

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
// CUBIC SPLINES
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------

// Cubic spline solution: see Wolfram MathWorld:
//     http://mathworld.wolfram.com/CubicSpline.html
// Uses solution for tridiagonal matrix, see: Wikipedia:
//     https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

// INPUT:
const int npoints=4; // 4 points
// matrix C [point] [dimension: 0=x, 1=y]
double C[npoints][2]={{0,0}, {10,10}, {20,9}, {15,0}};
// OUTPUT:
// solution matrix S:
// S [spline number: 0..(npoints-2)] [2 dimensions: x,y] [4 constants: a,b,c,d]
double S[npoints-1][2][4];

// calculate all cubic spline values:
void calcCubicSplines(void)
{
double w, b[npoints], d[npoints], x[npoints];
int i, dim;
int n = npoints-1; // number of splines

for(dim=0;dim<2;dim++){
  // define d[]:
  d[0]=3.0*(C[1][dim]-C[0][dim]);
  for(i=1;i<n;i++)
	d[i]=3.0*(C[i+1][dim]-C[i-1][dim]);
  d[n]=3.0*(C[n][dim]-C[n-1][dim]);
  // forward sweep: (simplified tridiagonal solution: all a[i]=1 and all c[i]=1)
  b[0]=2.0;
  for(i=1;i<npoints;i++){
	w = 1.0/b[i-1];
	b[i] = 4.0 - w;
	d[i] -= (w*d[i-1]);
	}
  b[n] -= 2.0;
  // calculate solution vector x[i] = D[i]:
  // (Wikipedia x[] = Wolfram D[])
  x[n]=d[n]/b[n];
  for(i=n-1;i>=0;i--)
	x[i]=(d[i]-x[i+1])/b[i];
  // calculate spline S[i][dim] a, b, c and d:
  for(i=0;i<n;i++){
	S[i][dim][0]=C[i][dim];
	S[i][dim][1]=x[i];
	S[i][dim][2]=3.0*(C[i+1][dim]-C[i][dim]) - 2.0*x[i] - x[i+1];
	S[i][dim][3]=2.0*(C[i][dim]-C[i+1][dim]) + x[i] + x[i+1];
	}
  }

}

// get x and y coordinates for spline i (i = 0..npoints-2) where t in interval [0, 1]:
void getSplinePoint(int i, double t, double *x, double *y)
{
(*x) = ((S[i][0][3]*t+S[i][0][2])*t+S[i][0][1])*t+S[i][0][0];
(*y) = ((S[i][1][3]*t+S[i][1][2])*t+S[i][1][1])*t+S[i][1][0];
}

Voor de 4 punten P1 t/m P4 uit jouw voorbeeld krijg ik hiermee het volgende resultaat:



Punt Q op spline 2 (als we bij nul beginnen met tellen) met t=0.5 levert via getSplinePoint(2, 0.5, &x, &y):
x = 18.2345...
y = 4.1404...

Kom je hiermee verder?

[michelw]

Beste Arie,

Dit is een fantastisch antwoord !!

Ik heb uw c code draaiend binnen 5 minuten in c++

Een minimalistische aanpassing naar c++ gemaakt in het programma QT. Sinds kort ben ik overgestapt van QT naar CodeBlocks. Wil je QT professioneel gebruiken kost het ca 4000-5000 euro per jaar. CodeBlocks is gratis en heeft bijna dezelfde functionaliteiten alleen is er minder debug (minder uitgebreide foutmeldingen rapporten).
Dus bij moelijke zaken zet ik vaak de code op in QT, als deze dan getest is, dan hevel ik de code over naar CodeBlocks.
Zelf ben ik ca 2 jaar geleden begonnen met low level c, sinds een jaar overgestapt naar c++.

Code: Selecteer alles

#include "mainwindow.h"
#include "ui_mainwindow.h"
#include <iostream>

//this can be deactived when activated in mainwindow header file..
void calcCubicSplines();
void getSplinePoint(int i, double t/*, double *x, double *y*/);

// INPUT:
const int npoints=4; // 4 points
// matrix C [point] [dimension: 0=x, 1=y]
double C[npoints][2]={{0,0}, {10,10}, {20,9}, {15,0}};
// OUTPUT:
// solution matrix S:
// S [spline number: 0..(npoints-2)] [2 dimensions: x,y] [4 constants: a,b,c,d]
double S[npoints-1][2][4];

MainWindow::MainWindow(QWidget *parent) :
    QMainWindow(parent),
    ui(new Ui::MainWindow)
{
    ui->setupUi(this);
}

MainWindow::~MainWindow()
{
    delete ui;
}

void MainWindow::on_pushButton_clicked()
{
    calcCubicSplines();
    getSplinePoint(2,0.5);
}

// calculate all cubic spline values:
void calcCubicSplines(){
    double w, b[npoints], d[npoints], x[npoints];
    int i, dim;
    int n = npoints-1; // number of splines

    for(dim=0;dim<2;dim++){
        // define d[]:
        d[0]=3.0*(C[1][dim]-C[0][dim]);
        for(i=1;i<n;i++)
            d[i]=3.0*(C[i+1][dim]-C[i-1][dim]);
        d[n]=3.0*(C[n][dim]-C[n-1][dim]);
        // forward sweep: (simplified tridiagonal solution: all a[i]=1 and all c[i]=1)
        b[0]=2.0;
        for(i=1;i<npoints;i++){
            w = 1.0/b[i-1];
            b[i] = 4.0 - w;
            d[i] -= (w*d[i-1]);
        }
        b[n] -= 2.0;
        // calculate solution vector x[i] = D[i]:
        // (Wikipedia x[] = Wolfram D[])
        x[n]=d[n]/b[n];
        for(i=n-1;i>=0;i--)
            x[i]=(d[i]-x[i+1])/b[i];
        // calculate spline S[i][dim] a, b, c and d:
        for(i=0;i<n;i++){
            S[i][dim][0]=C[i][dim];
            S[i][dim][1]=x[i];
            S[i][dim][2]=3.0*(C[i+1][dim]-C[i][dim]) - 2.0*x[i] - x[i+1];
            S[i][dim][3]=2.0*(C[i][dim]-C[i+1][dim]) + x[i] + x[i+1];
        }
    }

}

// get x and y coordinates for spline i (i = 0..npoints-2) where t in interval [0, 1]:
void getSplinePoint(int i, double t/*, double *x, double *y*/){
    /*(*x)*/ double x = ((S[i][0][3]*t+S[i][0][2])*t+S[i][0][1])*t+S[i][0][0];
    /*(*y)*/ double y = ((S[i][1][3]*t+S[i][1][2])*t+S[i][1][1])*t+S[i][1][0];
    std::cout<< "x:" << x << " y:" << y << std::endl;
}

Dit is het cad cam programma waar ik over sprak. Dit is gebouwd in Codeblocks. Het programma heeft een opengl interface.


Mijn eerste cadcam programma is geschreven in QT, had een cubic bezier spline, deze werkte dus niet goed voor
de start en eindpunten. Het programma genereerd offset contours en g_code met 1 click. Enkel de spline functie was niet correct en dat geeft gelijk problemen bijvoorbeeld voor het inladen van dxf files.
https://imgshare.io/images/2020/02/10/qt.png

Nu kan ik uw spline code implementeren (in de opengl loop zetten, koppelen aan de data array.. etc)
en verder uit testen. Dit was het laatste vraagstuk op tekengebied.
Nu is het overblijvende vraagstuk een nesting algoritme, maar dat heeft geen haast !!

Heeft u de grafiek gemaakt in matlab?

U heeft mij echt heel erg geholpen !! waarvoor veel dank. Ik kan u zeggen dat na meer dan een week proberen ik nog steeds niet verder kwam. Maar nu is het in 1 keer opgelost !!

[michelw]

Beste Arie,

Bijgaand een array setup die ik zelf vaak gebruik. De array memory is automatisch allocated door de vector<vector<double>> C;. Zodoende kan de array zijn eigen groote bepalen (memory allocation) tijdens het invoegen van de data.
(ter ondersteuning voor de lezers een 1d vector = 1d array, de naam vector is verwarrend en geen wiskunde vector oid).

Op deze manier kun je heel eenvoudig spline punten in de array toevoegen in de opengl loop.
Bijgaand de code die ik dan vandaag heb getest met c++:

Code: Selecteer alles

	
    //using namespace std;
    //vector<vector<double>> C; //2d array
    C.push_back({0,0});
    C.push_back({10,10});
    C.push_back({20,9}); //cout<< P[2].at(0)  <<  P[2].at(1) << endl;
    C.push_back({15,0});

    double S[C.size()-1][2][4];
    double w, b[C.size()], d[C.size()], x[C.size()];
    int i, dim;
    int n = C.size()-1; // number of splines

    for(dim=0;dim<2;dim++){
        // define d[]:
        d[0]=3.0*(C[1][dim]-C[0][dim]);
        for(i=1;i<n;i++)
            d[i]=3.0*(C[i+1][dim]-C[i-1][dim]);
        d[n]=3.0*(C[n][dim]-C[n-1][dim]);
        // forward sweep: (simplified tridiagonal solution: all a[i]=1 and all c[i]=1)
        b[0]=2.0;
        for(i=1;i<C.size();i++){
            w = 1.0/b[i-1];
            b[i] = 4.0 - w;
            d[i] -= (w*d[i-1]);
        }
        b[n] -= 2.0;
        // calculate solution vector x[i] = D[i]:
        // (Wikipedia x[] = Wolfram D[])
        x[n]=d[n]/b[n];
        for(i=n-1;i>=0;i--)
            x[i]=(d[i]-x[i+1])/b[i];
        // calculate spline S[i][dim] a, b, c and d:
        for(i=0;i<n;i++){
            S[i][dim][0]=C[i][dim];
            S[i][dim][1]=x[i];
            S[i][dim][2]=3.0*(C[i+1][dim]-C[i][dim]) - 2.0*x[i] - x[i+1];
            S[i][dim][3]=2.0*(C[i][dim]-C[i+1][dim]) + x[i] + x[i+1];
        }
    }

    glColor3f(1,1,1);
    glPointSize(5);
    glBegin(GL_POINTS);
    for(int p=0; p<=C.size()-1; p++){ //spline points
            glVertex2d(C[p][0],C[p][1]);

    }
    glEnd();

    glLineWidth(1);
    glBegin(GL_LINE_STRIP);
    for(int p=0; p<=C.size()-2; p++){ //spline points
        for(double t=0; t<1; t+=0.01){ //time 0-1 for each spline
            double px = ((S[p][0][3]*t+S[p][0][2])*t+S[p][0][1])*t+S[p][0][0];
            double py = ((S[p][1][3]*t+S[p][1][2])*t+S[p][1][1])*t+S[p][1][0];
            //std::cout<< "x:" << px << " y:" << py << std::endl;
            glVertex2d(px,py);
        }
    }
    glEnd();

    C.clear(); //empty 2d array

En het plaatje met de uitkomst in het codeblocks cadcam programma :