Déterminer facilement si un point appartient à un triangle

Deux solutions très simples pour déterminer si un point appartient à un triangle

Publié le 8 avril 2017

On voit parfois sur des forums des questions à propos de méthodes pour déterminer l'appartenance d'un point à un polygone. Dans ce billet, je vais donner deux solutions très faciles qui marchent pour un triangle. Nous verrons dans un autre billet comment généraliser ces méthodes à un polygone convexe quelconque.

Nous allons voir tout d'abord d'un point de vue mathématique comment faire pour déterminer l'appartenance ou non d'un point à un triangle, puis l'application dans un langage de programmation. Bien entendu, nous utiliserons le C++. Ainsi vous pouvez directement aller voir l'application informatique sans regarder la partie mathématiques si l'envie vous en prend.

Des mathématiques : les coordonnées barycentriques

Je ne vais pas donner un cours sur ce que sont les coordonnées barycentriques mais juste montrer comment cette notion va nous permettre de déterminer si un point appartient à un triangle.

En considérant un triangle \( ABC \) non plat, on peut dire que \( P \), un point quelconque du plan appartenant au triangle \( ABC \), a pour coordonnées les trois déterminants suivant : $$\alpha = \left| \overrightarrow{PB},\overrightarrow{PC} \right|,\ \beta=\left|\overrightarrow{PC},\overrightarrow{PA}\right|,\ \gamma=\left| \overrightarrow{PA},\overrightarrow{PB}\right| .$$ Ces trois déterminants sont relatifs à une base quelconque du plan, cela veut dire qu'ils sont vrais quel que soit le repère du plan dans lequel ils sont exprimés. L'expression de ces trois déterminants implique l'égalité suivante : $$\alpha \cdot \overrightarrow{PA} + \beta \cdot \overrightarrow{PB} + \gamma\cdot\overrightarrow{PC} =0$$ avec : $$\alpha+\beta+\gamma = \left| \overrightarrow{PB},\overrightarrow{PC} \right| + \left|\overrightarrow{PC},\overrightarrow{PA} \right| + \left| \overrightarrow{PA},\overrightarrow{PB} \right| = \left| \overrightarrow{AB},\overrightarrow{AC} \right| \neq 0 .$$

De là s'offrent à nous deux choix : calculer des aires de triangle ou résoudre un système d'équations linéaires.

Calculer des aires de triangle

Si la base choisie est \( \left( \overrightarrow{\scriptstyle AB}, \overrightarrow{\scriptstyle AC} \right) \), les coordonnées barycentriques de \( P \) sont alors \( \left( \alpha , \beta , \gamma \right) \) vérifiant \( \alpha + \beta + \gamma = 1 \).

Ces coordonnées s'interprètent comme des rapports d'aires algébriques de la façon suivante : $$\begin{cases} \alpha = \frac{\textrm{Aire}(PBC)}{\textrm{Aire}(ABC)} \\ \beta = \frac{\textrm{Aire}(APC)}{\textrm{Aire}(ABC)} \\ \gamma = \frac{\textrm{Aire}(ABP)}{\textrm{Aire}(ABC)} \end{cases}$$

La question est alors, comment calculer ces aires ? Et bien c'est très simple.

L'aire du parallélogramme défini par deux vecteurs \( \overrightarrow{u} \) et \( \overrightarrow{v} \) est la norme de leur produit vectoriel : \[ S_{p} = \left\| \vec{u} \wedge \vec{v} \right\| .\] On peut calculer l'aire d'un triangle à partir de cette formule : \[ S = \frac{1}{2} \left\| \overrightarrow{AB} \wedge \overrightarrow{AC} \right\| .\] Un repère orthonormé étant donné, l'aire du triangle \( ABC \) peut être calculée à partir des coordonnées des sommets.

Dans le plan, si les coordonnées de \( A \), \( B \) et \( C \) sont données par \(A\left(x_{A};y_{A}\right)\), \(B\left(x_{B};y_{B}\right)\) et \(C\left(x_{C};y_{C}\right)\), alors l'aire \(S\) est la moitié de la valeur absolue du déterminant : \[ S = \frac{1}{2} \left| \det \begin{pmatrix} x_B-x_A & x_C-x_A \\ y_B-y_A & y_C-y_A \end{pmatrix} \right| = \frac{1}{2} \left| \left(x_B-x_A\right)\left(y_C-y_A\right) - \left(x_C-x_A\right) \left(y_B-y_A\right) \right|. \]

