Transformations affines et perspective

1. Introduction

Ce document présente une implémentation Python des transformations affines et de la projection perspective.

Un point de l'espace ou un vecteur est représenté par ses coordonnées homogènes :

$$X=\left(\begin{array}{c}\hfill x\hfill \\ \hfill y\hfill \\ \hfill z\hfill \\ \hfill w\hfill \end{array}\right)\text{(1)}$$

Pour un point w=1, pour un vecteur w=0.

Effectuer une transformation affine (ou une transformation linéaire pour un vecteur) consiste à multiplier cette matrice colonne par une matrice 4x4.

Pour obtenir une projection perspective, une matrice de perspective doit être appliquée, qui modifie la valeur de w. Enfin la division perspective est effectuée : les coordonnées (x,y,z) sont divisées par w pour obtenir les coordonnées du point image dans le plan image.

2. Classe Python

Le programme comporte une classe Transformation3D, dont les principales données sont la matrice de transformation et son inverse. Une pile de transformation est aussi utilisée, pour permettre de sauvegarder une transformation. Les matrices sont des tableaux numpy.

import numpy
import math
            

Voici le constructeur. Les matrices sont initialisées à l'identité.

class Transformation3D:
    def __init__(self):
        self.M = numpy.identity(4,dtype=numpy.float64)
        self.invM = numpy.identity(4,dtype=numpy.float64)
        self.pile = []
        self.deg2rad = math.pi/180.0
            

Les deux fonctions suivantes permettent d'ajouter une transformation à la pile et de récupérer la dernière transformation :

    def save(self):
        self.pile.append((self.M,self.invM))
    def restore(self):
        (self.M,self.invM) = self.pile.pop()
            

La fonction suivante effectue une composition avec une autre transformation. La matrice M est multipliée à droite; son inverse est multipliée à gauche.

    def compose(self,transform):
        self.M = self.M.dot(transform.M)
        self.invM = transform.invM.dot(self.invM)
            

Voici les transformations affines élémentaires. Chaque fonction effectue une composition avec une transformation affine.

    def translate(self,b1,b2,b3):
        self.M = self.M.dot(numpy.array([[1.0,0,0,b1],[0,1.0,0,b2],[0,0,1.0,b3],[0,0,0,1.0]]))
        self.invM = numpy.array([[1.0,0,0,-b1],[0,1.0,0,-b2],[0,0,1.0,-b3],[0,0,0,1.0]]).dot(self.invM)
    def rotateZrad(self,a):
        c = math.cos(a)
        s = math.sin(a)
        self.M = self.M.dot(numpy.array([[c,-s,0,0],[s,c,0,0],[0,0,1.0,0],[0,0,0,1.0]]))
        self.invM = numpy.array([[c,s,0,0],[-s,c,0,0],[0,0,1.0,0],[0,0,0,1.0]]).dot(self.invM)
    def rotateZ(self,a):
        self.rotateZrad(a*self.deg2rad)
    def rotateXrad(self,a):
        c = math.cos(a)
        s = math.sin(a)
        self.M = self.M.dot(numpy.array([[1.0,0,0,0],[0,c,-s,0],[0,s,c,0],[0,0,0,1.0]]))
        self.invM = numpy.array([[1.0,0,0,0],[0,c,s,0],[0,-s,c,0],[0,0,0,1.0]]).dot(self.invM)
    def rotateX(self,a):
        self.rotateXrad(a*self.deg2rad)
    def rotateYrad(self,a):
        c = math.cos(a)
        s = math.sin(a)
        self.M = self.M.dot(numpy.array([[c,0,s,0],[0,1.0,0,0],[-s,0,c,0],[0,0,0,1.0]]))
        self.invM = numpy.array([[c,0,-s,0],[0,1.0,0,0],[s,0,c,0],[0,0,0,1.0]]).dot(self.invM)
    def rotateY(self,a):
        self.rotateYrad(a*self.deg2rad)
    def scale(self,sx,sy,sz):
        self.M = self.M.dot(numpy.array([[sx,0,0,0],[0,sy,0,0],[0,0,sz,0],[0,0,0,1.0]]))
        self.invM = numpy.array([[1.0/sx,0,0,0],[0,1.0/sy,0,0],[0,0,1.0/sz,0],[0,0,0,1.0]]).dot(self.invM)

            

