Voir version :
Pas de dépendances
Télécharger : 
#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; }
Voici un exemple de construction de courbe Spline. Une courbe Spline, c'est une belle courbe bien arrondie, bien douce, qui passe la ou on l'a choisi. Lancez donc le programme, vous pouvez constater qu'il y a une belle courbe rouge, qui passe par des points verts. L'idée, c'est que j'ai donné les points verts, et la spline c'est la courbe rouge. Les points verts, c'est le tableau "point" en haut du programme. Vous pouvez rajouter des points si vous voulez. Le tableau doit être pair, puisque ce sont des coordonnées x,y qui se suivent. La courbe Spline est définie entre 0 et 1. Avant de voir comment elle est construite, voyons l'utilisation : Regardez juste le main. J'appelle une fonction InitialiserSpline, je lui passe le tableau de points, le nombre de points, et un facteur 1.0 qu'on verra après. Puis je fais un for de t = 0 .. 1, en faisant des pas de 0.001. Je demande avec SplineGet le point x,y correspondant à t Et je les affiche. Très simple donc d'utilisation : P = f(t), avec f la fonction de la spline. Voila, pour l'utilisation, c'est tout ! ************************************************************* Pour les matheux, voici maintenant une explication. Entre chaque couple de points, qu'on appellera "segment", on va définir un polynôme de degré 3 : f(t) = f(t) = a + b*t + c*t^2 + d*t^3 On parlera donc de spline cubique (degré 3). On a donc nombre de segments = nombre de points -1. Pour chacun des points , on va donner un paramètre de passage. (entre 0 et 1) Ainsi, pour calculer un point de la spline de paramètre t, on choisira le segment entre deux points Pi de telle sorte que t_Pi < t < t_Pi+1 On calcule ces paramètres selon la distance chordale, c'est à dire qu'on va considérer la distance (euclidienne) entre deux points, et le paramètre choisi sera proportionnel à cette distance, par rapport à la distance globale. On ramènera tout entre 0 et 1. Ainsi, le premier point aura comme paramètre 0, le dernier 1, et les points intermédiaires des valeurs proportionnelles aux distances chordales. On aura besoin ensuite de tangentes. Nous voulons que notre spline soit C1, c'est à dire qu'à chaque point de passage, la dérivée soit continue. Note : On peut calculer des splines cubiques C2, mais le calcul est plus complexe. Un bon algorithme pour avoir de bonnes tangentes est l'algorithme de Catmul-Rom, qui, pour un point P donné, se sert du point précédent, et du suivant, pour calculer la tangente. La Tangente est parallèle au segment formé par les points suivents et précédent et s'appuie sur le paramétrage. Notez qu'un facteur peut être appliqué. Dans l'exemple, j'ai mis comme facteur 1.0 Essayez de changer, et de mettre 1.5, ou alors 2, ou d'autres valeurs. Ce facteur permettra d'arrondir davantage la courbe, ou de la raidir. Si vous mettez 0.0, la courbe devient une ligne brisée. (un polyline) L'inconvénient de Catmull Rom est qu'il ne permet pas de calculer les tangentes des premiers et derniers points. Nous devrons donc faire des cas particuliers pour les segments extrêmes. Notez que si vous considérez une spline circulaire, vous n'aurez pas ce problème. Ensuite, nous allons calculer, pour chaque segment, les facteurs a,b,c,d du polynôme local. Pour cela, dans le cas des segments milieux, nous avons 4 contraintes : les 2 points de passage, et les 2 tangentes. ça tombe bien, nous avons 4 inconnus. La fonction CalculMidSegment détaille les 4 équations, et 4 inconnus. J'ai résolu le système à part, et posé les solutions directement dans le code. C'est précalculé, c'est rapide. Pour les segments extremes, je n'ai donc pas la tangente extrême. Il me manque donc une contrainte. Je vais poser tout simplement que la dérivée seconde à l'extrémité de la courbe est nulle. ça me fait une 4e contrainte. Ainsi, je calcule tous mes polynomes pour chaque segment. Pour obtenir un point, avec la fonction SplineGet, c'est très simple : je regarde, en fonction de mon paramètre t, dans quel segment je suis (une recherche dichotomique avec FindSegment), puis j'applique simplement le polynôme local. Et j'ai donc mon point au paramètre t voulu. Notez juste que je considère, pour simplifier les calculs, que chaque polynôme local part de t = 0.0, d'ou la soustraction du paramètre à la 2e ligne de SplineGet. Notez que je pourrais aussi simplement avoir la dérivée (la vitesse), en appliquant, pour t donné, la dérivée du polynôme local. A quoi ça servirai ? Et bien tout simplement si je veux faire une voiture qui suit la courbe, la dérivée m'indiquerait tout simplement dans quelle direction est l'avant de la voiture. Et ainsi, avec une fonction de rotation d'image, je pourrais avoir ma voiture qui suit la courbe en tournant bien comme il faut... Notez également que, même si la spline est définie entre 0 et 1, on peut l'étendre : si vous donnez un nombre hors limite, la fonction FindSegment trouvera le premier ou le dernier segment, et comme un polynôme est défini sur R, vous aurez un beau prolonongement C2.
