Spline de Catmull-Rom

Principe

Dans la Spline de Hermite on a vu comment définir une spline en fournissant la liste des points et la liste des tangentes.

Cela peut-être ennuyeux de devoir fournir les tangentes. Dans certains cas, on voudrait seulement fournir la liste des points et obtenir une spline passant par tous les points.

On utilise alors l'astuce suivante :

Quand on considère le point $P_i$, on lui attribue un vecteur tangent égal à $\overrightarrow{P_{i-1}P_{i+1}}$.

À l'usage, on réalise que cette méthode donne de bons résultats mais les courbes obtenues sont un peu trop tordues. On bricole un peu et on décide de choisir comme vecteur tangent $\tau\overrightarrow{P_{i-1}P_{i+1}}$ où $\tau$ est un réglage à choisir entre $0$ et $1$. $\tau = 0,5$ est par exemple un bon choix.

$\tau$ est appelé tension.

Problème des extrémités

On est ennuyé pour $P_0$ et $P_{n-1}$.

Pour ces cas on décide de construite un point fictif par symétrie. Par exemple, pour $P_0$, on construit un point $P_{-1}$ symétrique de $P_1$ par rapport à $P_0$. Ainsi, même pour $i = 0$ on aura un $P_{i-1}$ et un $P_{i+1}$.

Formules

Considérons un arc de spline délimité par les points $P_i$ et $P_{i+1}$. On utilise les points $P_{i-1}$ et $P_{i+2}$ éventuellement construits par symétrie comme expliqué ci-dessus.

On veut construire l'arc de spline de degré 3 tel que :

$M(0) = P_i$
$M(1) = P_{i+1}$
$M'(0) = \tau\cdot(P_{i+1} - P_{i-1})$
$M'(1) = \tau\cdot(P_{i+2} - P_{i})$

Tous calculs faits on obtient $M(t) = c_0 + c_1\cdot t + c_2 \cdot t^2 + c_3 \cdot t^3$ avec

$c_0 = P_i$ ;
$c_1 = \tau(P_{i+1} - P_{i-1})$
$c_2 = -\tau(P_{i+2} + (3 - 2\tau) P_{i+1} + (\tau -3) P_{i} + 2\tau P_{i-1}$ ;
$c_3 = \tau P_{i+2} + (\tau - 2) P_{i+1} + (2 - \tau) P_{i} - \tau P_{i-1}$.

Implémentation

On souhaite faire une implémentation semblable à celle faite dans la fiche Courbe de Bézier.

Une classe Point

On reprend la classe Point identique à celle faite dans Courbe de Bézier.

La classe Point sert surtout à faciliter les calculs sur les paires (x,y), car tous les calculs que l'on fait sur x doivent être faits à l'identique sur y. On pourra utiliser la classe Point pour les vecteurs, le fonctionnement est le même.

Les données fournies

Pour produire la spline de Catmull-Rom, on devra lui fournir une liste de paires (x,y) pour les points. Rien de plus.

Exemple :

pts = [(10,50), (200, 400), (300, 300), (500, 200)]

Une classe Catmullrom

Cette classe aura la responsabilité des calculs liés à la spline de Catmull-Rom.

Ajoutez le module suivant à votre projet et complétez les fonctions.

# carmullrom.py

from point import Point

class Catmullrom:
    def __init__(self, pts, tau = 0.5):
        '''
        pts: liste de paires (x,y) représentant les points
        tau: tension de la spline. Par défaut 0.5
        '''
        assert len(pts) >= 2, "Il faut au moins deux points dans la spline."
        self.points = [Point(x,y) for x,y in pts]
        assert 0 <= tau <= 1
        self.tau = tau
    
    def get_point(self, i):
        '''
        renvoie le point d'indice i. Si i =0 ou len(pts), calcule un symétrique
        '''
        n = len(self.pts)
        assert -1 <= i <= n
        if i == -1:
            return self.pts[0] * 2 - self.pts[1]
        if i == n:
            return self.pts[n-1] * 2 - self.pts[n-2]
        return self.pts[i] 
       

    def calc_coeffs(self, i:int):
        '''
        i:indice de l'arc considéré
        la courbe de Bézier est calculée avec les points P(i-1), P(i), P(i+1), P(i+2)
        on calcule 4 coefficients :
        c0 = P(i)
        c1 = tau (P(i+1) - P(i-1))
        c2 = -tau P(i+2) + (3 - 2 tau) P(i+1) + (tau -3) P(i) + 2 tau P(i-1)
        c3 = tau P(i+2) + (tau - 2) P(i+1) +  (2 - tau) P(i) - tau P(i-1)
        La fonction renvoie c0, c1, c2, c3
        remarque: chaque coefficient est de type Point.
        '''
        # vérifie que i est bien une valeur valide : 0 <= i <= n-2
        # récupérer les P(i-1), P(i), P(i+1) et P(i+2)
        # calculer les coefficients et les renvoie
        
    def M(self, t:float, i:int):
        '''
        t: flottant entre 0 et 1, compris
        i: rang de l'arc considéré
        renvoie c0 + c1*t + c2*t**2 + c3*t**3 avec c0, c1, c2, c3
        les coeffs de la fonction précédente.
        '''
        # obtenir les valeurs des coefficients
        # renvoie le résultat du calcul

    def tangent(self, t:float, i:int):
        t: flottant entre 0 et 1, compris
        i: rang de l'arc considéré
        renvoie c1 + 2*c2*t + 3*c3*t**2 avec c0, c1, c2, c3
        '''
        # obtenir les valeurs des coefficients
        # renvoie le résultat du calcul
    
    def plot_points(self, N=100):
        '''
        renvoie la liste des points (type Point) constituant la courbe.
        N représente le nombre de points à calculer pour chaque arc.
        '''
        # préparer une liste vide
        # pour chaque arc i,
        #    pour chaque t allant de 0/N à (N-1)/N
        #        calculer M(t, i) et l'ajouter à la liste
        # compléter la liste avec une copie du dernier point de self.points
        # renvoyer la liste

Lors de l'utilisation, pour un arc donné, on calculera M(t) une bonne centaine de fois. Il serait donc bienvenu de ne pas devoir recalculer les coefficients c0, c1, c2, c3 à chaque fois. Le mieux serait de calculer les coefficients de chaque arc dès le début (__init__) et les placer dans un tableau de façon à ne plus avoir à refaire le calcul ensuite.

Utilisation avec matplotlib

On se propose d'expérimenter le module précédent avec un tracer dans matplotlib.

# democatmullrom.py
# démonstration en utilisant matplotlib

import matplotlib.pyplot as plt
from catmullrom import Catmullrom

pts = [(10,50), (200, 400), (300, 300), (500, 200)]

spline = Catmullrom(pts)

# tracé de la courbe
plot_pts = spline.plot_points()
# dans matplotlib, il faut placer les x et les y dans deux listes séparées
x = [pt.x for pt in plot_pts ]
y = [pt.y for pt in plot_pts ]
plt.plot(x, y)

# tracé des points :
x = [x for x, y in pts]
y = [y for x, y in pts]
plt.scatter(x, y, c='b')

plt.show()