La fonction suivante compose avec la matrice de perspective. Voir Projection perspective pour la définition de cette matrice. d est la distance algébrique (non nulle) entre le centre de projection et le plan image.

    def perspective(self,d,a=0,b=1):
        self.M = self.M.dot(numpy.array([[d,0,0,0],[0,d,0,0],[0,0,a,b],[0,0,1,0]]))
        self.invM = numpy.array([[1.0/d,0,0,0],[0,1.0/d,0,0],[0,0,0,1],[0,0,1.0/b,-a/b]]).dot(self.invM)
            

La fonction suivante effectue la multiplication d'une matrice colonne contenant des coordonnées homogènes (sans division perspective).

    def transform(self,X):
        return self.M.dot(X) 
            

La fonction suivante effectue la transformation d'un tableau. Chaque ligne comporte les 4 coordonnées à transformer.

    def transformArray(self,array):
        n = array.shape[0]    
        a = numpy.zeros((n,4),dtype=numpy.float64)
        for k in range(n):
            a[k] = self.M.dot(array[k].T).T
        return a
            

Les fonctions suivantes transforment un point, un vecteur et un vecteur normal. Seule la multiplication matricielle est appliquée (pas de division perspective). Les coordonnées du point (ou les composantes du vecteur) sont fournies sous forme d'une liste de 3 valeurs ou d'un tableau numpy à une dimension. Le résultat est sous forme d'une matrice ligne à 4 éléments, avec la coordonnées w.

    def transformPoint(self,p):
        return self.M.dot(numpy.array([[p[0]],[p[1]],[p[2]],[1.0]])).T
    def transformVect(self,v):
        return self.M.dot(numpy.array([[v[0]],[v[1]],[v[2]],[0.0]])).T
    def transformNormal(self,n):
        return self.invM.T.dot(numpy.array([[n[0]],[n[1]],[n[2]],[0.0]])).T
            

Les fonctions suivantes effectuent la transformation d'un tableau de points ou de vecteurs.

    def transformPointArray(self,array):
        n = array.shape[0]    
        a = numpy.zeros((n,4),dtype=numpy.float64)
        for k in range(n):
            a[k] = self.transformPoint(array[k])
        return a
    def transformVectArray(self,array):
        n = array.shape[0]    
        a = numpy.zeros((n,4),dtype=numpy.float64)
        for k in range(n):
            a[k] = self.transformVect(array[k])
        return a
    def transformNormalArray(self,array):
        n = array.shape[0]    
        a = numpy.zeros((n,4),dtype=numpy.float64)
        for k in range(n):
            a[k] = self.transformNormal(array[k])
        return a

            

La fonction suivante effectue la division perspective pour un point :

    def division(self,X):
        return numpy.array([X[0]/X[3],X[1]/X[3],X[2]/X[3]])
            

La même chose pour un tableau de points :

    def divisionArray(self,array):
        n = array.shape[0]    
        a = numpy.zeros((n,3),dtype=numpy.float64)
        for k in range(n):
            a[k] = self.division(array[k])
        return a
            

Exemple d'utilisation sans perspective (projection orthoscopique) :

from Transformation3D import *
import numpy
from matplotlib.pyplot import *
t = Transformation3D()
c = numpy.array([[-1,-1,0],[1,-1,0],[1,1,0],[-1,1,0],[-1,-1,0]],dtype=numpy.float64)
figure(figsize=(6,6))
for n in range(18):
    t.rotateZ(10.0)
    t.save()
    t.translate(2.0,0,0)
    cp = t.transformPointArray(c)
    t.restore()
    plot(cp.T[0],cp.T[1])
axis([-5,5,-5,5])
grid()

                

Figure pleine page

3. Application : effet de perspective

L'objet considéré est une grille carrée dont les mailles sont carrées. Elle est définie dans le plan Oxy, centrée sur l'origine. Les tranformations affines et la projection perspective transforment une droite en droite, éventuellement réduite à un point. Il suffit donc de calculer les extrémités des segments définissant la grille. Les coordonnées sont stockées dans une matrice dont chaque ligne contient les coordonnées homogènes du point (avec w=1).

L = 1.0 # demi-longueur d'une ligne
n = 5
nt = (2*n+1) # nombre de lignes dans chaque dimension
grille = numpy.zeros((nt*nt*2,4),dtype=numpy.float64)
dL = 2*L/(nt-1)
k = 0
for i in range(-n,n+1):
    grille[k] = numpy.array([i*dL,-L,0.0,1.0])
    k += 1
    grille[k] = numpy.array([i*dL,L,0.0,1.0])
    k += 1
