123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625 |
- /*************************************/
- /* Auteur : Rémi Synave */
- /* Date de création : 01/03/07 */
- /* Date de modification : 15/03/15 */
- /* 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 "boundingbox.h"
- /********** INTERMEDIATE TYPES AND FUNCTIONS **********/
- /* Les fonctions intermédiaires sont préfixées de IF */
- /* et les types intermédiaires de IT */
- /********** MAIN FUNCTIONS **********/
- /**
- Calcul des boites englobantes (parallélépipèdes) des faces de la liste
- @param m le modèle
- @param faces tableau contenant les numéros des faces pour lesquels on veut les boites englobantes
- @param nbfaces taille du tableau
- @param ptmin pointeur sur le tableau contenant les points (xmin,ymin,zmin) des boites englobantes
- @param ptmax pointeur sur le tableau contenant les points (xmax,ymax,zmax) des boites englobantes
- **/
- void
- a2ri_vf_axis_aligned_bounding_box (
- const vf_model * const m,
- const int * const faces,
- int nbfaces,
- point3d ** ptmin,
- point3d ** ptmax)
- {
- //allocation de l'espace pour les tableaux ptmin et ptmax
- (*ptmin) = (point3d *) malloc (nbfaces * sizeof (point3d));
- a2ri_erreur_critique_si (*ptmin == NULL,
- "erreur allocation memoire pour ptmin\na2ri_vf_axis_aligned_bounding_box");
- (*ptmax) = (point3d *) malloc (nbfaces * sizeof (point3d));
- a2ri_erreur_critique_si (*ptmax == NULL,
- "erreur allocation memoire pour ptmax\na2ri_vf_axis_aligned_bounding_box");
- for (int i = 0; i < nbfaces; i++)
- {
- //pour la face i : renseignement des coins ptmin et ptmax avec le sommet ve1
- (*ptmin)[faces[i]].x = m->ve[m->fa[faces[i]].ve1].x;
- (*ptmin)[faces[i]].y = m->ve[m->fa[faces[i]].ve1].y;
- (*ptmin)[faces[i]].z = m->ve[m->fa[faces[i]].ve1].z;
- (*ptmax)[faces[i]].x = m->ve[m->fa[faces[i]].ve1].x;
- (*ptmax)[faces[i]].y = m->ve[m->fa[faces[i]].ve1].y;
- (*ptmax)[faces[i]].z = m->ve[m->fa[faces[i]].ve1].z;
- //mise à jour si nécessaire avec ve2 et ve3
- if (m->ve[m->fa[faces[i]].ve2].x < (*ptmin)[faces[i]].x)
- (*ptmin)[faces[i]].x = m->ve[m->fa[faces[i]].ve2].x;
- if (m->ve[m->fa[faces[i]].ve2].y < (*ptmin)[faces[i]].y)
- (*ptmin)[faces[i]].y = m->ve[m->fa[faces[i]].ve2].y;
- if (m->ve[m->fa[faces[i]].ve2].z < (*ptmin)[faces[i]].z)
- (*ptmin)[faces[i]].z = m->ve[m->fa[faces[i]].ve2].z;
- if (m->ve[m->fa[faces[i]].ve2].x > (*ptmax)[faces[i]].x)
- (*ptmax)[faces[i]].x = m->ve[m->fa[faces[i]].ve2].x;
- if (m->ve[m->fa[faces[i]].ve2].y > (*ptmax)[faces[i]].y)
- (*ptmax)[faces[i]].y = m->ve[m->fa[faces[i]].ve2].y;
- if (m->ve[m->fa[faces[i]].ve2].z > (*ptmax)[faces[i]].z)
- (*ptmax)[faces[i]].z = m->ve[m->fa[faces[i]].ve2].z;
- if (m->ve[m->fa[faces[i]].ve3].x < (*ptmin)[faces[i]].x)
- (*ptmin)[faces[i]].x = m->ve[m->fa[faces[i]].ve3].x;
- if (m->ve[m->fa[faces[i]].ve3].y < (*ptmin)[faces[i]].y)
- (*ptmin)[faces[i]].y = m->ve[m->fa[faces[i]].ve3].y;
- if (m->ve[m->fa[faces[i]].ve3].z < (*ptmin)[faces[i]].z)
- (*ptmin)[faces[i]].z = m->ve[m->fa[faces[i]].ve3].z;
- if (m->ve[m->fa[faces[i]].ve3].x > (*ptmax)[faces[i]].x)
- (*ptmax)[faces[i]].x = m->ve[m->fa[faces[i]].ve3].x;
- if (m->ve[m->fa[faces[i]].ve3].y > (*ptmax)[faces[i]].y)
- (*ptmax)[faces[i]].y = m->ve[m->fa[faces[i]].ve3].y;
- if (m->ve[m->fa[faces[i]].ve3].z > (*ptmax)[faces[i]].z)
- (*ptmax)[faces[i]].z = m->ve[m->fa[faces[i]].ve3].z;
- }
- }
- /**
- Calcul des sphères englobantes des faces de la liste en utilisant les fastball de Ritter
- @param m le modèle
- @param faces tableau contenant les numéros des faces pour lesquels on veut les sphères englobantes
- @param nbfaces taille du tableau
- @param listcentre pointeur sur le tableau de point3d contenant les centres de sphères englobantes
- @param listradius pointeur sur le tableau contenant les rayons des sphères englobantes
- **/
- void
- a2ri_vf_bounding_ball_ritter (
- const vf_model * const m,
- const int * const faces,
- int nbfaces,
- point3d ** listcentre,
- double **listradius)
- {
- int nbelt = 0;
- //allocation de l'espae pour les centres des sphères
- *listcentre = (point3d *) malloc (nbfaces * sizeof (point3d));
- a2ri_erreur_critique_si (*listcentre == NULL,
- "erreur allocation memoire pour listcentre\na2ri_vf_bounding_ball_ritter");
- for (int i = 0; i < nbfaces; i++)
- {
- //pour la face faces[i]
- int ve1 = m->fa[faces[i]].ve1;
- int ve2 = m->fa[faces[i]].ve2;
- int ve3 = m->fa[faces[i]].ve3;
- //calcul de la boite englobgante
- double xmin = m->ve[ve1].x;
- double xmax = xmin;
- double ymin = m->ve[ve1].y;
- double ymax = ymin;
- double zmin = m->ve[ve1].z;
- double zmax = zmin;
- if (m->ve[ve2].x < xmin)
- xmin = m->ve[ve2].x;
- if (m->ve[ve2].x > xmax)
- xmax = m->ve[ve2].x;
- if (m->ve[ve3].x < xmin)
- xmin = m->ve[ve3].x;
- if (m->ve[ve3].x > xmax)
- xmax = m->ve[ve3].x;
- if (m->ve[ve2].y < ymin)
- ymin = m->ve[ve2].y;
- if (m->ve[ve2].y > ymax)
- ymax = m->ve[ve2].y;
- if (m->ve[ve3].y < ymin)
- ymin = m->ve[ve3].y;
- if (m->ve[ve3].y > ymax)
- ymax = m->ve[ve3].y;
- if (m->ve[ve2].z < zmin)
- zmin = m->ve[ve2].z;
- if (m->ve[ve2].z > zmax)
- zmax = m->ve[ve2].z;
- if (m->ve[ve3].z < zmin)
- zmin = m->ve[ve3].z;
- if (m->ve[ve3].z > zmax)
- zmax = m->ve[ve3].z;
- //calcul du centre de la boite englobante
- double centrex = (xmin + xmax) / 2.0;
- double centrey = (ymin + ymax) / 2.0;
- double centrez = (zmin + zmax) / 2.0;
- (*listcentre)[i].x = centrex;
- (*listcentre)[i].y = centrey;
- (*listcentre)[i].z = centrez;
- //calcul des trois rayons -> centre-ve1, centre-ve2, centre-ve3
- double radius1 =
- sqrt ((centrex - m->ve[ve1].x) * (centrex - m->ve[ve1].x) +
- (centrey - m->ve[ve1].y) * (centrey - m->ve[ve1].y) + (centrez -
- m->ve[ve1].
- z) *
- (centrez - m->ve[ve1].z));
- double radius2 =
- sqrt ((centrex - m->ve[ve2].x) * (centrex - m->ve[ve2].x) +
- (centrey - m->ve[ve2].y) * (centrey - m->ve[ve2].y) + (centrez -
- m->ve[ve2].
- z) *
- (centrez - m->ve[ve2].z));
- double radius3 =
- sqrt ((centrex - m->ve[ve3].x) * (centrex - m->ve[ve3].x) +
- (centrey - m->ve[ve3].y) * (centrey - m->ve[ve3].y) + (centrez -
- m->ve[ve3].
- z) *
- (centrez - m->ve[ve3].z));
- if (radius2 > radius1)
- radius1 = radius2;
- if (radius3 > radius1)
- radius1 = radius3;
- //ajout du rayon maximal à la liste des rayons
- list_double_add (listradius, &nbelt, radius1, WITH_REDUNDANCE);
- }
- }
- /**
- Calcul des sphères englobantes des faces de la liste
- @param m le modèle
- @param faces tableau contenant les numéros des faces pour lesquels on veut les sphères englobantes
- @param nbfaces taille du tableau
- @param listcentre pointeur sur le tableau de point3d contenant les centres de sphères englobantes
- @param listradius pointeur sur le tableau contenant les rayons des sphères englobantes
- **/
- void
- a2ri_vf_bounding_ball_minimale (
- const vf_model * const m,
- const int * const faces,
- int nbfaces,
- point3d ** listcentre,
- double **listradius)
- {
- double radius,
- x,
- y,
- z;
- int ve1,
- ve2,
- ve3,
- nbelt = 0;
- point3d p3d_1,
- p3d_2,
- p3d_3,
- p3d_centre;
- double longar[3],
- longp1p2,
- longp1p3,
- longp2p3;
- free (*listcentre);
- *listcentre = (point3d *) malloc (nbfaces * sizeof (point3d));
- a2ri_erreur_critique_si (*listcentre == NULL,
- "erreur allocation memoire pour listcentre\na2ri_vf_bounding_ball_minimale");
- for (int alpha = 0; alpha < nbfaces; alpha++)
- {
- ve1 = m->fa[faces[alpha]].ve1;
- ve2 = m->fa[faces[alpha]].ve2;
- ve3 = m->fa[faces[alpha]].ve3;
- point3d_init (&p3d_1, m->ve[ve1].x, m->ve[ve1].y, m->ve[ve1].z);
- point3d_init (&p3d_2, m->ve[ve2].x, m->ve[ve2].y, m->ve[ve2].z);
- point3d_init (&p3d_3, m->ve[ve3].x, m->ve[ve3].y, m->ve[ve3].z);
- longar[0] = point3d_square_length (&p3d_1, &p3d_2);
- longp1p2 = longar[0];
- longar[1] = point3d_square_length (&p3d_1, &p3d_3);
- longp1p3 = longar[1];
- longar[2] = point3d_square_length (&p3d_2, &p3d_3);
- longp2p3 = longar[2];
- list_double_sort (longar, 3, ASC);
- if (longar[2] > longar[0] + longar[1])
- {
- if (longar[2] == longp1p2)
- {
- (*listcentre)[alpha].x = ((p3d_1.x + p3d_2.x) * 1.0) / 2.0;
- (*listcentre)[alpha].y = ((p3d_1.y + p3d_2.y) * 1.0) / 2.0;
- (*listcentre)[alpha].z = ((p3d_1.z + p3d_2.z) * 1.0) / 2.0;
- }
- if (longar[2] == longp1p3)
- {
- (*listcentre)[alpha].x = ((p3d_1.x + p3d_3.x) * 1.0) / 2.0;
- (*listcentre)[alpha].y = ((p3d_1.y + p3d_3.y) * 1.0) / 2.0;
- (*listcentre)[alpha].z = ((p3d_1.z + p3d_3.z) * 1.0) / 2.0;
- }
- if (longar[2] == longp2p3)
- {
- (*listcentre)[alpha].x = ((p3d_2.x + p3d_3.x) * 1.0) / 2.0;
- (*listcentre)[alpha].y = ((p3d_2.y + p3d_3.y) * 1.0) / 2.0;
- (*listcentre)[alpha].z = ((p3d_2.z + p3d_3.z) * 1.0) / 2.0;
- }
- list_double_add (listradius, &nbelt, sqrt (longar[2]) / 2.0,
- WITH_REDUNDANCE);
- }
- else
- {
- p3d_centre = circumcircle_center (&p3d_1, &p3d_2, &p3d_3);
- x = p3d_centre.x;
- y = p3d_centre.y;
- z = p3d_centre.z;
- radius = sqrt ((m->ve[ve1].x - x) * (m->ve[ve1].x - x) +
- (m->ve[ve1].y - y) * (m->ve[ve1].y - y) +
- (m->ve[ve1].z - z) * (m->ve[ve1].z - z));
- (*listcentre)[alpha].x = x;
- (*listcentre)[alpha].y = y;
- (*listcentre)[alpha].z = z;
- list_double_add (listradius, &nbelt, radius, WITH_REDUNDANCE);
- }
- }
- }
- /**
- Calcul l'oriented bounding box du maillage
- @param m le modèle
- @param newbasex pointeur sur le vecteur "X" de la nouvelle base dans laquelle l'oriented bounding box est représenté
- @param newbasey pointeur sur le vecteur "Y" de la nouvelle base dans laquelle l'oriented bounding box est représenté
- @param newbasez pointeur sur le vecteur "Z" de la nouvelle base dans laquelle l'oriented bounding box est représenté
- @param ptmin pointeur sur le point3d contenant le coin inférieur de l'oriented bounding box
- @param ptmax pointeur sur le point3d contenant le coin supérieur de l'oriented bounding box
- **/
- void
- a2ri_vf_oriented_bounding_box (
- const vf_model * const m,
- vector3d * newbasex,
- vector3d * newbasey,
- vector3d * newbasez,
- point3d * ptmin,
- point3d * ptmax)
- {
- point3d *listsommet = NULL;
- gsl_matrix *covariance,
- *matpass = NULL,
- *point,
- *res = NULL;
- double data[9];
- gsl_vector *eval = gsl_vector_alloc (3);
- a2ri_erreur_critique_si (eval == NULL,
- "erreur allocation memoire pour eval\na2ri_vf_oriented_bounding_box");
- gsl_matrix *evec = gsl_matrix_alloc (3, 3);
- a2ri_erreur_critique_si (evec == NULL,
- "erreur allocation memoire pour evec\na2ri_vf_oriented_bounding_box");
- gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (3);
- a2ri_erreur_critique_si (w == NULL,
- "erreur allocation memoire pour w\na2ri_vf_oriented_bounding_box");
- gsl_vector_view evec_i;
- point = gsl_matrix_alloc (1, 3);
- a2ri_erreur_critique_si (point == NULL,
- "erreur allocation memoire pour point\na2ri_vf_oriented_bounding_box");
- //calcul de la matrice de covariance
- listsommet = (point3d *) malloc (m->nbvertex * sizeof (point3d));
- a2ri_erreur_critique_si (listsommet == NULL,
- "erreur allocation memoire pour listsommet\na2ri_vf_oriented_bounding_box");
- for (int i = 0; i < m->nbvertex; i++)
- point3d_init (&(listsommet[i]), m->ve[i].x, m->ve[i].y, m->ve[i].z);
- covariance = cross_variance (listsommet, listsommet, m->nbvertex);
- //calcul des vecteurs propres qui seront les directions principales de la faces
- gsl_eigen_symmv (covariance, eval, evec, w);
- //récupération du premier vecteur propre
- evec_i = gsl_matrix_column (evec, 0);
- //affectation du nouveau vecteur "X" dans newbasex
- newbasex->dx = gsl_vector_get (&evec_i.vector, 0);
- newbasex->dy = gsl_vector_get (&evec_i.vector, 1);
- newbasex->dz = gsl_vector_get (&evec_i.vector, 2);
- data[0] = gsl_vector_get (&evec_i.vector, 0);
- data[3] = gsl_vector_get (&evec_i.vector, 1);
- data[6] = gsl_vector_get (&evec_i.vector, 2);
- //récupération du second vecteur propre
- evec_i = gsl_matrix_column (evec, 1);
- //affectation du nouveau vecteur "Y" dans newbasey
- newbasey->dx = gsl_vector_get (&evec_i.vector, 0);
- newbasey->dy = gsl_vector_get (&evec_i.vector, 1);
- newbasey->dz = gsl_vector_get (&evec_i.vector, 2);
- data[1] = gsl_vector_get (&evec_i.vector, 0);
- data[4] = gsl_vector_get (&evec_i.vector, 1);
- data[7] = gsl_vector_get (&evec_i.vector, 2);
- //récupération du troisième vecteur propre
- evec_i = gsl_matrix_column (evec, 2);
- //affectation du nouveau vecteur "Z" dans newbaseZ
- newbasez->dx = gsl_vector_get (&evec_i.vector, 0);
- newbasez->dy = gsl_vector_get (&evec_i.vector, 1);
- newbasez->dz = gsl_vector_get (&evec_i.vector, 2);
- data[2] = gsl_vector_get (&evec_i.vector, 0);
- data[5] = gsl_vector_get (&evec_i.vector, 1);
- data[8] = gsl_vector_get (&evec_i.vector, 2);
- //construction de la matrice de passage de l'ancienne à la nouvelle base
- matpass = matrix_init (data, 3, 3);
- //calcul des nouvelles coordonnées du premier sommet et mise à jour des coins inférieurs et supérieurs dans ptmin et ptmax
- gsl_matrix_set (point, 0, 0, m->ve[0].x);
- gsl_matrix_set (point, 0, 1, m->ve[0].y);
- gsl_matrix_set (point, 0, 2, m->ve[0].z);
- res = matrix_mul (point, matpass);
- ptmin->x = gsl_matrix_get (res, 0, 0);
- ptmax->x = gsl_matrix_get (res, 0, 0);
- ptmin->y = gsl_matrix_get (res, 0, 1);
- ptmax->y = gsl_matrix_get (res, 0, 1);
- ptmin->z = gsl_matrix_get (res, 0, 2);
- ptmax->z = gsl_matrix_get (res, 0, 2);
- //calcul des nouvelles coordonnées des autres sommets et mise à jour des coins inférieurs et supérieurs dans ptmin et ptmax
- for (int i = 1; i < m->nbvertex; i++)
- {
- gsl_matrix_set (point, 0, 0, m->ve[i].x);
- gsl_matrix_set (point, 0, 1, m->ve[i].y);
- gsl_matrix_set (point, 0, 2, m->ve[i].z);
- res = matrix_mul (point, matpass);
- if (ptmin->x > gsl_matrix_get (res, 0, 0))
- ptmin->x = gsl_matrix_get (res, 0, 0);
- if (ptmax->x < gsl_matrix_get (res, 0, 0))
- ptmax->x = gsl_matrix_get (res, 0, 0);
- if (ptmin->y > gsl_matrix_get (res, 0, 1))
- ptmin->y = gsl_matrix_get (res, 0, 1);
- if (ptmax->y < gsl_matrix_get (res, 0, 1))
- ptmax->y = gsl_matrix_get (res, 0, 1);
- if (ptmin->z > gsl_matrix_get (res, 0, 2))
- ptmin->z = gsl_matrix_get (res, 0, 2);
- if (ptmax->z < gsl_matrix_get (res, 0, 2))
- ptmax->z = gsl_matrix_get (res, 0, 2);
- }
- gsl_matrix_free (covariance);
- gsl_matrix_free (matpass);
- gsl_matrix_free (point);
- gsl_matrix_free (res);
- gsl_eigen_symmv_free (w);
- }
- /**
- Calcul des sphères englobantes des faces de la liste
- @param m le modèle
- @param faces tableau contenant les numéros des faces pour lesquels on veut les sphères englobantes
- @param nbfaces taille du tableau
- @param listcentre pointeur sur le tableau de point3d contenant les centres de sphères englobantes
- @param listradius pointeur sur le tableau contenant les rayons des sphères englobantes
- **/
- void
- a2ri_vef_bounding_ball_minimale (
- const vef_model * const m,
- const int * const faces,
- int nbfaces,
- point3d ** listcentre,
- double **listradius)
- {
- double radius,
- x,
- y,
- z;
- int ve1,
- ve2,
- ve3,
- nbelt = 0;
- point3d p3d_1,
- p3d_2,
- p3d_3,
- p3d_centre;
- double longar[3],
- longp1p2,
- longp1p3,
- longp2p3;
- free (*listcentre);
- *listcentre = (point3d *) malloc (nbfaces * sizeof (point3d));
- a2ri_erreur_critique_si (*listcentre == NULL,
- "erreur allocation memoire pour listcentre\na2ri_vef_bounding_ball_minimale");
- for (int alpha = 0; alpha < nbfaces; alpha++)
- {
- vef_face_get_vertices (&(m->fa[faces[alpha]]), m->ed, &ve1, &ve2, &ve3);
- point3d_init (&p3d_1, m->ve[ve1].x, m->ve[ve1].y, m->ve[ve1].z);
- point3d_init (&p3d_2, m->ve[ve2].x, m->ve[ve2].y, m->ve[ve2].z);
- point3d_init (&p3d_3, m->ve[ve3].x, m->ve[ve3].y, m->ve[ve3].z);
- longar[0] = point3d_square_length (&p3d_1, &p3d_2);
- longp1p2 = longar[0];
- longar[1] = point3d_square_length (&p3d_1, &p3d_3);
- longp1p3 = longar[1];
- longar[2] = point3d_square_length (&p3d_2, &p3d_3);
- longp2p3 = longar[2];
- list_double_sort (longar, 3, ASC);
- if (longar[2] > longar[0] + longar[1])
- {
- if (longar[2] == longp1p2)
- {
- (*listcentre)[alpha].x = ((p3d_1.x + p3d_2.x) * 1.0) / 2.0;
- (*listcentre)[alpha].y = ((p3d_1.y + p3d_2.y) * 1.0) / 2.0;
- (*listcentre)[alpha].z = ((p3d_1.z + p3d_2.z) * 1.0) / 2.0;
- }
- if (longar[2] == longp1p3)
- {
- (*listcentre)[alpha].x = ((p3d_1.x + p3d_3.x) * 1.0) / 2.0;
- (*listcentre)[alpha].y = ((p3d_1.y + p3d_3.y) * 1.0) / 2.0;
- (*listcentre)[alpha].z = ((p3d_1.z + p3d_3.z) * 1.0) / 2.0;
- }
- if (longar[2] == longp2p3)
- {
- (*listcentre)[alpha].x = ((p3d_2.x + p3d_3.x) * 1.0) / 2.0;
- (*listcentre)[alpha].y = ((p3d_2.y + p3d_3.y) * 1.0) / 2.0;
- (*listcentre)[alpha].z = ((p3d_2.z + p3d_3.z) * 1.0) / 2.0;
- }
- list_double_add (listradius, &nbelt, sqrt (longar[2]) / 2.0,
- WITH_REDUNDANCE);
- }
- else
- {
- p3d_centre = circumcircle_center (&p3d_1, &p3d_2, &p3d_3);
- x = p3d_centre.x;
- y = p3d_centre.y;
- z = p3d_centre.z;
- radius = sqrt ((m->ve[ve1].x - x) * (m->ve[ve1].x - x) +
- (m->ve[ve1].y - y) * (m->ve[ve1].y - y) +
- (m->ve[ve1].z - z) * (m->ve[ve1].z - z));
- (*listcentre)[alpha].x = x;
- (*listcentre)[alpha].y = y;
- (*listcentre)[alpha].z = z;
- list_double_add (listradius, &nbelt, radius, WITH_REDUNDANCE);
- }
- }
- }
- /**
- Calcul des sphères englobantes des faces de la liste
- @param m le modèle
- @param faces tableau contenant les numéros des faces pour lesquels on veut les sphères englobantes
- @param nbfaces taille du tableau
- @param listcentre pointeur sur le tableau de point3d contenant les centres de sphères englobantes
- @param listradius pointeur sur le tableau contenant les rayons des sphères englobantes
- **/
- void
- point3d_bounding_ball_minimale (
- const point3d * const faces,
- int nbfaces,
- point3d ** listcentre,
- double **listradius)
- {
- double radius,
- x,
- y,
- z;
- int nbelt = 0;
- point3d p3d_1,
- p3d_2,
- p3d_3,
- p3d_centre;
- double longar[3],
- longp1p2,
- longp1p3,
- longp2p3;
- free (*listcentre);
- *listcentre = (point3d *) malloc (nbfaces * sizeof (point3d));
- a2ri_erreur_critique_si (*listcentre == NULL,
- "erreur allocation memoire pour lsitcentre\nvef_bounding_ball_minimale");
- for (int alpha = 0; alpha < nbfaces; alpha++)
- {
- point3d_init (&p3d_1, faces[alpha * 3].x, faces[alpha * 3].y,
- faces[alpha * 3].z);
- point3d_init (&p3d_2, faces[alpha * 3 + 1].x, faces[alpha * 3 + 1].y,
- faces[alpha * 3 + 1].z);
- point3d_init (&p3d_3, faces[alpha * 3 + 2].x, faces[alpha * 3 + 2].y,
- faces[alpha * 3 + 2].z);
- longar[0] = point3d_square_length (&p3d_1, &p3d_2);
- longp1p2 = longar[0];
- longar[1] = point3d_square_length (&p3d_1, &p3d_3);
- longp1p3 = longar[1];
- longar[2] = point3d_square_length (&p3d_2, &p3d_3);
- longp2p3 = longar[2];
- list_double_sort (longar, 3, ASC);
- if (longar[2] > longar[0] + longar[1])
- {
- if (longar[2] == longp1p2)
- {
- (*listcentre)[alpha].x = ((p3d_1.x + p3d_2.x) * 1.0) / 2.0;
- (*listcentre)[alpha].y = ((p3d_1.y + p3d_2.y) * 1.0) / 2.0;
- (*listcentre)[alpha].z = ((p3d_1.z + p3d_2.z) * 1.0) / 2.0;
- }
- if (longar[2] == longp1p3)
- {
- (*listcentre)[alpha].x = ((p3d_1.x + p3d_3.x) * 1.0) / 2.0;
- (*listcentre)[alpha].y = ((p3d_1.y + p3d_3.y) * 1.0) / 2.0;
- (*listcentre)[alpha].z = ((p3d_1.z + p3d_3.z) * 1.0) / 2.0;
- }
- if (longar[2] == longp2p3)
- {
- (*listcentre)[alpha].x = ((p3d_2.x + p3d_3.x) * 1.0) / 2.0;
- (*listcentre)[alpha].y = ((p3d_2.y + p3d_3.y) * 1.0) / 2.0;
- (*listcentre)[alpha].z = ((p3d_2.z + p3d_3.z) * 1.0) / 2.0;
- }
- list_double_add (listradius, &nbelt, sqrt (longar[2]) / 2.0,
- WITH_REDUNDANCE);
- }
- else
- {
- p3d_centre = circumcircle_center (&p3d_1, &p3d_2, &p3d_3);
- x = p3d_centre.x;
- y = p3d_centre.y;
- z = p3d_centre.z;
- radius = sqrt ((p3d_1.x - x) * (p3d_1.x - x) +
- (p3d_1.y - y) * (p3d_1.y - y) +
- (p3d_1.z - z) * (p3d_1.z - z));
- (*listcentre)[alpha].x = x;
- (*listcentre)[alpha].y = y;
- (*listcentre)[alpha].z = z;
- list_double_add (listradius, &nbelt, radius, WITH_REDUNDANCE);
- }
- }
- }
|