/*************************************/
/* 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 . */
/***************************************************************************/
#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
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;inbvertex;i++)
for(int j=0;jnbvertex;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(cptnbface && !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(cptnbface && !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;inbvertex;i++)
{
point3d_init(&s1,m->ve[i].x,m->ve[i].y,m->ve[i].z);
for(int j=0;jnbvertex;j++)
{
point3d_init(&s2,s->ve[j].x,s->ve[j].y,s->ve[j].z);
int nb_intersection=0;
int k=0;
while(knbface && 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;
}