Dans un repère \( \left(A~;~\vec{i}~,~\vec{j}\right) \) où \( x_A = y_A = 0 \), on a alors \[ S = \frac{1}{2} \left| x_B y_C - x_C y_B \right| . \]

Dès lors, nous avons notre première méthode applicable en langage de programmation, que nous pouvons énoncer ainsi :
Dans un repère d'origine \( A \), si le triplet \( \left( \alpha , \beta , \gamma \right) \) défini par $$\begin{cases} \alpha = \frac{\textrm{Aire}(PBC)}{\textrm{Aire}(ABC)} \\ \beta = \frac{\textrm{Aire}(APC)}{\textrm{Aire}(ABC)} \\ \gamma = \frac{\textrm{Aire}(ABP)}{\textrm{Aire}(ABC)} \end{cases}$$ est tel que \( \alpha + \beta + \gamma = 1 \), alors \( P \in ABC \).

Résoudre un système d'équations linéaires

Cette méthode est plus rapide à expliquer et plus concise, elle demande cependant de savoir résoudre un système de trois équations linéaires à trois inconnues. En effet, nous avons vu que \( \alpha \cdot \overrightarrow{PA} + \beta \cdot \overrightarrow{PB} + \gamma\cdot\overrightarrow{PC} =0 \), donc on en déduit que si \( P \) est à l'intérieur du triangle \( ABC \), cela signifie qu'il existe trois réels positifs \(\alpha\), \(\beta\) et \(\gamma\) tels que \( \alpha + \beta + \gamma = 1 \) et \( \alpha A + \beta B + \gamma C = P \).

Nous avons alors notre seconde méthode :
Si les trois coefficients \( \alpha\), \(\beta\) et \(\gamma\) définis par \[\begin{cases} \alpha x_A + \beta x_B + \gamma x_C= x_P \\ \alpha y_A + \beta y_B + \gamma y_C= y_P \\ \alpha + \beta + \gamma = 1 \\ \end{cases} \] sont positifs, alors \( P \in ABC \).

Pour résoudre ce système, il faut utiliser la méthode d'élimination de Gauss-Jordan. Je ne vais pas l'expliquer, je donnerai simplement le code en C++ qui permet de résoudre ce système et qui utilise cette méthode.

Nous avons nos deux méthodes, on peut maintenant passer à l'application en C++ !

De la programmation : application en C++

Méthode des aires

On rappelle l'expression de la première méthode :
Dans un repère d'origine \( A \), si le triplet \( \left( \alpha , \beta , \gamma \right) \) défini par $$\begin{cases} \alpha = \frac{\textrm{Aire}(PBC)}{\textrm{Aire}(ABC)} \\ \beta = \frac{\textrm{Aire}(APC)}{\textrm{Aire}(ABC)} \\ \gamma = \frac{\textrm{Aire}(ABP)}{\textrm{Aire}(ABC)} \end{cases}$$ est tel que \( \alpha + \beta + \gamma = 1 \), alors \( P \in ABC \).

C'est la méthode la plus simple du point de vue de la programmation. Elle nécessite de calculer simplement trois rapports d'aires qui se calculent en fonction des coordonnées des points considérés et de vérifier une condition sur la somme de ces rapports d'aires.

Voici une proposition d'application en C++ :

#include 
#include 
#include 

using namespace std;

//On définit la strucuture "Point"
struct Point {
    float x;
    float y;
};

//Fonction qui renvoie la lettre de l'alphabet correspondant à son rang
char alphabet(int n);

//Fonction qui calcule l'aire d'un triangle en fonction de ses trois sommets
float aire(Point A, Point B, Point C);

