1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555 |
- /*************************************/
- /* Auteur : Rémi Synave */
- /* Date de création : 03/04/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 "quality.h"
- /********** INTERMEDIATE TYPES AND FUNCTIONS **********/
- /* Les fonctions intermédiaires sont préfixées de IF */
- /* et les types intermédiaires de IT */
- #define NB_SAMPLE_HAUSDORFF 50
- #define NB_SAMPLE_EMN 50
- typedef struct
- {
- vf_model *m;
- double *list;
- int size;
- } ITcouple;
- //ajout de la longueur de l'arete value dans la liste contenu dans la variable user_data
- void
- IFlongueur_arete (
- int key,
- vf_edge * value,
- void * user_data)
- {
- ITcouple *c = user_data;
-
- list_double_add (&((*c).list), &((*c).size), value->att_double,
- WITH_REDUNDANCE);
- }
- /********** MAIN FUNCTIONS **********/
- /**
- Calcul de l'erreur mean ratio metric sur un modèle pour une face
- @param m le modèle
- @param numface numéro de la face à évaluer
- @return évaluation de la face
- */
- double
- a2ri_vf_mean_ratio_metric_for_a_face (
- const vf_model * const m,
- int numface)
- {
- int ve1 = m->fa[numface].ve1;
- int ve2 = m->fa[numface].ve2;
- int ve3 = m->fa[numface].ve3;
- point3d p[3];
- point2d pnvbase[3];
- vector3d AB,
- AC,
- U;
- double a,
- b,
- c,
- d;
- //initialisation des triangles
- point3d_init (p, m->ve[ve1].x, m->ve[ve1].y, m->ve[ve1].z);
- point3d_init (p + 1, m->ve[ve2].x, m->ve[ve2].y, m->ve[ve2].z);
- point3d_init (p + 2, m->ve[ve3].x, m->ve[ve3].y, m->ve[ve3].z);
- //recherche equation du plan
- equation_plan (p, 3, &a, &b, &c, &d);
- //changement de base
- //initialisation du premier vecteur de la nouvelle base
- vector3d_init (&AB, p[1].x - p[0].x, p[1].y - p[0].y, p[1].z - p[0].z);
- vector3d_init (&AC, p[2].x - p[0].x, p[2].y - p[0].y, p[2].z - p[0].z);
- //recherche du second vecteur de la nouvelle base
- find_second_base_vector (&AB, &AC, &U);
- //normalisation des vecteurs
- vector3d_normalize (&AB);
- vector3d_normalize (&U);
- //recuperation des nouvelles coordonnées du point dans la nouvelle base
- base_modification_3d_to_2d (&(p[0]), &(p[0]), &AB, &U, &pnvbase[0]);
- base_modification_3d_to_2d (&(p[1]), &(p[0]), &AB, &U, &pnvbase[1]);
- base_modification_3d_to_2d (&(p[2]), &(p[0]), &AB, &U, &pnvbase[2]);
- //application du calcul mean ratio metric
- //-> article : A comparison of two optimization methods for mesh quality improvement
- gsl_matrix *A = gsl_matrix_alloc (2, 2);
- a2ri_erreur_critique_si (A == NULL,
- "erreur allocation memoire pour A\na2ri_vf_mean_ratio_metric_for_a_face");
- gsl_matrix_set (A, 0, 0, pnvbase[1].x - pnvbase[0].x);
- gsl_matrix_set (A, 0, 1, pnvbase[2].x - pnvbase[0].x);
- gsl_matrix_set (A, 1, 0, pnvbase[1].y - pnvbase[0].y);
- gsl_matrix_set (A, 1, 1, pnvbase[2].y - pnvbase[0].y);
- gsl_matrix *Winv = gsl_matrix_alloc (2, 2);
- a2ri_erreur_critique_si (Winv == NULL,
- "erreur allocation memoire pour Winv\na2ri_vf_mean_ratio_metric_for_a_face");
- gsl_matrix_set (Winv, 0, 0, 1);
- gsl_matrix_set (Winv, 0, 1, -1.0 / sqrt (3.0));
- gsl_matrix_set (Winv, 1, 0, 0);
- gsl_matrix_set (Winv, 1, 1, 2.0 / sqrt (3.0));
- gsl_matrix *S = matrix_mul (A, Winv);
- double erreurface =
- ((2 * matrix_determinant (S)) /
- (matrix_frobenius_norm (S) * matrix_frobenius_norm (S)));
- if (erreurface < 0)
- erreurface *= -1;
- gsl_matrix_free (A);
- gsl_matrix_free (Winv);
- gsl_matrix_free (S);
- return erreurface;
- }
- /**
- Calcul de l'erreur mean ratio metric sur un modèle complet
- @param m le modèle
- @return évaluation du modèle
- */
- double
- a2ri_vf_mean_ratio_metric (
- const vf_model * const m)
- {
- double somme = 0.0;
- for (int i = 0; i < m->nbface; i++)
- somme += a2ri_vf_mean_ratio_metric_for_a_face (m, i);
- somme /= (m->nbface) * 1.0;
- return somme;
- }
- /**
- Calcule la liste des angles des triangles composant le modèle
- @param m le modèle
- @param list pointeur sur le tableau contenant la liste des angles
- @param size pointeur sur la taille du tableau
- @return aucun
- */
- void
- a2ri_vf_list_angle (
- const vf_model * const m,
- double ** list,
- int * size)
- {
- *size = 0;
- *list = NULL;
- int ve1,
- ve2,
- ve3;
- double dx,
- dy,
- dz,
- angledegre;
- vector3d v1,
- v2;
- for (int i = 0; i < m->nbface; i++)
- {
- ve1 = m->fa[i].ve1;
- ve2 = m->fa[i].ve2;
- ve3 = m->fa[i].ve3;
- //angle ve1-ve2 ve1-ve3
- dx = m->ve[ve2].x - m->ve[ve1].x;
- dy = m->ve[ve2].y - m->ve[ve1].y;
- dz = m->ve[ve2].z - m->ve[ve1].z;
- vector3d_init (&v1, dx, dy, dz);
- dx = m->ve[ve3].x - m->ve[ve1].x;
- dy = m->ve[ve3].y - m->ve[ve1].y;
- dz = m->ve[ve3].z - m->ve[ve1].z;
- vector3d_init (&v2, dx, dy, dz);
- angledegre = vector3d_angle_degre (&v1, &v2);
- list_double_add (list, size, angledegre, WITH_REDUNDANCE);
- //angle ve2-ve1 ve2-ve3
- dx = m->ve[ve1].x - m->ve[ve2].x;
- dy = m->ve[ve1].y - m->ve[ve2].y;
- dz = m->ve[ve1].z - m->ve[ve2].z;
- vector3d_init (&v1, dx, dy, dz);
- dx = m->ve[ve3].x - m->ve[ve2].x;
- dy = m->ve[ve3].y - m->ve[ve2].y;
- dz = m->ve[ve3].z - m->ve[ve2].z;
- vector3d_init (&v2, dx, dy, dz);
- angledegre = vector3d_angle_degre (&v1, &v2);
- list_double_add (list, size, angledegre, WITH_REDUNDANCE);
- //angle ve3-ve2 ve3-ve1
- dx = m->ve[ve2].x - m->ve[ve3].x;
- dy = m->ve[ve2].y - m->ve[ve3].y;
- dz = m->ve[ve2].z - m->ve[ve3].z;
- vector3d_init (&v1, dx, dy, dz);
- dx = m->ve[ve1].x - m->ve[ve3].x;
- dy = m->ve[ve1].y - m->ve[ve3].y;
- dz = m->ve[ve1].z - m->ve[ve3].z;
- vector3d_init (&v2, dx, dy, dz);
- angledegre = vector3d_angle_degre (&v1, &v2);
- list_double_add (list, size, angledegre, WITH_REDUNDANCE);
- }
- }
- /**
- Calcule la liste des aires des triangles composant le modèle
- @param m le modèle
- @param list pointeur sur le tableau contenant la liste des aires
- @param size pointeur sur la taille du tableau
- @return aucun
- */
- void
- a2ri_vf_list_area (
- const vf_model * const m,
- double ** list,
- int * size)
- {
- int ve1,
- ve2,
- ve3;
- point3d p1,
- p2,
- p3;
- *list = NULL;
- *size = 0;
- double aire;
- for (int i = 0; i < m->nbface; i++)
- {
- ve1 = m->fa[i].ve1;
- ve2 = m->fa[i].ve2;
- ve3 = m->fa[i].ve3;
- point3d_init (&p1, m->ve[ve1].x, m->ve[ve1].y, m->ve[ve1].z);
- point3d_init (&p2, m->ve[ve2].x, m->ve[ve2].y, m->ve[ve2].z);
- point3d_init (&p3, m->ve[ve3].x, m->ve[ve3].y, m->ve[ve3].z);
- aire = point3d_area (&p1, &p2, &p3);
- list_double_add (list, size, aire, WITH_REDUNDANCE);
- }
- }
- /**
- Calcule la liste des valences des sommets composant le modèle
- @param m le modèle
- @param list pointeur sur le tableau contenant la liste des valences
- @param size pointeur sur la taille du tableau
- @return aucun
- */
- void
- a2ri_vf_list_valence (
- const vf_model * const m,
- int ** list,
- int * size)
- {
- *list = NULL;
- *size = 0;
- for (int i = 0; i < m->nbvertex; i++)
- list_int_add (list, size, m->ve[i].nbincidentvertices, WITH_REDUNDANCE);
- }
- /**
- Calcule la liste des longueurs des aretes composant le modèle
- @param m le modèle
- @param list pointeur sur le tableau contenant la liste des longueurs des aretes
- @param size pointeur sur la taille du tableau
- @return aucun
- */
- //On ne peut pas garantir le const du vf_model ici
- void
- a2ri_vf_list_edge_length (
- vf_model * m,
- double ** list,
- int * size)
- {
- ITcouple c;
- c.m = m;
- c.list = NULL;
- c.size = 0;
- ptf_func_hashtable func[1];
- func[0] = a2ri_vf_update_length_edge;
- //construction de la table de hachage des aretes
- hashtable *table = a2ri_vf_construction_edge_table (m, func, 1);
- //pour toute les aretes, on ajoute leur longueur dans la liste
- hashtable_foreach (table, IFlongueur_arete, &c);
- *list = c.list;
- *size = c.size;
- hashtable_free (table);
- free (table);
- }
- /**
- Calcule la liste des longueurs des hauteurs des triangles composant le modèle
- @param m le modèle
- @param list pointeur sur le tableau contenant la liste des hauteurs
- @param size pointeur sur la taille du tableau
- @return aucun
- */
- //On ne peut pas garantir le const du vf_model ici
- void
- a2ri_vf_list_height_length (
- vf_model * m,
- double ** list,
- int * size)
- {
- *list = NULL;
- *size = 0;
- point3d p1,
- p2,
- p3;
- vector3d v1,
- v2,
- vp;
- for (int i = 0; i < m->nbface; i++)
- {
- //initialisation des trois points correpondant aux trois sommets de la face
- point3d_init (&p1, m->ve[m->fa[i].ve1].x, m->ve[m->fa[i].ve1].y,
- m->ve[m->fa[i].ve1].z);
- point3d_init (&p2, m->ve[m->fa[i].ve2].x, m->ve[m->fa[i].ve2].y,
- m->ve[m->fa[i].ve2].z);
- point3d_init (&p3, m->ve[m->fa[i].ve3].x, m->ve[m->fa[i].ve3].y,
- m->ve[m->fa[i].ve3].z);
- //ajout de la première hauteur du triangle
- vector3d_init (&v1, p3.x - p2.x, p3.y - p2.y, p3.z - p2.z);
- vector3d_init (&v2, p2.x - p1.x, p2.y - p1.y, p2.z - p1.z);
- vp = vector3d_vectorialproduct (&v1, &v2);
- list_double_add (list, size, vector3d_size (&vp) / vector3d_size (&v1),
- WITH_REDUNDANCE);
- //ajout de la seconde hauteur du triangle
- vector3d_init (&v1, p1.x - p3.x, p1.y - p3.y, p1.z - p3.z);
- vector3d_init (&v2, p3.x - p2.x, p3.y - p2.y, p3.z - p2.z);
- vp = vector3d_vectorialproduct (&v1, &v2);
- list_double_add (list, size, vector3d_size (&vp) / vector3d_size (&v1),
- WITH_REDUNDANCE);
- //ajout de la troisieme hauteur du triangle
- vector3d_init (&v1, p2.x - p1.x, p2.y - p1.y, p2.z - p1.z);
- vector3d_init (&v2, p1.x - p3.x, p1.y - p3.y, p1.z - p3.z);
- vp = vector3d_vectorialproduct (&v1, &v2);
- list_double_add (list, size, vector3d_size (&vp) / vector3d_size (&v1),
- WITH_REDUNDANCE);
- }
- }
- /**
- Compte le nombre d'angles du modèle par intervalle de 10 degres
- @param m le modèle
- @return un tableau de 18 entiers contenant le nombre d'angle compris entre 0-10, 10-20, ... , 170-180
- */
- int *
- a2ri_vf_angle_measure (
- const vf_model * const m)
- {
- //|0-10|10-20|20-30|...|160-170|170-180|
- int *angle = (int *) malloc (18 * sizeof (int));
- a2ri_erreur_critique_si (angle == NULL,
- "erreur allocation memoire pour angle\na2ri_vf_angle_measure");
- int ve1,
- ve2,
- ve3,
- tabindice;
- double dx,
- dy,
- dz,
- angledegre;
- vector3d v1,
- v2;
- for (int i = 0; i < 18; i++)
- angle[i] = 0;
- for (int i = 0; i < m->nbface; i++)
- {
- ve1 = m->fa[i].ve1;
- ve2 = m->fa[i].ve2;
- ve3 = m->fa[i].ve3;
- //angle ve1-ve2 ve1-ve3
- dx = m->ve[ve2].x - m->ve[ve1].x;
- dy = m->ve[ve2].y - m->ve[ve1].y;
- dz = m->ve[ve2].z - m->ve[ve1].z;
- vector3d_init (&v1, dx, dy, dz);
- dx = m->ve[ve3].x - m->ve[ve1].x;
- dy = m->ve[ve3].y - m->ve[ve1].y;
- dz = m->ve[ve3].z - m->ve[ve1].z;
- vector3d_init (&v2, dx, dy, dz);
- angledegre = vector3d_angle_degre (&v1, &v2);
- angledegre *= 10;
- tabindice = (int) angledegre;
- if (tabindice % 10 >= 5)
- tabindice += 5;
- tabindice /= 100;
- angle[tabindice] = angle[tabindice] + 1;
- //angle ve2-ve1 ve2-ve3
- dx = m->ve[ve1].x - m->ve[ve2].x;
- dy = m->ve[ve1].y - m->ve[ve2].y;
- dz = m->ve[ve1].z - m->ve[ve2].z;
- vector3d_init (&v1, dx, dy, dz);
- dx = m->ve[ve3].x - m->ve[ve2].x;
- dy = m->ve[ve3].y - m->ve[ve2].y;
- dz = m->ve[ve3].z - m->ve[ve2].z;
- vector3d_init (&v2, dx, dy, dz);
- angledegre = vector3d_angle_degre (&v1, &v2);
- angledegre *= 10;
- tabindice = (int) angledegre;
- if (tabindice % 10 >= 5)
- tabindice += 5;
- tabindice /= 100;
- angle[tabindice] = angle[tabindice] + 1;
- //angle ve3-ve2 ve3-ve1
- dx = m->ve[ve2].x - m->ve[ve3].x;
- dy = m->ve[ve2].y - m->ve[ve3].y;
- dz = m->ve[ve2].z - m->ve[ve3].z;
- vector3d_init (&v1, dx, dy, dz);
- dx = m->ve[ve1].x - m->ve[ve3].x;
- dy = m->ve[ve1].y - m->ve[ve3].y;
- dz = m->ve[ve1].z - m->ve[ve3].z;
- vector3d_init (&v2, dx, dy, dz);
- angledegre = vector3d_angle_degre (&v1, &v2);
- angledegre *= 10;
- tabindice = (int) angledegre;
- if (tabindice % 10 >= 5)
- tabindice += 5;
- tabindice /= 100;
- angle[tabindice] = angle[tabindice] + 1;
- }
- return angle;
- }
- /**
- Calcul de l'erreur de Hausdorff
- @param m1 premier modèle
- @param m2 second modèle
- @return erreur de Hausdorff
- */
- double
- a2ri_vf_hausdorff (
- const vf_model * const m1,
- const vf_model * const m2,
- int sampling)
- {
- point3d *liste = NULL;
- point3d *listesamplem1 = NULL;
- point3d *listesamplem2 = NULL;
- double max = 0.0;
- point3d p,
- A,
- B,
- C;
- //premier cas : pas d'echantillonnage
- if (sampling == NO_SAMPLING)
- {
- //transformation des sommets de m2 en lsite de points
- liste = (point3d *) malloc (m2->nbvertex * sizeof (point3d));
- a2ri_erreur_critique_si (liste == NULL,
- "erreur allocation memoire pour liste\na2ri_vf_hausdorff");
- for (int i = 0; i < m2->nbvertex; i++)
- {
- liste[i].x = m2->ve[i].x;
- liste[i].y = m2->ve[i].y;
- liste[i].z = m2->ve[i].z;
- }
- //on cherche la longueur maximale entre les sommets de m1 et les somemts de m2
- for (int i = 0; i < m1->nbvertex; i++)
- {
- point3d_init (&p, m1->ve[i].x, m1->ve[i].y, m1->ve[i].z);
- double temp = DPP (&p, liste, m2->nbvertex);
- if (max < temp)
- max = temp;
- }
- free (liste);
- liste = NULL;
- //transformation des sommets de m1 en lsite de points
- liste = (point3d *) malloc (m1->nbvertex * sizeof (point3d));
- a2ri_erreur_critique_si (liste == NULL,
- "erreur allocation memoire pour liste\na2ri_vf_hausdorff");
- for (int i = 0; i < m1->nbvertex; i++)
- {
- liste[i].x = m1->ve[i].x;
- liste[i].y = m1->ve[i].y;
- liste[i].z = m1->ve[i].z;
- }
- //on cherche la longueur maximale entre les sommets de m2 et les somemts de m1
- for (int i = 0; i < m2->nbvertex; i++)
- {
- point3d_init (&p, m2->ve[i].x, m2->ve[i].y, m2->ve[i].z);
- double temp = DPP (&p, liste, m1->nbvertex);
- if (max < temp)
- max = temp;
- }
-
- free (liste);
- }
- else
- {
- listesamplem1 =
- (point3d *) malloc ((m1->nbvertex + m1->nbface * NB_SAMPLE_HAUSDORFF)
- * sizeof (point3d));
- a2ri_erreur_critique_si (listesamplem1 == NULL,
- "erreur allocation memoire pour listesamplem1\na2ri_vf_hausdorff");
- listesamplem2 =
- (point3d *) malloc ((m2->nbvertex + m2->nbface * NB_SAMPLE_HAUSDORFF)
- * sizeof (point3d));
- a2ri_erreur_critique_si (listesamplem2 == NULL,
- "erreur allocation memoire pour listesamplem2\na2ri_vf_hausdorff");
- //echantillonage de m1 et m2
- for (int i = 0; i < m1->nbvertex; i++)
- {
- listesamplem1[i].x = m1->ve[i].x;
- listesamplem1[i].y = m1->ve[i].y;
- listesamplem1[i].z = m1->ve[i].z;
- }
- for (int i = 0; i < m2->nbvertex; i++)
- {
- listesamplem2[i].x = m2->ve[i].x;
- listesamplem2[i].y = m2->ve[i].y;
- listesamplem2[i].z = m2->ve[i].z;
- }
- for (int i = 0; i < m1->nbface; i++)
- {
- point3d_init (&A, m1->ve[m1->fa[i].ve1].x, m1->ve[m1->fa[i].ve1].y,
- m1->ve[m1->fa[i].ve1].z);
- point3d_init (&B, m1->ve[m1->fa[i].ve2].x, m1->ve[m1->fa[i].ve2].y,
- m1->ve[m1->fa[i].ve2].z);
- point3d_init (&C, m1->ve[m1->fa[i].ve3].x, m1->ve[m1->fa[i].ve3].y,
- m1->ve[m1->fa[i].ve3].z);
- sample_triangle (&A, &B, &C, NB_SAMPLE_HAUSDORFF, &liste);
- for (int j = 0; j < NB_SAMPLE_HAUSDORFF; j++)
- {
- listesamplem1[(i * NB_SAMPLE_HAUSDORFF) + m1->nbvertex + j].x =
- liste[i].x;
- listesamplem1[(i * NB_SAMPLE_HAUSDORFF) + m1->nbvertex + j].y =
- liste[i].y;
- listesamplem1[(i * NB_SAMPLE_HAUSDORFF) + m1->nbvertex + j].z =
- liste[i].z;
- }
- }
- for (int i = 0; i < m2->nbface; i++)
- {
- point3d_init (&A, m2->ve[m2->fa[i].ve1].x, m2->ve[m2->fa[i].ve1].y,
- m2->ve[m2->fa[i].ve1].z);
- point3d_init (&B, m2->ve[m2->fa[i].ve2].x, m2->ve[m2->fa[i].ve2].y,
- m2->ve[m2->fa[i].ve2].z);
- point3d_init (&C, m2->ve[m2->fa[i].ve3].x, m2->ve[m2->fa[i].ve3].y,
- m2->ve[m2->fa[i].ve3].z);
- sample_triangle (&A, &B, &C, NB_SAMPLE_HAUSDORFF, &liste);
- for (int j = 0; j < NB_SAMPLE_HAUSDORFF; j++)
- {
- listesamplem2[(i * NB_SAMPLE_HAUSDORFF) + m2->nbvertex + j].x =
- liste[i].x;
- listesamplem2[(i * NB_SAMPLE_HAUSDORFF) + m2->nbvertex + j].y =
- liste[i].y;
- listesamplem2[(i * NB_SAMPLE_HAUSDORFF) + m2->nbvertex + j].z =
- liste[i].z;
- }
- }
- //recherche de la longueur max entre m1 et m2
- for (int i = 0; i < m1->nbvertex + m1->nbface * NB_SAMPLE_HAUSDORFF;
- i++)
- {
- double temp =
- DPP (&(listesamplem1[i]), listesamplem2,
- m2->nbvertex + m2->nbface * NB_SAMPLE_HAUSDORFF);
- if (max < temp)
- max = temp;
- }
- //recherche de la distance max entre m2 et m1
- for (int i = 0; i < m2->nbvertex + m2->nbface * NB_SAMPLE_HAUSDORFF;
- i++)
- {
- double temp =
- DPP (&(listesamplem2[i]), listesamplem1,
- m1->nbvertex + m1->nbface * NB_SAMPLE_HAUSDORFF);
- if (max < temp)
- max = temp;
- }
- free (listesamplem1);
- free (listesamplem2);
- }
- return max;
- }
- /**
- Calcul de la distance entre un point et un ensemble de points en se basant sur la métrique Point-Point
- @param p le point
- @param points tableau de réel contenant les coordoonées des sommets du maillage
- @param size taille du tableau
- @return distance entre le point et le modèle
- */
- double
- DPP (
- const point3d * const p,
- const point3d *const points,
- int size)
- {
- double min = 0;
- //on cherche la distance minimale entre un point p et un ensemble de points
- for (int i = 0; i < size; i++)
- {
- double temp = point3d_length (p, &(points[i]));
- if (i == 0)
- min = temp;
- else if (temp < min)
- min = temp;
- }
- return min;
- }
- /**
- Calcul de l'erreur moyenne "Emean" entre deux modèles
- @param m1 le premier modèle
- @param m2 le second modèle
- @param sampling SAMPLING ou NO_SAMPLING <BR> indique s'il faut ou non échantillonner les faces.
- @return l'erreur "Emean"
- */
- double
- a2ri_vf_Emn_DPP (
- const vf_model * const m1,
- const vf_model * const m2,
- int sampling)
- {
- point3d *liste = NULL;
- double erreur = 0.0;
- point3d p;
- //premier cas : pas d'echantillonnage
- if (sampling == NO_SAMPLING)
- {
- //transformation de la liste de sommets de m2 en liste de points
- liste = (point3d *) malloc (m2->nbvertex * sizeof (point3d));
- a2ri_erreur_critique_si (liste == NULL,
- "erreur allocation memoire pour liste\na2ri_vf_Emn_DPP");
- for (int i = 0; i < m2->nbvertex; i++)
- {
- liste[i].x = m2->ve[i].x;
- liste[i].y = m2->ve[i].y;
- liste[i].z = m2->ve[i].z;
- }
- //somme DPP(m1,m2)*DPP(m1,m2)
- for (int i = 0; i < m1->nbvertex; i++)
- {
- point3d_init (&p, m1->ve[i].x, m1->ve[i].y, m1->ve[i].z);
- double temp = DPP (&p, liste, m2->nbvertex);
- erreur += temp * temp;
- }
- free (liste);
- liste = NULL;
- //transformation de la liste de sommets de m1 en liste de points
- liste = (point3d *) malloc (m1->nbvertex * sizeof (point3d));
- a2ri_erreur_critique_si (liste == NULL,
- "erreur allocation memoire pour liste\na2ri_vf_Emn_DPP");
- for (int i = 0; i < m1->nbvertex; i++)
- {
- liste[i].x = m1->ve[i].x;
- liste[i].y = m1->ve[i].y;
- liste[i].z = m1->ve[i].z;
- }
- //somme DPP(m2,m1)*DPP(m2,m1)
- for (int i = 0; i < m2->nbvertex; i++)
- {
- point3d_init (&p, m2->ve[i].x, m2->ve[i].y, m2->ve[i].z);
- double temp = DPP (&p, liste, m1->nbvertex);
- erreur += temp * temp;
- }
- free (liste);
- }
- else
- {
- /*TODO meme algo avec echantillonnage*/
- erreur = -1;
- }
- return erreur / ((m1->nbvertex + m2->nbvertex) * 1.0);
- }
- /**
- Evaluation de l'aspect ratio d'un triangle
- @param m le modèle
- @param numface numéro de la face à évaluer
- @return évaluation de l'aspect ratio
- */
- double
- a2ri_vf_triangle_aspect_ratio (
- const vf_model * const m,
- int numface)
- {
- double distmax = 0,
- hauteur;
- int index = 0;
- point3d p1,
- p2,
- p3;
- vector3d v1,
- v2,
- vp;
- point3d_init (&p1, m->ve[m->fa[numface].ve1].x, m->ve[m->fa[numface].ve1].y,
- m->ve[m->fa[numface].ve1].z);
- point3d_init (&p2, m->ve[m->fa[numface].ve2].x, m->ve[m->fa[numface].ve2].y,
- m->ve[m->fa[numface].ve2].z);
- point3d_init (&p3, m->ve[m->fa[numface].ve3].x, m->ve[m->fa[numface].ve3].y,
- m->ve[m->fa[numface].ve3].z);
- //on cherche le plus grand cote
- distmax = point3d_length (&p1, &p2);
- if (distmax < point3d_length (&p1, &p3))
- {
- index = 1;
- distmax = point3d_length (&p1, &p3);
- }
- if (distmax < point3d_length (&p2, &p3))
- {
- index = 2;
- distmax = point3d_length (&p2, &p3);
- }
- //on cherche la hauteur du coté le plus grand
- switch (index)
- {
- case 0:
- vector3d_init (&v1, p1.x - p2.x, p1.y - p2.y, p1.z - p2.z);
- vector3d_init (&v2, p3.x - p2.x, p3.y - p2.y, p3.z - p2.z);
- break;
- case 1:
- vector3d_init (&v1, p3.x - p1.x, p3.y - p1.y, p3.z - p1.z);
- vector3d_init (&v2, p2.x - p1.x, p2.y - p1.y, p2.z - p1.z);
- break;
- case 2:
- vector3d_init (&v1, p2.x - p3.x, p2.y - p3.y, p2.z - p3.z);
- vector3d_init (&v2, p1.x - p3.x, p1.y - p3.y, p1.z - p3.z);
- break;
- }
- vp = vector3d_vectorialproduct (&v1, &v2);
- hauteur = vector3d_size (&vp) / vector3d_size (&v1);
- return (vector3d_size (&v1) / hauteur);
- }
- /**
- Evaluation de la qualité d'un triangle avec la méthode de Rypl
- @param m le modèle
- @param numface numéro de la face à évaluer
- @return évaluation de l'aspect ratio
- */
- double
- a2ri_vf_triangle_rypl (
- const vf_model * const m,
- int numface)
- {
- double a,
- b,
- c,
- aire;
- point3d p1,
- p2,
- p3;
- point3d_init (&p1, m->ve[m->fa[numface].ve1].x, m->ve[m->fa[numface].ve1].y,
- m->ve[m->fa[numface].ve1].z);
- point3d_init (&p2, m->ve[m->fa[numface].ve2].x, m->ve[m->fa[numface].ve2].y,
- m->ve[m->fa[numface].ve2].z);
- point3d_init (&p3, m->ve[m->fa[numface].ve3].x, m->ve[m->fa[numface].ve3].y,
- m->ve[m->fa[numface].ve3].z);
- aire = point3d_area (&p1, &p2, &p3);
- a = point3d_length (&p2, &p3);
- a = a * a;
- b = point3d_length (&p1, &p3);
- b = b * b;
- c = point3d_length (&p1, &p2);
- c = c * c;
- return (4 * sqrt (3) * aire / (a + b + c));
- }
- /**
- Calcul de la liste d'aspect ratio d'un modèle
- @param m le modèle
- @param list pointeur sur un tableau contenant tous les aspect ratio des triangles composant le modèle
- @param siuze pointeur sur la taille du tableau
- @return aucun
- */
- void
- a2ri_vf_list_triangle_aspect_ratio (
- const vf_model * const m,
- double ** list,
- int * size)
- {
- *list = NULL;
- *size = 0;
- for (int i = 0; i < m->nbface; i++)
- list_double_add (list, size, a2ri_vf_triangle_aspect_ratio (m, i),
- WITH_REDUNDANCE);
- }
- /**
- Calcul de la variation de taille (d'aire) entre deux faces du modèle
- @param m le modèle
- @param face1 numéro de la première face servant de base
- @param face2 numéro de la seconde face
- @return variation de la taille : Aire_face1/Aire_face2
- */
- double
- a2ri_vf_face_size_variation (
- const vf_model * const m,
- int face1,
- int face2)
- {
- point3d A,
- B,
- C,
- D,
- E,
- F;
- point3d_init (&A, m->ve[m->fa[face1].ve1].x, m->ve[m->fa[face1].ve1].y,
- m->ve[m->fa[face1].ve1].z);
- point3d_init (&B, m->ve[m->fa[face1].ve2].x, m->ve[m->fa[face1].ve2].y,
- m->ve[m->fa[face1].ve2].z);
- point3d_init (&C, m->ve[m->fa[face1].ve3].x, m->ve[m->fa[face1].ve3].y,
- m->ve[m->fa[face1].ve3].z);
- point3d_init (&D, m->ve[m->fa[face2].ve1].x, m->ve[m->fa[face2].ve1].y,
- m->ve[m->fa[face2].ve1].z);
- point3d_init (&E, m->ve[m->fa[face2].ve2].x, m->ve[m->fa[face2].ve2].y,
- m->ve[m->fa[face2].ve2].z);
- point3d_init (&F, m->ve[m->fa[face2].ve3].x, m->ve[m->fa[face2].ve3].y,
- m->ve[m->fa[face2].ve3].z);
- return (point3d_area (&A, &B, &C) / point3d_area (&D, &E, &F));
- }
- /**
- Calcul de la liste des rayons des cercles inscrit aux faces du modèles
- @param m le modèle
- @param list pointeur sur le tableau contenant la liste des rayons des cercles inscrits
- @param size pointeur sur la taille du tableau
- @return aucun
- */
- void
- a2ri_vf_list_radius_incircle (
- const vf_model * const m,
- double ** list,
- int * size)
- {
- point3d A,
- B,
- C,
- center;
- for (int i = 0; i < m->nbface; i++)
- {
- point3d_init (&A, m->ve[m->fa[i].ve1].x, m->ve[m->fa[i].ve1].y,
- m->ve[m->fa[i].ve1].z);
- point3d_init (&B, m->ve[m->fa[i].ve2].x, m->ve[m->fa[i].ve2].y,
- m->ve[m->fa[i].ve2].z);
- point3d_init (&C, m->ve[m->fa[i].ve3].x, m->ve[m->fa[i].ve3].y,
- m->ve[m->fa[i].ve3].z);
- center = incircle_center (&A, &B, &C);
- list_double_add (list, size,
- distance_point_straight_line (¢er, &A, &B),
- WITH_REDUNDANCE);
- }
- }
- /**
- Calcul de la courbure de gauss pour un sommet
- @param m le modele
- @param ve numéro du sommet
- @return le courbure de gauss en ce sommet
- */
- double
- a2ri_vf_gauss_curvature (
- const vf_model * const m,
- int ve)
- {
- //TODO vérifier le bon fonctionnement de la fonction
- vector3d v1,
- v2;
- int *star = NULL,
- size = 0,
- ve1,
- ve2,
- ve3;
- double sommeangle = 0;
- a2ri_vf_faces_next_vertex (m, ve, &star, &size);
- for (int i = 0; i < size; i++)
- {
- ve1 = m->fa[star[i]].ve1;
- ve2 = m->fa[star[i]].ve2;
- ve3 = m->fa[star[i]].ve3;
- if (ve == ve2)
- {
- ve2 = ve1;
- ve1 = ve;
- }
- if (ve == ve3)
- {
- ve3 = ve1;
- ve1 = ve;
- }
- vector3d_init (&v1, m->ve[ve2].x - m->ve[ve1].x,
- m->ve[ve2].y - m->ve[ve1].y,
- m->ve[ve2].z - m->ve[ve1].z);
- vector3d_init (&v2, m->ve[ve3].x - m->ve[ve1].x,
- m->ve[ve3].y - m->ve[ve1].y,
- m->ve[ve3].z - m->ve[ve1].z);
- sommeangle += vector3d_angle_radian (&v1, &v2);
- }
- free (star);
- return (2.0 * PI - sommeangle);
- }
- /**
- Calcul de la courbure de gauss pour une liste de sommet
- @param m le modele
- @param ve liste de sommets
- @param size taille de la liste
- @param list liste des courbures associés à chaque sommet
- @return aucun
- */
- void
- a2ri_vf_gauss_curvature_list_of_vertex (
- const vf_model * const m,
- const int * const ve,
- int size,
- double ** list)
- {
- //TODO vérifier le bon fonctionnement de la fonction
- int newsize = 0;
- *list = NULL;
- for (int i = 0; i < size; i++)
- list_double_add (list, &newsize, a2ri_vf_gauss_curvature (m, ve[i]),
- WITH_REDUNDANCE);
- }
- /**
- calcule de la mean curvature pour une arete
- @param m le modele
- @param ve1 premier sommet de l'arete
- @param ve2 second sommet de l'arete
- @return la mean curvature pour cette arete
- */
- vector3d
- a2ri_vf_mean_curvature (
- const vf_model * const m,
- int ve1,
- int ve2)
- {
- int *list1 = NULL,
- *list2 = NULL,
- *list3 = NULL;
- int size1 = 0,
- size2 = 0,
- size3 = 0;
- vector3d retour,
- e,
- mu1,
- mu2;
- a2ri_vf_faces_next_vertex (m, ve1, &list1, &size1);
- a2ri_vf_faces_next_vertex (m, ve2, &list2, &size2);
- list_int_intersection (list1, size1, list2, size2, &list3, &size3);
- if (size3 != 2)
- {
- vector3d_init (&retour, 0, 0, 0);
- return retour;
- }
- vector3d_init (&e, m->ve[ve1].x - m->ve[ve2].x, m->ve[ve1].y - m->ve[ve2].y,
- m->ve[ve1].z - m->ve[ve2].z);
- mu1 = a2ri_vf_normal_face (m, list3[0]);
- mu2 = a2ri_vf_normal_face (m, list3[1]);
- vector3d_normalize (&mu1);
- vector3d_normalize (&mu2);
- vector3d vp1=vector3d_vectorialproduct (&e, &mu1);
- vector3d vp2=vector3d_vectorialproduct (&e, &mu2);
- return (vector3d_sub(&vp1,&vp2));
- }
- /**
- Calcul de la mean curvature pour un chemin
- @param m le modele
- @param list chemin décrit par ses sommets
- @param size taille du tableau list
- @param list_vector liste des vector3d retourner pour chaque arete
- @param size_list_vector taille de la liste list_vector
- @return aucun
- */
- void
- a2ri_vf_mean_curvature_path (
- const vf_model * const m,
- const int * const list,
- int size,
- vector3d ** list_vector,
- int * size_list_vector)
- {
- *size_list_vector = 0;
- *list_vector = NULL;
- for (int i = 0; i < size - 1; i++)
- {
- vector3d vtemp=a2ri_vf_mean_curvature (m, list[i], list[i + 1]);
- list_vector3d_add (list_vector, size_list_vector,
- &vtemp,
- WITH_REDUNDANCE);
- }
- }
- /**
- Calcul de la courbure de Levi-Civita
- @param m le modele
- @param numve numéro du sommet
- @return la courbure de Levi-Civita
- */
- double
- a2ri_vf_levi_civita (
- const vf_model * const m,
- int numve)
- {
- double sommeangle = 0;
- vector3d v1,
- v2;
- int index,
- fa1,
- fa2,
- *dejavu = NULL,
- size = 0;
- if (m->ve[numve].nbincidentvertices == 0)
- return 0;
- int prem = m->ve[numve].incidentvertices[0];
- index = prem++;
- list_int_add (&dejavu, &size, prem, WITH_REDUNDANCE);
- list_int_add (&dejavu, &size, numve, WITH_REDUNDANCE);
- hashtable *table = a2ri_vf_construction_edge_table (m, NULL, 0);
- vf_edge *temp = hashtable_look_for (table, numve, index);
- while (index != -1)
- {
- if (temp->nbsharedfaces != 2)
- return 0;
- fa1 = temp->sharedfaces[0];
- fa2 = temp->sharedfaces[1];
- v1 = a2ri_vf_normal_face (m, fa1);
- v2 = a2ri_vf_normal_face (m, fa2);
- sommeangle += vector3d_angle_radian (&v1, &v2);
- index = -1;
- if (list_int_contains (dejavu, size, m->fa[fa1].ve1) == -1)
- index = m->fa[fa1].ve1;
- if (list_int_contains (dejavu, size, m->fa[fa1].ve2) == -1)
- index = m->fa[fa1].ve2;
- if (list_int_contains (dejavu, size, m->fa[fa1].ve3) == -1)
- index = m->fa[fa1].ve3;
- if (list_int_contains (dejavu, size, m->fa[fa2].ve1) == -1)
- index = m->fa[fa2].ve1;
- if (list_int_contains (dejavu, size, m->fa[fa2].ve2) == -1)
- index = m->fa[fa2].ve2;
- if (list_int_contains (dejavu, size, m->fa[fa2].ve3) == -1)
- index = m->fa[fa2].ve3;
- list_int_add (&dejavu, &size, index, WITH_REDUNDANCE);
- if (index != -1)
- temp = hashtable_look_for (table, numve, index);
- }
- free (dejavu);
- hashtable_free (table);
- free (table);
- return (2.0 * PI - sommeangle);
- }
- /**
- Calcul de la courbure de Levi-Civita
- @param m le modele
- @param ve liste des numéros de sommets
- @param size taille de la liste
- @param list liste des courbures de levi-civita
- @return aucun
- */
- void
- a2ri_vf_levi_civita_list_of_vertex (
- const vf_model * const m,
- const int * const ve,
- int size,
- double ** list)
- {
- int newsize = 0;
- *list = NULL;
- for (int i = 0; i < size; i++)
- list_double_add (list, &newsize, a2ri_vf_levi_civita (m, ve[i]),
- WITH_REDUNDANCE);
- }
- /**
- Calcul des courbures avec la méthode garimella
- @param m le modele
- @param numve numéro du sommet
- @param K gaussian curvature
- @param H mean curvature
- @return aucun
- **/
- void
- a2ri_vf_garimella (
- const vf_model * const m,
- int numve,
- double * K,
- double * H)
- {
- //TODO vérifier le bon fonctionnement de la fonction
- vector3d X,
- Y,
- Z,
- temp;
- int *listface = NULL,
- size_listface = 0;
- int *listfacestar = NULL,
- size_listfacestar = 0;
- int *listvertex = NULL,
- size_listvertex = 0;
- point3d *sommet = NULL,
- aajouter;
- int size_sommet = 0;
- double data2[3],
- *data3,
- *data4,
- data3_temp[3],
- data9_temp[9];
- double a,
- b,
- c,
- d,
- e;
- gsl_matrix *mat_identity,
- *mat_normale,
- *mat_normale_trans,
- *mat_normale_carre,
- *mat_i,
- *mat_sous,
- *mat_mul,
- *mat_R;
- mat_identity = gsl_matrix_calloc (3, 3);
- a2ri_erreur_critique_si (mat_identity == NULL,
- "erreur allocation memoire pour mat_identity\na2ri_vf_garimella");
- gsl_matrix_set_identity (mat_identity);
- //evaluer la normale au sommet numve qui est l'axe Z du nouveau repere
- a2ri_vf_faces_next_vertex (m, numve, &listface, &size_listface);
- vector3d_init (&Z, 0, 0, 0);
- for (int i = 0; i < size_listface; i++)
- {
- temp = a2ri_vf_normal_face (m, listface[i]);
- Z.dx += temp.dx;
- Z.dy += temp.dy;
- Z.dz += temp.dz;
- }
- vector3d_normalize (&Z);
- //vector3d_reverse(&Z);
- data3_temp[0] = Z.dx;
- data3_temp[1] = Z.dy;
- data3_temp[2] = Z.dz;
- mat_normale = matrix_init (data3_temp, 3, 1);
- mat_normale_trans = matrix_transpose (mat_normale);
- mat_normale_carre = matrix_mul (mat_normale, mat_normale_trans);
- if ((Z.dx == 1 || Z.dx == -1) && Z.dy == 0 && Z.dz == 0)
- {
- //cas particulier
- data3_temp[0] = 0;
- data3_temp[1] = 1;
- data3_temp[2] = 0;
- }
- else
- {
- data3_temp[0] = 1;
- data3_temp[1] = 0;
- data3_temp[2] = 0;
- }
- mat_i = matrix_init (data3_temp, 3, 1);
- mat_sous = matrix_sub (mat_identity, mat_normale_carre);
- mat_mul = matrix_mul (mat_sous, mat_i);
- X.dx = gsl_matrix_get (mat_mul, 0, 0);
- X.dy = gsl_matrix_get (mat_mul, 1, 0);
- X.dz = gsl_matrix_get (mat_mul, 2, 0);
- Y = vector3d_vectorialproduct (&X, &Z);
- vector3d_normalize (&X);
- vector3d_normalize (&Y);
- gsl_matrix_free (mat_mul);
- gsl_matrix_free (mat_sous);
- gsl_matrix_free (mat_normale_carre);
- gsl_matrix_free (mat_normale);
- gsl_matrix_free (mat_normale_trans);
- gsl_matrix_free (mat_identity);
- gsl_matrix_free (mat_i);
- data9_temp[0] = X.dx;
- data9_temp[1] = X.dy;
- data9_temp[2] = X.dz;
- data9_temp[3] = Y.dx;
- data9_temp[4] = Y.dy;
- data9_temp[5] = Y.dz;
- data9_temp[6] = Z.dx;
- data9_temp[7] = Z.dy;
- data9_temp[8] = Z.dz;
- mat_R = matrix_init (data9_temp, 3, 3);
- //selectionner au moins 5 voisins
- for (int i = 0; i < size_listface; i++)
- {
- if (m->fa[listface[i]].ve1 != numve)
- list_int_add (&listvertex, &size_listvertex, m->fa[listface[i]].ve1,
- WITHOUT_REDUNDANCE);
- if (m->fa[listface[i]].ve2 != numve)
- list_int_add (&listvertex, &size_listvertex, m->fa[listface[i]].ve2,
- WITHOUT_REDUNDANCE);
- if (m->fa[listface[i]].ve3 != numve)
- list_int_add (&listvertex, &size_listvertex, m->fa[listface[i]].ve3,
- WITHOUT_REDUNDANCE);
- }
- if (size_listvertex < 5)
- {
- free (listvertex);
- listvertex = NULL;
- size_listvertex = 0;
- a2ri_vf_star (BYVERTEX, m, listface, size_listface, &listfacestar,
- &size_listfacestar, 2);
- for (int i = 0; i < size_listfacestar; i++)
- {
- if (m->fa[listfacestar[i]].ve1 != numve)
- list_int_add (&listvertex, &size_listvertex,
- m->fa[listfacestar[i]].ve1, WITHOUT_REDUNDANCE);
- if (m->fa[listfacestar[i]].ve2 != numve)
- list_int_add (&listvertex, &size_listvertex,
- m->fa[listfacestar[i]].ve2, WITHOUT_REDUNDANCE);
- if (m->fa[listfacestar[i]].ve3 != numve)
- list_int_add (&listvertex, &size_listvertex,
- m->fa[listfacestar[i]].ve3, WITHOUT_REDUNDANCE);
- }
- }
- free (listface);
- //changement de repère pour les voisins
- //translation des sommets
- for (int i = 0; i < size_listvertex; i++)
- {
- point3d_init (&aajouter,
- m->ve[listvertex[i]].x - m->ve[numve].x,
- m->ve[listvertex[i]].y - m->ve[numve].y,
- m->ve[listvertex[i]].z - m->ve[numve].z);
- list_point3d_add (&sommet, &size_sommet, &aajouter, WITH_REDUNDANCE);
- }
- free (listvertex);
- //changement de repere avec la matrice R
- for (int i = 0; i < size_sommet; i++)
- {
- data2[0] = sommet[i].x;
- data2[1] = sommet[i].y;
- data2[2] = sommet[i].z;
- gsl_matrix *m = matrix_init (data2, 3, 1);
- gsl_matrix *mul = matrix_mul (mat_R, m);
- sommet[i].x = gsl_matrix_get (mul, 0, 0);
- sommet[i].y = gsl_matrix_get (mul, 1, 0);
- sommet[i].z = gsl_matrix_get (mul, 2, 0);
- gsl_matrix_free (m);
- gsl_matrix_free (mul);
- }
- //resolution avec la methodes des moindres carres
- data3 = (double *) malloc (5 * size_sommet * sizeof (double));
- a2ri_erreur_critique_si (data3 == NULL,
- "erreur allocation memoire pour data3\na2ri_vf_garimella");
- data4 = (double *) malloc (size_sommet * sizeof (double));
- a2ri_erreur_critique_si (data4 == NULL,
- "erreur allocation memoire pour data4\na2ri_vf_garimella");
- for (int i = 0; i < size_sommet; i++)
- {
- data3[i * 5] = sommet[i].x * sommet[i].x;
- data3[i * 5 + 1] = sommet[i].x * sommet[i].y;
- data3[i * 5 + 2] = sommet[i].y * sommet[i].y;
- data3[i * 5 + 3] = sommet[i].x;
- data3[i * 5 + 4] = sommet[i].y;
- data4[i] = sommet[i].z;
- }
- gsl_matrix *A = matrix_init (data3, size_sommet, 5);
- gsl_matrix *B = matrix_init (data4, size_sommet, 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 *x = matrix_mul (pseudoinverse, B);
- a = gsl_matrix_get (x, 0, 0);
- b = gsl_matrix_get (x, 1, 0);
- c = gsl_matrix_get (x, 2, 0);
- d = gsl_matrix_get (x, 3, 0);
- e = gsl_matrix_get (x, 4, 0);
- gsl_matrix_free (A);
- gsl_matrix_free (B);
- gsl_matrix_free (pseudoinverse);
- gsl_matrix_free (x);
- gsl_matrix_free (transA);
- gsl_matrix_free (multransAA);
- gsl_matrix_free (invmultransAA);
- free (sommet);
- //calcul des courbures
- *K = (4 * a * c * -b * b) / ((1 + d * d + e * e) * (1 + d * d + e * e));
- *H =
- (a + c + a * e * e + c * d * d -
- b * d * e) / (pow (1 + d * d + e * e, 1.5));
- }
- /**
- Calcul de la précision, justesse d'un squelette par rapport à un modèle numérique
- @param m le modèle
- @param s le squelette
- @return un nombre définissant la précision du squelette par rapport au modèle. Plus ce nombre sera élevé, moins le squelette sera ajusté au modèle.
- **/
- int a2ri_vf_model_skeleton_accuracy(
- const vf_model * const m,
- const skeleton * const s)
- {
- /*int non_visible=0;
- for(int i=0;i<m->nbvertex;i++)
- for(int j=0;j<s->nbvertex;j++)
- {
- int cpt=0;
- point3d s1,s2,t1,t2,t3;
- double t;
- point3d_init(&s1,m->ve[i].x,m->ve[i].y,m->ve[i].z);
- point3d_init(&s2,s->ve[j].x,s->ve[j].y,s->ve[j].z);
- point3d_init(&t1,m->ve[m->fa[cpt].ve1].x,m->ve[m->fa[cpt].ve1].y,m->ve[m->fa[cpt].ve1].z);
- point3d_init(&t2,m->ve[m->fa[cpt].ve2].x,m->ve[m->fa[cpt].ve2].y,m->ve[m->fa[cpt].ve2].z);
- point3d_init(&t3,m->ve[m->fa[cpt].ve3].x,m->ve[m->fa[cpt].ve3].y,m->ve[m->fa[cpt].ve3].z);
-
- //while(cpt<m->nbface && !intersection_segment_triangle(&s1,&s2,&t1,&t2,&t3,&t))
- // {
- // cpt++;
- // if(cpt!=m->nbface)
- // {
- //point3d_init(&t1,m->ve[m->fa[cpt].ve1].x,m->ve[m->fa[cpt].ve1].y,m->ve[m->fa[cpt].ve1].z);
- // point3d_init(&t2,m->ve[m->fa[cpt].ve2].x,m->ve[m->fa[cpt].ve2].y,m->ve[m->fa[cpt].ve2].z);
- // point3d_init(&t3,m->ve[m->fa[cpt].ve3].x,m->ve[m->fa[cpt].ve3].y,m->ve[m->fa[cpt].ve3].z);
- // }
- // }
- int fin=0;
- while(cpt<m->nbface && !fin)
- {
- printf("-----------------\n");
- printf("point du modele numerique : (%lf , %lf , %lf)\n",s1.x,s1.y,s1.z);
- printf("point du squelette : (%lf , %lf , %lf)\n",s2.x,s2.y,s2.z);
- printf("test face : %d - (%lf , %lf , %lf) (%lf , %lf , %lf) (%lf , %lf , %lf)\n",cpt,t1.x,t1.y,t1.z,t2.x,t2.y,t2.z,t3.x,t3.y,t3.z);
- if(intersection_segment_triangle(&s1,&s2,&t1,&t2,&t3,&t))
- printf("********************************************intersection en t = %lf - (%lf , %lf , %lf)\n",t,s1.x+t*(s2.x-s1.x),s1.y+t*(s2.y-s1.y),s1.z+t*(s2.z-s1.z));
- else
- printf("Pas d'intersection\n");
- printf("-----------------\n");
- if(!intersection_segment_triangle(&s1,&s2,&t1,&t2,&t3,&t))
- fin=1;
-
- cpt++;
- if(cpt!=m->nbface)
- {
- point3d_init(&t1,m->ve[m->fa[cpt].ve1].x,m->ve[m->fa[cpt].ve1].y,m->ve[m->fa[cpt].ve1].z);
- point3d_init(&t2,m->ve[m->fa[cpt].ve2].x,m->ve[m->fa[cpt].ve2].y,m->ve[m->fa[cpt].ve2].z);
- point3d_init(&t3,m->ve[m->fa[cpt].ve3].x,m->ve[m->fa[cpt].ve3].y,m->ve[m->fa[cpt].ve3].z);
- }
- }
-
- if(cpt==m->nbface)
- non_visible++;
- }
- return non_visible;*/
- int non_visible=0;
- point3d s1,s2,t1,t2,t3;
- double t;
-
- for(int i=0;i<m->nbvertex;i++)
- {
- point3d_init(&s1,m->ve[i].x,m->ve[i].y,m->ve[i].z);
- for(int j=0;j<s->nbvertex;j++)
- {
- point3d_init(&s2,s->ve[j].x,s->ve[j].y,s->ve[j].z);
- int nb_intersection=0;
- int k=0;
- while(k<m->nbface && nb_intersection==0)
- {
- point3d_init(&t1,m->ve[m->fa[k].ve1].x,m->ve[m->fa[k].ve1].y,m->ve[m->fa[k].ve1].z);
- point3d_init(&t2,m->ve[m->fa[k].ve2].x,m->ve[m->fa[k].ve2].y,m->ve[m->fa[k].ve2].z);
- point3d_init(&t3,m->ve[m->fa[k].ve3].x,m->ve[m->fa[k].ve3].y,m->ve[m->fa[k].ve3].z);
- if(intersection_segment_triangle(&s1,&s2,&t1,&t2,&t3,&t))
- nb_intersection++;
- k++;
- }
- if(nb_intersection>0)
- non_visible++;
- }
- }
- return non_visible;
- }
|