Spline curves

Spline curves

See version :

Pas de dépendances

Download :

#include <math.h>
#include <sdl/sdl.h>
#include <stdio.h>
#include <stdlib.h>

#ifdef _WIN32
#pragma comment(lib,"sdl.lib")
#pragma comment(lib,"sdlmain.lib")
#endif

#define XRES 400
#define YRES 300

float points[] = {10,10,100,70,200,50,150,250,380,200,300,10};

/* ------------------ Partie Spline --------------------------- */

typedef struct
{
    float a[2],b[2],c[2],d[2];
    float t;
} SplineSegment;


typedef struct
{
    SplineSegment* segs;
    int nb;
} Spline;

float* AutoParamDistChordal(float* tabpoints,int nb)
{
    float* res;
    int i;
    float max;
    res = malloc(nb*sizeof(float));
    res[0] = 0.0f;
    for(i=1;i<nb;i++)
    {
        float vx,vy;
        vx = tabpoints[2*i] - tabpoints[2*(i-1)];
        vy = tabpoints[2*i+1] - tabpoints[2*(i-1)+1];
        res[i] = res[i-1] + (float)sqrt(vx*vx+vy*vy);
    }
    max = res[nb-1];
    for(i=0;i<nb;i++)
        res[i]/=max;
    return res;
}

float* CatmullRom(float* tabpoints,float* t,int nb,float facteur)
{
    int i;
    float* res = malloc(2*nb*sizeof(float));
    /* extremités impossibles */
    for(i=1;i<nb-1;i++)
    {
        float dt = t[i+1] - t[i-1];
        res[2*i] = facteur * (tabpoints[2*(i+1)] - tabpoints[2*(i-1)])/dt;
        res[2*i+1] = facteur * (tabpoints[2*(i+1)+1] - tabpoints[2*(i-1)+1])/dt;
    }
    return res;
}

void CalculFirstSegment(SplineSegment* seg,float* P,float* T,float dt)
{
    /*
        forme a trouver : f(t) = a + b*t + c*t^2 + d*t^3
        4 inconnus.
        4 equations : 
        f(0) = P0
        f(dt) = P1
        f''(0) = 0
        f'(dt) = T1
        --> résolution : 
        a = P0, c = 0, d = 1/2/dt^3*(P0+dt*T1-P1), b = -1/2*(3*P0+dt*T1-3*P1)/dt
    */
    int i;
    for(i=0;i<2;i++)
    {
        float P0,P1,T1;
        P0 = P[i];
        P1 = P[2+i];
        T1 = T[2+i];
        seg->a[i] = P0;
        seg->b[i] = (float)(-1./2*(3*P0+dt*T1-3*P1)/dt);
        seg->c[i] = 0;
        seg->d[i] = (float)(1./2./(dt*dt*dt)*(P0+dt*T1-P1));
    }
}

void CalculMidSegment(SplineSegment* seg,float* P,float* T,float dt)
{
    /*
        forme a trouver : f(t) = a + b*t + c*t^2 + d*t^3
        4 inconnus.
        4 equations : 
        f(0) = P0
        f(dt) = P1
        f'(0) = T0
        f'(dt) = T1
        --> résolution : 
        b = T0, a = P0, d = 1/dt^3*(2*P0+T0*dt-2*P1+dt*T1), c = -(2*T0*dt+3*P0-3*P1+dt*T1)/dt^2
    */
    int i;
    for(i=0;i<2;i++)
    {
        float P0,P1,T0,T1;
        P0 = P[i];
        P1 = P[2+i];
        T0 = T[i];
        T1 = T[2+i];
        seg->a[i] = P0;
        seg->b[i] = T0;
        seg->c[i] = -(2*T0*dt+3*P0-3*P1+dt*T1)/(dt*dt);
        seg->d[i] = (2*P0+T0*dt-2*P1+dt*T1)/(dt*dt*dt);
    }
}