int main() {
    //On définit le tableau "sommets" dans lequel on va placer les sommets du triangle
    vector sommets(3);
    sommets[0].x = 0;
    sommets[0].y = 0;
    cout << "Coordonnées de A : (0;0)." << endl;
    for (int i = 1; i < 3; i++) {
        char nom_point = alphabet(i);
        cout << "Coordonnées de " << nom_point << " :" << endl;
        cout << "x_" << nom_point << " = ";
        cin >> sommets[i].x;
        cout << "y_" << nom_point << " = ";
        cin >> sommets[i].y;
    }

    //On définit A, B et C comme étant les sommets du triangle
    Point A = sommets[0];
    Point B = sommets[1];
    Point C = sommets[2];

    //Soit un point P
    Point P;
    cout << "Coordonnées de P :" << endl;
    cout << "x_P = ";
    cin >> P.x;
    cout << "y_P = ";
    cin >> P.y;

    //On calcule les coordonnées barycentriques α, β et γ du point P
    float alpha = aire(P, B, C) / aire(A, B, C);
    float beta = aire(A, P, C) / aire(A, B, C);
    float gamma = aire(A, B, P) / aire(A, B, C);

    //On affiche les coordonnées barycentriques
    cout << endl << "(α ; β ; γ) = (" << alpha << " ; " << beta << " ; " << gamma << ")" << endl;

    //Si la somme des coordonnées barycentriques vaut 1, alors P appartient à ABC, sinon
    //il n'appartient pas à ABC.
    if (alpha + beta + gamma == 1) {
        cout << "Le point P appartient au triangle ABC." << endl;
    } else {
        cout << "Le point P n'appartient pas au triangle ABC." << endl;
    }

    return EXIT_SUCCESS;
}

float aire(Point A, Point B, Point C) {
    //L'aire S d'un triangle ABC avec A(0;0) peut se calculer comme suit :
    // S = 0.5 * partie_entière(xB * yC − xC * yB).
    return 0.5 * abs(B.x * C.y - C.x * B.y);
}

char alphabet(int n) {
    string alphabet = {'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S',
                       'T', 'U', 'V', 'W', 'X', 'Y', 'Z'};
    return alphabet[n];
}

Méthode du système d'équation

On rappelle l'expression de la seconde méthode :
Si les trois coefficients \( \alpha\), \(\beta\) et \(\gamma\) définis par \[\begin{cases} \alpha x_A + \beta x_B + \gamma x_C= x_P \\ \alpha y_A + \beta y_B + \gamma y_C= y_P \\ \alpha + \beta + \gamma = 1 \\ \end{cases} \] sont positifs, alors \( P \in ABC \).

Nous allons utiliser la méthode d'élimination de Gauss-Jordan afin de déterminer la valeur de ces trois coefficients. S'ils sont supérieurs ou égaux à zéro, alors \( P \) appartient à \( ABC \).

Voici une proposition d'application en C++, dont l'application de la méthode d'élimination de Gauss-Jordan est directement prise de cette page :

#include 
#include 
#include 
#include 
#include 

using namespace std;

//On définit la strucuture "Point"
struct Point {
    float x;
    float y;
};

float * * malloc2 ( int rows , int cols ) ;

//Fonction qui renvoie la lettre de l'alphabet correspondant à son rang
char alphabet(int n);

//Méthode d'élimination de Gauss-Jordan (merci à http://www.student.montefiore.ulg.ac.be/~ggilles/C/gauss.c)
void gauss(float **A, float *b, float *x, int n);