for j in range(-n,n+1):
    grille[k] = numpy.array([-L,j*dL,0.0,1.0])
    k += 1
    grille[k] = numpy.array([L,j*dL,0.0,1.0])
    k+= 1
            

La fonction suivante trace la grille :

def tracerGrille(grille):
    for k in range(nt*2):
        i = 2*k
        plot([grille[i][0],grille[i+1][0]],[grille[i][1],grille[i+1][1]],color="b")
            

figure(figsize=(6,6))
axis([-2,2,-2,2])
tracerGrille(grille)
            

Figure pleine page

Comme premier exemple de transformation, on effectue une rotation autour de l'axe Oz :

trans = Transformation3D()
trans.rotateZ(45)
grille2 = trans.transformArray(grille)
figure(figsize=(6,6))
axis([-2,2,-2,2])
tracerGrille(grille2)
            

Figure pleine page

Faisons à présent une rotation de -45 degrés autour de l'axe Ox suivie d'une translation pour amener le centre de la grille en z=-10. Pour faire ces transformations affines, il faut se rappeler que la dernière transformation définie s'applique la première.

trans = Transformation3D()
trans.translate(0,0,-10)
trans.rotateX(-45)
grille2 = trans.transformArray(grille)
figure(figsize=(6,6))
axis([-2,2,-2,2])
tracerGrille(grille2)
            

Figure pleine page

Il s'agit là d'une projection orthoscopique : les lignes de la grille restent parallèles sur l'image. Pour obtenir une perspective, on définit tout d'abord une matrice de perspective, que l'on compose avec la transformation affine précédente. La distance entre le centre de projection et le plan image est d=5.

Figure pleine page

projection = Transformation3D()
projection.perspective(-5.0)
projection.compose(trans)
             

La transformation des points doit inclure la division perspective :

grille2 = projection.transformArray(grille)
grille3 = projection.divisionArray(grille2)
figure(figsize=(6,6))
axis([-2,2,-2,2])
tracerGrille(grille3)
             

Figure pleine page

Voyons l'effet d'une réduction de d (sans changer la distance de l'objet) :

projection = Transformation3D()
projection.perspective(-2.0)
projection.compose(trans)
grille2 = projection.transformArray(grille)
grille3 = projection.divisionArray(grille2)
figure(figsize=(6,6))
axis([-2,2,-2,2])
tracerGrille(grille3)
             

Figure pleine page

L'image de la grille est plus petite : la réduction de d revient à utiliser un objectif de focale plus petite, dont le champ est plus large. Pour agrandir l'image, il faut rapprocher la caméra de la grille, ce qui revient à rapprocher la grille :

trans = Transformation3D()
trans.translate(0,0,-3)
trans.rotateX(-45)
projection = Transformation3D()
projection.perspective(-2.0)
projection.compose(trans)
grille2 = projection.transformArray(grille)
grille3 = projection.divisionArray(grille2)
figure(figsize=(6,6))
axis([-2,2,-2,2])
tracerGrille(grille3)
             

Figure pleine page

L'effet de perspective est plus marquée, car l'objet est plus proche de la caméra.

Pour s'approcher d'une projection orthoscopique, il faut au contraire éloigner l'objet. Comme cet éloignement conduit à une réduction de son image, il faut aussi augmenter la distance d (augmenter la focale de l'objectif photo).

trans = Transformation3D()
trans.translate(0,0,-50)
trans.rotateX(-45)
projection = Transformation3D()
projection.perspective(-30.0)
projection.compose(trans)
grille2 = projection.transformArray(grille)
grille3 = projection.divisionArray(grille2)
figure(figsize=(6,6))
axis([-2,2,-2,2])
tracerGrille(grille3)
             

Figure pleine page

Pour réduire l'effet de perspective dans une prise de vue photographique, il faut donc s'éloigner de l'objet, ce qui nécessite d'augmenter la focale de l'objectif pour que l'objet soit assez grand sur l'image. C'est la distance de l'objet qui déterminer l'amplitude de l'effet de perspective, pas la focale de l'objectif. Cette propriété se comprend aisément si l'on remarque que l'effet de perspective est d'autant plus grand que les rayons sont inclinés par rapport à l'axe. Lorsqu'on éloigne l'objet, on réduit l'inclinaison des rayons; la distance d n'a aucun effet sur cette inclinaison.