void CalculLastSegment(SplineSegment* seg,float* P,float* T,float dt)
{
    /*
        forme a trouver : f(t) = a + b*t + c*t^2 + d*t^3
        4 inconnus.
        4 equations : 
        f(0) = P0
        f(dt) = P1
        f'(0) = T0
        f''(dt) = 0
        --> résolution : 
        b = T0, a = P0, d = 1/2/dt^3*(P0+T0*dt-P1), c = -3/2/dt^2*(P0+T0*dt-P1)
    */
    int i;
    for(i=0;i<2;i++)
    {
        float P0,P1,T0;
        P0 = P[i];
        P1 = P[2+i];
        T0 = T[i];
        seg->a[i] = P0;
        seg->b[i] = T0;
        seg->c[i] = (float)(-3./2./(dt*dt)*(P0+T0*dt-P1));
        seg->d[i] = (float)(1./2./(dt*dt*dt)*(P0+T0*dt-P1));
    }
}

Spline* InitialiserSpline(float* tabpoints,int nb,float catmulromfacteur)
{
    float* params;
    float* tang;
    Spline* spl;
    int i;
    spl = malloc(sizeof(Spline));
    spl->nb = nb-1;
    spl->segs = malloc(spl->nb*sizeof(SplineSegment));
    params = AutoParamDistChordal(tabpoints,nb);
    for(i=0;i<nb-1;i++)
        spl->segs[i].t = params[i];
    tang = CatmullRom(tabpoints,params,nb,catmulromfacteur);
    CalculFirstSegment(&spl->segs[0],tabpoints,tang,params[1]-params[0]);
    for(i=1;i<nb-2;i++)
        CalculMidSegment(&spl->segs[i],tabpoints+2*i,tang+2*i,params[i+1]-params[i]);
    CalculLastSegment(&spl->segs[nb-2],tabpoints+2*(nb-2),tang+2*(nb-2),params[nb-1]-params[nb-2]);
    free(tang);
    free(params);
    return spl;
}

void FreeSpline(Spline* spl)
{
    free(spl);
}

SplineSegment* FindSegment(Spline* spl,float t)
{
    int min,max,med;
    min = 0;
    max = spl->nb;
    while(max-min>1)
    {
        med = (max+min)/2;
        if (t>spl->segs[med].t)
            min = med;
        else
            max = med;
    }
    return &spl->segs[min];
}


void SplineGet(Spline* spl,float t,float* outX,float* outY)
{
    SplineSegment* seg = FindSegment(spl,t);
    t-= seg->t;
    *outX = seg->a[0] + seg->b[0]*t + seg->c[0]*t*t + seg->d[0]*t*t*t;
    *outY = seg->a[1] + seg->b[1]*t + seg->c[1]*t*t + seg->d[1]*t*t*t;
}



/* ----------------- Fin Partie Spline ------------------------------------- */

void SDL_PutPixel32(SDL_Surface *surface, int x, int y, Uint32 pixel)
{
    Uint8 *p;
    if (x<0 || y<0 || x>surface->w-1 || y>surface->h-1)
        return;
    p = (Uint8*)surface->pixels + y * surface->pitch + x * 4;
    *(Uint32*)p = pixel;
}

int KeyHit()
{
    SDL_Event e;
    if (SDL_PollEvent(&e))
        if (e.type == SDL_KEYDOWN)
            return 1;
    return 0;
}

int main(int argc,char** argv)
{
    SDL_Surface* screen;
    Spline* spl;
    float t;
    int i;
    SDL_Init(SDL_INIT_VIDEO);
    screen=SDL_SetVideoMode(XRES,YRES,32,SDL_SWSURFACE|SDL_DOUBLEBUF);  
    spl = InitialiserSpline(points,sizeof(points)/(2*sizeof(int)),1.0);
    /* dessine la spline a 0.001 de pas */
    for(t=0.0f;t<=1.0f;t+=0.001f)
    {
        float x,y;
        SplineGet(spl,t,&x,&y);
        SDL_PutPixel32(screen,(int)x,(int)y,0xFF0000);
    }
    /* dessine les points de passage */
    for(i=0;i<sizeof(points)/(2*sizeof(int));i++)
        SDL_PutPixel32(screen,(int)points[2*i],(int)points[2*i+1],0x00FF00);
    SDL_Flip(screen);
    while(!KeyHit())
        SDL_Delay(10);
    FreeSpline(spl);
    SDL_Quit();
    return 0;
}



Explanations

	No explanations yet.