123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855 |
- /*************************************/
- /* Auteur : Rémi Synave */
- /* Date de création : 01/03/07 */
- /* Date de modification : 18/09/16 */
- /* Version : 0.4 */
- /*************************************/
- /***************************************************************************/
- /* This file is part of a2ri. */
- /* */
- /* a2ri is free software: you can redistribute it and/or modify it */
- /* under the terms of the GNU Lesser General Public License as published */
- /* by the Free Software Foundation, either version 3 of the License, or */
- /* (at your option) any later version. */
- /* */
- /* a2ri is distributed in the hope that it will be useful, */
- /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
- /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
- /* GNU Lesser General Public License for more details. */
- /* */
- /* You should have received a copy of the GNU Lesser General Public */
- /* License along with a2ri. */
- /* If not, see <http://www.gnu.org/licenses/>. */
- /***************************************************************************/
- #include "geometry.h"
- /********** INTERMEDIATE TYPES AND FUNCTIONS **********/
- /* Les fonctions intermédiaires sont préfixées de IF */
- /* et les types intermédiaires de IT */
- void
- IFplanecoeff (
- const point3d * const A,
- const point3d * const B,
- const point3d * const C,
- vector3d * N,
- double *D)
- {
- vector3d AB,
- AC;
- vector3d_init (&AB, B->x - A->x, B->y - A->y, B->z - A->z);
- vector3d_init (&AC, C->x - A->x, C->y - A->y, C->z - A->z);
- *N = vector3d_vectorialproduct (&AB, &AC);
- vector3d_init (&AB, A->x, A->y, A->z);
- *D = vector3d_scalarproduct (&AB, N);
- }
- /********** MAIN FUNCTIONS **********/
- /**
- Calcul de la longueur entre deux point3d
- @param p1 premier point3d
- @param p2 second point3d
- @return longueur \f$p_1p_2\f$
- */
- double
- point3d_length (
- const point3d * const p1,
- const point3d * const p2)
- {
- return sqrt (point3d_square_length (p1, p2));
- }
- /**
- Calcul de la longueur au carré entre deux point3d
- @param p1 premier point3d
- @param p2 second point3d
- @return longueur \f$p_1p_2\f$
- */
- double
- point3d_square_length (
- const point3d * const p1,
- const point3d * const p2)
- {
- return (p2->x - p1->x) * (p2->x - p1->x) + (p2->y - p1->y) * (p2->y - p1->y) +
- (p2->z - p1->z) * (p2->z - p1->z);
- }
- /**
- Calcul de la longueur entre deux point2d
- @param p1 premier point2d
- @param p2 second point2d
- @return longueur \f$p_1p_2\f$
- */
- double
- point2d_length (
- const point2d * const p1,
- const point2d * const p2)
- {
- return sqrt ((p2->x - p1->x) * (p2->x - p1->x) + (p2->y - p1->y) * (p2->y - p1->y));
- }
- /**
- Calcule l'équation d'un plan à partir d'une liste de point sous la forme ax+by+cz+d=0
- @param p liste de points
- @param size nombre de points
- @param aeq valeur de "a"
- @param beq valeur de "b"
- @param ceq valeur de "c"
- @param deq valeur de "d"
- @return aucun
- Exemple de code
- @code
- double a,b,c,d;
- point3d p[3];
- point3d_init(&(p[0]),0.0,0.0,0.0);
- point3d_init(&(p[1]),1.0,1.0,-2.0);
- point3d_init(&(p[2]),-1.0,1.0,0.0);
- equation_plan(p,3,&a,&b,&c,&d);
- printf("%dx+%dy+%dz+%d=0\n",a,b,c,d);
- @endcode
- donne le résultat
- @code
- 1.0x+1.0y+1.0z+0.0=0
- @endcode
- */
- void
- equation_plan (
- const point3d * const p,
- int size,
- double *aeq,
- double *beq,
- double *ceq,
- double *deq)
- {
- double a = 0;
- double b = 0;
- double c = 0;
- double d = 0;
- double xav = 0;
- double yav = 0;
- double zav = 0;
- for (int i = 0; i < size; i++)
- {
- a += (p[i].y - p[(i + 1) % size].y) * (p[i].z + p[(i + 1) % size].z);
- b += (p[i].z - p[(i + 1) % size].z) * (p[i].x + p[(i + 1) % size].x);
- c += (p[i].x - p[(i + 1) % size].x) * (p[i].y + p[(i + 1) % size].y);
- xav += p[i].x;
- yav += p[i].y;
- zav += p[i].z;
- }
- xav /= (size * 1.0);
- yav /= (size * 1.0);
- zav /= (size * 1.0);
- d = -(xav * a + yav * b + zav * c);
- *aeq = a;
- *beq = b;
- *ceq = c;
- *deq = d;
- }
- /**
- Calcul de l'aire du triangle défini par les trois points
- @param p1 premier point du triangle
- @param p2 second point
- @param p3 troisième point
- @return aire du triangle
- */
- double
- point3d_area (
- const point3d * const p1,
- const point3d * const p2,
- const point3d * const p3)
- {
- vector3d AB,
- AC,
- U;
- vector3d_init (&AB, p2->x - p1->x, p2->y - p1->y, p2->z - p1->z);
- vector3d_init (&AC, p3->x - p1->x, p3->y - p1->y, p3->z - p1->z);
- U = vector3d_vectorialproduct (&AB, &AC);
- //aire du triangle = 0.5 * norme du produit vectoriel AB AC
- return (0.5 * vector3d_size (&U));
- }
- /**
- Transformation des coordonnées d'un point dans un espace à trois dimension en coordonnées dans un espace à deux dimensions dont la base est le couple (base1,base2) et d'origine le point3d "origin"
- @param p point dont les coordonnées vont etre recalculées
- @param origin point origine du nouveau repére
- @param base1 premier vecteur formant la nouvelle base
- @param base2 second vecteur formant la nouvelle base
- @param newp point contenant les coordonnées dans la nouvelle base
- @return aucun
- */
- void
- base_modification_3d_to_2d (
- const point3d * const p,
- const point3d * const origin,
- const vector3d * const base1,
- const vector3d * const base2,
- point2d * newp)
- {
- //résolution du système de trois équations à deux inconnues grace à la méthode des moindres carrés
- point3d ptemp;
- point3d_init(&ptemp,p->x-origin->x,p->y-origin->y,p->z-origin->z);
- double left_data[6];
- double right_data[3];
- left_data[0] = base1->dx;
- left_data[1] = base2->dx;
- left_data[2] = base1->dy;
- left_data[3] = base2->dy;
- left_data[4] = base1->dz;
- left_data[5] = base2->dz;
- right_data[0] = ptemp.x;
- right_data[1] = ptemp.y;
- right_data[2] = ptemp.z;
- gsl_matrix *A = matrix_init (left_data, 3, 2);
- gsl_matrix *Y = matrix_init (right_data, 3, 1);
- gsl_matrix *transA = matrix_transpose (A);
- gsl_matrix *multransAA = matrix_mul (transA, A);
- gsl_matrix *invmultransAA = matrix_inverse (multransAA);
- gsl_matrix *pseudoinverse = matrix_mul (invmultransAA, transA);
- gsl_matrix *inc = matrix_mul (pseudoinverse, Y);
- newp->x = gsl_matrix_get (inc, 0, 0);
- newp->y = gsl_matrix_get (inc, 1, 0);
- gsl_matrix_free (A);
- gsl_matrix_free (Y);
- gsl_matrix_free (pseudoinverse);
- gsl_matrix_free (inc);
- gsl_matrix_free (transA);
- gsl_matrix_free (multransAA);
- gsl_matrix_free (invmultransAA);
- }
- /**
- Calcul du vecteur normal à AB et contenu dans le meme plan que celui formé par les vecteurs (AB,AC). Le résultat final est une base orthogonal du repére.
- @param AB premier vecteur de la base
- @param AC second vecteur du repére
- @param U pointeur sur le nouveau vecteur formant la base orthogonal
- */
- void
- find_second_base_vector (
- const vector3d * const AB,
- const vector3d * const AC,
- vector3d * U)
- {
- vector3d W,
- Utemp;
- W = vector3d_vectorialproduct (AB, AC);
- Utemp = vector3d_vectorialproduct (&W, AB);
- U->dx = Utemp.dx;
- U->dy = Utemp.dy;
- U->dz = Utemp.dz;
- }
- /**
- Vérification qu'un point M se trouve bien dans le triangle ABC
- @param M le point à tester
- @param A premier sommet du triangle
- @param B second sommet du triangle
- @param C troisième sommet du triangle
- @return 1 si le point M appartient au triangle ABC, 0 sinon
- */
- int
- point_in_triangle (
- const point3d * const M,
- const point3d * const A,
- const point3d * const B,
- const point3d * const C)
- {
- double sommeangle = 0;
- vector3d MA,
- MB,
- MC;
- vector3d_init (&MA, A->x - M->x, A->y - M->y, A->z - M->z);
- vector3d_init (&MB, B->x - M->x, B->y - M->y, B->z - M->z);
- vector3d_init (&MC, C->x - M->x, C->y - M->y, C->z - M->z);
- sommeangle += vector3d_angle_degre (&MA, &MB);
- sommeangle += vector3d_angle_degre (&MB, &MC);
- sommeangle += vector3d_angle_degre (&MC, &MA);
- //le point est contenu dans le triangle si la somme des angles est égale à 360
- return (egalite (sommeangle, 360));
- }
- /**
- Calcule l'angle en radian formé par les deux vector3d
- @param v1 le premier vector3d
- @param v2 le second vector3d
- @return angle en radian
- */
- double
- vector3d_angle_radian (
- const vector3d * const v1,
- const vector3d * const v2)
- {
- vector3d v1temp,v2temp;
- vector3d_init(&v1temp,v1->dx,v1->dy,v1->dz);
- vector3d_init(&v2temp,v2->dx,v2->dy,v2->dz);
- vector3d_normalize (&v1temp);
- vector3d_normalize (&v2temp);
- double prodscal = vector3d_scalarproduct (&v1temp, &v2temp);
- if (prodscal < -1)
- prodscal = -1;
- if (prodscal > 1)
- prodscal = 1;
- return acos (prodscal);
- }
- /**
- Calcule l'angle en degré formé par les deux vector3d
- @param v1 le premier vector3d
- @param v2 le second vector3d
- @return angle en degré
- */
- double
- vector3d_angle_degre (
- const vector3d * const v1,
- const vector3d * const v2)
- {
- vector3d v1temp,v2temp;
- vector3d_init(&v1temp,v1->dx,v1->dy,v1->dz);
- vector3d_init(&v2temp,v2->dx,v2->dy,v2->dz);
- vector3d_normalize (&v1temp);
- vector3d_normalize (&v2temp);
- double prodscal = vector3d_scalarproduct (&v1temp, &v2temp);
- if (prodscal < -1)
- prodscal = -1;
- if (prodscal > 1)
- prodscal = 1;
- return (acos (prodscal)) * (360 / (2 * PI));
- }
- /**
- Calcule le centre du cercle inscrit à un triangle
- @param A premier point du triangle
- @param B second point
- @param C troisème point
- @return le point centre du cercle inscrit
- */
- point3d
- incircle_center (
- const point3d * const A,
- const point3d * const B,
- const point3d * const C)
- {
- point3d center;
- double a = point3d_length (B, C);
- double b = point3d_length (A, C);
- double c = point3d_length (A, B);
- center.x = (a * A->x + b * B->x + c * C->x) / (a + b + c);
- center.y = (a * A->y + b * B->y + c * C->y) / (a + b + c);
- center.z = (a * A->z + b * B->z + c * C->z) / (a + b + c);
- center.att_int = 0;
- center.att_double = 0;
- return center;
- }
- /**
- Calcule le centre du cercle circonscrit à un triangle
- @param A premier point du triangle
- @param B second point
- @param C troisème point
- @return le point centre du cercle inscrit
- */
- point3d
- circumcircle_center (
- point3d * A,
- point3d * B,
- point3d * C)
- {
- /** méthode du bouquin : 3D Math Pimer for Graphics and Game Development **/
- vector3d e1,
- e2,
- e3,
- me1,
- me2,
- me3;
- double d1,
- d2,
- d3,
- c1,
- c2,
- c3,
- c,
- alpha,
- beta,
- gamma;
- point3d centre;
- vector3d_init (&e1, C->x - B->x, C->y - B->y, C->z - B->z);
- vector3d_init (&e2, A->x - C->x, A->y - C->y, A->z - C->z);
- vector3d_init (&e3, B->x - A->x, B->y - A->y, B->z - A->z);
- vector3d_init (&me1, B->x - C->x, B->y - C->y, B->z - C->z);
- vector3d_init (&me2, C->x - A->x, C->y - A->y, C->z - A->z);
- vector3d_init (&me3, A->x - B->x, A->y - B->y, A->z - B->z);
- d1 = vector3d_scalarproduct (&me2, &e3);
- d2 = vector3d_scalarproduct (&me3, &e1);
- d3 = vector3d_scalarproduct (&me1, &e2);
- c1 = d2 * d3;
- c2 = d3 * d1;
- c3 = d1 * d2;
- c = c1 + c2 + c3;
- alpha = (c2 + c3) / (2.0 * c);
- beta = (c3 + c1) / (2.0 * c);
- gamma = (c1 + c2) / (2.0 * c);
- centre.x = (A->x * alpha + B->x * beta + C->x * gamma);
- centre.y = (A->y * alpha + B->y * beta + C->y * gamma);
- centre.z = (A->z * alpha + B->z * beta + C->z * gamma);
- centre.att_double = 0;
- centre.att_int = 0;
- return centre;
- }
- /**
- Calcule le centre de sphère circonscrite au tétraédre
- @param A premier point du tétraédre
- @param B second point
- @param C troisème point
- @param D quatrième point
- @return le point centre de la sphère circonscrite
- */
- point3d
- circumsphere_center (
- point3d * A,
- point3d * B,
- point3d * C,
- point3d * D)
- {
- double dataL[9] = { 2 * (B->x - A->x), 2 * (B->y - A->y), 2 * (B->z - A->z),
- 2 * (C->x - A->x), 2 * (C->y - A->y), 2 * (C->z - A->z),
- 2 * (D->x - A->x), 2 * (D->y - A->y), 2 * (D->z - A->z)
- };
- double dataR[3] = { point3d_length (A, B) * point3d_length (A, B),
- point3d_length (A, C) * point3d_length (A, C),
- point3d_length (A, D) * point3d_length (A, D)
- };
- gsl_matrix *matA = matrix_init (dataL, 3, 3);
- gsl_matrix *Y = matrix_init (dataR, 3, 1);
- gsl_matrix *pseudoinverse =
- matrix_mul (matrix_inverse (matrix_mul (matrix_transpose (matA), matA)),
- matrix_transpose (matA));
- gsl_matrix *inc = matrix_mul (pseudoinverse, Y);
- point3d centre;
- point3d_init (¢re, A->x + gsl_matrix_get (inc, 0, 0),
- A->y + gsl_matrix_get (inc, 1, 0), A->z + gsl_matrix_get (inc,
- 2,
- 0));
- return centre;
- }
- /**
- Calcule la distance entre un point M et une droite définie par deux points A et B
- @param M point
- @param A premier point sur la droite
- @param B second point sur la droite
- @return distance entre le point M et la droite AB
- */
- double
- distance_point_straight_line (
- point3d * M,
- point3d * A,
- point3d * B)
- {
- double aire = point3d_area (M, A, B);
- return (2.0 * aire) / point3d_length (A, B);
- }
- /**
- Calcule la distance entre un point M et un plan définie par trois points A,b et C
- @param M point
- @param A premier point du plan
- @param B second point du plan
- @param C troisieme point du plan
- @return distance entre le point M et le plan ABC
- */
- double
- distance_point_plane (
- point3d * M,
- point3d * A,
- point3d * B,
- point3d * C)
- {
- double a,
- b,
- c,
- d;
- point3d p[3];
- point3d_init (&p[0], A->x, A->y, A->z);
- point3d_init (&p[1], B->x, B->y, B->z);
- point3d_init (&p[2], C->x, C->y, C->z);
- equation_plan (p, 3, &a, &b, &c, &d);
- return (fabs (a * M->x + b * M->y + c * M->z + d) /
- (sqrt (a * a + b * b + c * c)));
- }
- /**
- Calcule la distance entre un point M et un triangle défini par trois points A,b et C
- @param M point
- @param A premier point du triangle
- @param B second point du triangle
- @param C troisieme point du triangle
- @return distance entre le point M et le triangle ABC
- */
- double
- distance_point_triangle (
- point3d * M,
- point3d * A,
- point3d * B,
- point3d * C)
- {
- vector3d AB,
- AC,
- U,
- V;
- point2d Mb,
- Ab,
- Bb,
- Cb;
- point3d Mc,
- Ac,
- Bc,
- Cc;
- double d[3];
- vector3d_init (&AB, B->x - A->x, B->y - A->y, B->z - A->z);
- vector3d_init (&AC, C->x - A->x, C->y - A->y, C->z - A->z);
- vector3d_normalize (&AB);
- U = AB;
- find_second_base_vector (&AB, &AC, &V);
- base_modification_3d_to_2d (M, A, &U, &V, &Mb);
- base_modification_3d_to_2d (A, A, &U, &V, &Ab);
- base_modification_3d_to_2d (B, A, &U, &V, &Bb);
- base_modification_3d_to_2d (C, A, &U, &V, &Cb);
- point3d_init (&Mc, Mb.x, Mb.y, 0);
- point3d_init (&Ac, Ab.x, Ab.y, 0);
- point3d_init (&Bc, Bb.x, Bb.y, 0);
- point3d_init (&Cc, Cb.x, Cb.y, 0);
- if (point_in_triangle (&Mc, &Ac, &Bc, &Cc))
- return distance_point_plane (M, A, B, C);
- else
- {
- d[0] = point3d_length (M, A);
- d[1] = point3d_length (M, B);
- d[2] = point3d_length (M, C);
- list_double_sort (d, 3, ASC);
- return d[0];
- }
- }
- /**
- Echantillonne le triangle ABC
- @param A point A du triangle
- @param B point B
- @param C point C
- @param nb_sample nombre d'échantillon à extraire
- @param list tableau de réel représentant les coordonées des échantillons
- @param size taille du tableau
- @return aucun
- */
- void
- sample_triangle (
- const point3d * const A,
- const point3d * const B,
- const point3d * const C,
- int nb_sample,
- point3d ** list)
- {
- vector3d AB,
- AC;
- double alea1,
- alea2,
- alea3,
- somme,
- echx,
- echy,
- echz;
- *list = (point3d *) malloc (nb_sample * sizeof (point3d));
- a2ri_erreur_critique_si (*list == NULL,
- "erreur allocation memoire pour list\nsample_triangle");
- srand (time (NULL));
- vector3d_init (&AB, B->x - A->x, B->y - A->y, B->z - A->z);
- vector3d_init (&AC, C->x - A->x, C->y - A->y, C->z - A->z);
- for (int i = 0; i < nb_sample; i++)
- {
- alea1 = (rand () % 1000) / 1000.0;
- alea2 = (rand () % 1000) / 1000.0;
- alea3 = (rand () % 1000) / 1000.0;
- somme = alea1 + alea2 + alea3;
- alea1 /= somme;
- alea2 /= somme;
- alea3 /= somme;
- echx = A->x + alea1 * AB.dx + alea1 * AC.dx;
- echy = A->y + alea2 * AB.dy + alea2 * AC.dy;
- echz = A->z + alea3 * AB.dz + alea3 * AC.dz;
- (*list)[i].x = echx;
- (*list)[i].y = echy;
- (*list)[i].z = echz;
- }
- }
- /**
- Test de l'intersection entre une droite et un plan
- @param d1 premier point définissant la droite
- @param d2 second point définissant la droite
- @param p1 premier point du plan
- @param p2 second point du plan
- @param p3 troisieme point du plan
- @param t paramètre d'intersection
- @return 1 si intersection, 0 sinon
- */
- int
- intersection_droite_plan (
- point3d * d1,
- point3d * d2,
- point3d * p1,
- point3d * p2,
- point3d * p3,
- double *t)
- {
- vector3d N,
- rq;
- double D,
- num,
- denom;
- IFplanecoeff (p1, p2, p3, &N, &D);
- vector3d_init (&rq, d1->x, d1->y, d1->z);
- num = D - vector3d_scalarproduct (&rq, &N);
- vector3d_init (&rq, d2->x - d1->x, d2->y - d1->y, d2->z - d1->z);
- denom = vector3d_scalarproduct (&rq, &N);
- if (denom == 0.0)
- {
- if (num == 0.0)
- {
- *t = 0.0;
- return 1;
- }
- else
- return 0;
- }
- *t = num / denom;
- return 1;
- }
- /**
- Test de l'intersection entre une droite et un triangle
- @param d1 premier point définissant la droite
- @param d2 second point définissant la droite
- @param p1 premier point du triangle
- @param p2 second point du triangle
- @param p3 troisieme point du triangle
- @return 1 si intersection, 0 sinon
- */
- int
- intersection_droite_triangle (
- point3d * d1,
- point3d * d2,
- point3d * p1,
- point3d * p2,
- point3d * p3,
- double *t)
- {
- point3d inter;
- intersection_droite_plan (d1, d2, p1, p2, p3, t);
- point3d_init (&inter, d1->x + (*t) * (d2->x - d1->x),
- d1->y + (*t) * (d2->y - d1->y), d1->z + (*t) * (d2->z - d1->z));
- return point_in_triangle (&inter, p1, p2, p3);
- }
- /**
- Test de l'intersection entre un segment et un triangle
- @param s1 premier point définissant le segment
- @param s2 second point définissant le segment
- @param p1 premier point du triangle
- @param p2 second point du triangle
- @param p3 troisieme point du triangle
- @return 1 si intersection, 0 sinon
- */
- int
- intersection_segment_triangle (
- point3d * s1,
- point3d * s2,
- point3d * p1,
- point3d * p2,
- point3d * p3,
- double *t)
- {
- point3d inter;
- intersection_droite_plan (s1, s2, p1, p2, p3, t);
- if(*t>=0 && *t<=1.0){
- point3d_init (&inter, s1->x + (*t) * (s2->x - s1->x),
- s1->y + (*t) * (s2->y - s1->y), s1->z + (*t) * (s2->z - s1->z));
-
- return point_in_triangle (&inter, p1, p2, p3);
- }
- return 0;
- }
- /**
- Test de l'intersection entre deux triangles
- @param t1 premier point définissant le premier triangle
- @param t2 second point définissant le premier triangle
- @param t3 troisieme point définissant le premier triangle
- @param p1 premier point du second triangle
- @param p2 second point du second triangle
- @param p3 troisieme point du second triangle
- @return 1 si intersection, 0 sinon
- */
- int
- intersection_triangle_triangle (
- point3d * t1,
- point3d * t2,
- point3d * t3,
- point3d * p1,
- point3d * p2,
- point3d * p3)
- {
- double t;
- if (intersection_droite_triangle (t1, t2, p1, p2, p3, &t))
- {
- if (t >= 0 && t <= 1)
- return 1;
- }
- if (intersection_droite_triangle (t1, t3, p1, p2, p3, &t))
- {
- if (t >= 0 && t <= 1)
- return 1;
- }
- if (intersection_droite_triangle (t2, t3, p1, p2, p3, &t))
- {
- if (t >= 0 && t <= 1)
- return 1;
- }
- if (intersection_droite_triangle (p1, p2, t1, t2, t3, &t))
- {
- if (t >= 0 && t <= 1)
- return 1;
- }
- if (intersection_droite_triangle (p1, p3, t1, t2, t3, &t))
- {
- if (t >= 0 && t <= 1)
- return 1;
- }
- if (intersection_droite_triangle (p2, p3, t1, t2, t3, &t))
- {
- if (t >= 0 && t <= 1)
- return 1;
- }
-
- return 0;
- }
- /**
- Teste si les trois points A, B et C sont alignés
- @param A premier point
- @param B second point
- @param C troisieme point
- @return 1 si les trois points sont alignés, 0 sinon
- */
- int
- trois_points_alignes (
- point3d * A,
- point3d * B,
- point3d * C)
- {
- double t;
- if (point3d_equal (A, B))
- return 1;
- if (point3d_equal (A, C))
- return 1;
- if (point3d_equal (B, C))
- return 1;
- if (B->x - A->x != 0)
- t = (C->x - A->x) / (B->x - A->x);
- else
- {
- if (B->y - A->y != 0)
- t = (C->y - A->y) / (B->y - A->y);
- else
- t = (C->z - A->z) / (B->z - A->z);
- }
- if (!egalite (A->x + t * (B->x - A->x), C->x))
- return 0;
- if (!egalite (A->y + t * (B->y - A->y), C->y))
- return 0;
- if (!egalite (A->z + t * (B->z - A->z), C->z))
- return 0;
- return 1;
- }
|