int main() {
    int n ;        // taille du système
    int i, j ;
    float * * A ;  // matrice A
    float * b ;    // vecteur b
    float * x ;    // vecteur des inconnues x

    //On définit le tableau "sommets" dans lequel on va placer les sommets du triangle
    vector sommets(3);
    sommets[0].x = 0;
    sommets[0].y = 0;
    cout << "Coordonnées de A : (0;0)." << endl;
    for (int k = 1; k < 3; k++) {
        char nom_point = alphabet(k);
        cout << "Coordonnées de " << nom_point << " :" << endl;
        cout << "x_" << nom_point << " = ";
        cin >> sommets[k].x;
        cout << "y_" << nom_point << " = ";
        cin >> sommets[k].y;
    }

    //On définit X, Y et Z comme étant les sommets du triangle
    Point X = sommets[0];
    Point Y = sommets[1];
    Point Z = sommets[2];

    //Soit un point P
    Point P;
    cout << "Coordonnées de P :" << endl;
    cout << "x_P = ";
    cin >> P.x;
    cout << "y_P = ";
    cin >> P.y;

    /* Saisie de la taille du système */
    n = 3;

    /* Allocation de mémoire pour A, b et x */
    A = malloc2(n,n) ;
    b = (float *) malloc (sizeof (float *) * n) ;
    x = (float *) malloc (sizeof (float *) * n) ;

    A[0][0] = X.x;
    A[0][1] = Y.x;
    A[0][2] = Z.x;

    A[1][0] = X.y;
    A[1][1] = Y.y;
    A[1][2] = Z.y;

    A[2][0] = 1;
    A[2][1] = 1;
    A[2][2] = 1;

    b[0] = P.x;
    b[1] = P.y;
    b[2] = 1;

    /* Résolution du système par la méthode d'élimination de Gauss */
    gauss(A,b,x,n) ;

    float alpha = x[0];
    float beta = x[1];
    float gamma = x[2];

    cout << endl << "(α ; β ; γ) = (" << alpha << " ; " << beta << " ; " << gamma << ")" << endl;

    if(alpha >=0 && beta >= 0 && gamma >= 0){
        cout << "Le point P appartient au triangle ABC." << endl;
    } else {
        cout << "Le point P appartient au triangle ABC." << endl;
    }

    return EXIT_SUCCESS;
}


char alphabet(int n) {
    string alphabet = {'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S',
                       'T', 'U', 'V', 'W', 'X', 'Y', 'Z'};
    return alphabet[n];
}

float * * malloc2 ( int rows , int cols )
{
    float * * ptr = ( float * * ) malloc ( sizeof ( float * ) * rows + sizeof ( float ) * cols * rows ) ;
    float *   dat = ( float *   ) ( ptr + rows ) ;
    int i ;
    if ( ptr == NULL ) exit( EXIT_FAILURE ) ;
    for ( i = 0 ; i < rows ; ++ i , dat += cols ) ptr [ i ] = dat ;
    return ptr ;
}

/* Méthode d'élimination de Gauss */
void gauss(float **A, float *b, float *x, int n)
{
    int i, j, k ;
    int imin ;
    float p ;
    float sum, valmin, tump1, tump2 ;

    for(k = 0 ; k < n-1 ; k++)
    {
        /* Dans un premier temps, on cherche l'élément minimum (non */
        /* nul) en valeur absolue dans la colonne k et d'indice i   */
        /* supérieur ou égal à k.                                   */

        valmin = A[k][k] ; imin = k ;
        for(i = k+1 ; i < n ; i++)
        {
            if (valmin != 0)
            {
                if (abs(A[i][k]) < abs(valmin) && A[i][k] != 0)
                {
                    valmin = A[i][k] ;
                    imin = i ;
                }
            }
            else
            {
                valmin = A[i][k] ;
                imin = i ;
            }
        }

        /* Si la matrice n'est pas singulière, on inverse    */
        /* les éléments de la ligne imax avec les éléments   */
        /* de la ligne k. On fait de même avec le vecteur b. */

        for(j = 0 ; j < n ; j++)
        {
            tump1 = A[imin][j] ;
            A[imin][j] = A[k][j] ;
            A[k][j] = tump1 ;
        }

        tump2 = b[imin] ;
        b[imin] = b[k] ;
        b[k] = tump2 ;


        /* On procède à la réduction de la matrice par la */
        /* méthode d'élimination de Gauss. */

        for(i = k+1 ; i < n ; i++)
        {
            p = A[i][k]/A[k][k] ;

            for(j = 0 ; j < n ; j++)
            {
                A[i][j] = A[i][j] - p*A[k][j] ;
            }

            b[i] = b[i] - p*b[k] ;
        }
    }

    /* Une fois le système réduit, on obtient une matrice triangulaire */
    /* supérieure et la résolution du système se fait très simplement. */

    x[n-1] = b[n-1]/A[n-1][n-1] ;

    for(i = n-2 ; i > -1 ; i--)
    {
        sum = 0 ;

        for(j = n-1 ; j > i ; j--)
        {
            sum = sum + A[i][j]*x[j] ;
        }
        x[i] = (b[i] - sum)/A[i][i] ;
    }
}

Ce billet touche à sa fin. J'espère qu'il vous aura plu, et peut-être permis de résoudre votre problème.

Travailler avec la SFML dans CLion
15 mai 2017