1 #ifndef _ISO_STATISTICS
2 #define _ISO_STATISTICS
3
4 #include <iso_parametrization.h>
5
6 template <class ScalarType>
geomAverage(const ScalarType & val0,const ScalarType & val1,const ScalarType & exp0,const ScalarType & exp1)7 ScalarType geomAverage(const ScalarType &val0,
8 const ScalarType &val1,
9 const ScalarType &exp0,
10 const ScalarType &exp1)
11 {
12 ScalarType v0=pow(val0,exp0);
13 ScalarType v1=pow(val1,exp1);
14 ScalarType ret=pow((ScalarType)(v0*v1),(ScalarType)(1.0/(exp0+exp1)));
15 return (ret);
16 }
17
18 //template <class MeshType>
19 class Statistic
20 {
21 IsoParametrization *Isoparam;
22 typedef IsoParametrization::PScalarType PScalarType;
23 typedef IsoParametrization::CoordType CoordType;
24 PScalarType Area3d;
25 PScalarType AbstractArea;
26
27 public:
28
AreaDistorsion(ParamFace * f,PScalarType & distorsion,PScalarType & area_3d)29 void AreaDistorsion(ParamFace *f,PScalarType &distorsion,PScalarType &area_3d)
30 {
31 const float epsilon=0.000001f;
32 const float maxRatio=10.f;
33 int indexDomain;
34 vcg::Point2f UV0,UV1,UV2;
35 /*int domainType=*/Isoparam->InterpolationSpace(f,UV0,UV1,UV2,indexDomain);
36
37
38 area_3d=((f->P(1)-f->P(0))^(f->P(2)-f->P(0))).Norm()/(PScalarType)2.0;
39 PScalarType area_2d=((UV1-UV0)^(UV2-UV0))/(PScalarType)2.0;
40 area_3d/=Area3d;
41 area_2d/=AbstractArea;
42
43 if (fabs(area_2d)<epsilon)area_2d=epsilon;
44 if (fabs(area_3d)<epsilon)area_3d=epsilon;
45
46 PScalarType r0=area_3d/area_2d;
47 PScalarType r1=area_2d/area_3d;
48 if (r0>maxRatio)r0=maxRatio;
49 if (r1>maxRatio)r1=maxRatio;
50
51 distorsion=((r0+r1)/(PScalarType)2.0)-(PScalarType)1.0;
52 }
53
AngleDistorsion(ParamFace * f,PScalarType & distortion)54 void AngleDistorsion(ParamFace *f,PScalarType &distortion)
55 {
56 const float epsilon=0.000001f;
57 PScalarType area_3d=((f->P(1)-f->P(0))^(f->P(2)-f->P(0))).Norm();
58 int indexDomain;
59 vcg::Point2f UV0,UV1,UV2;
60 /*int domainType=*/Isoparam->InterpolationSpace(f,UV0,UV1,UV2,indexDomain);
61 PScalarType area_2d=fabs((UV1-UV0)^(UV2-UV0));
62 PScalarType a2=(f->P(1)-f->P(0)).SquaredNorm();
63 PScalarType b2=(f->P(2)-f->P(1)).SquaredNorm();
64 PScalarType c2=(f->P(0)-f->P(2)).SquaredNorm();
65 PScalarType cot_a=((UV2-UV1)*(UV0-UV2));
66 PScalarType cot_b=((UV0-UV2)*(UV1-UV0));
67 PScalarType cot_c=((UV1-UV0)*(UV2-UV1));
68
69 PScalarType num;
70 if ((fabs(area_2d)<epsilon)||(fabs(area_3d)<epsilon))
71 {
72 distortion=0;
73 num=0;
74 }
75 else
76 {
77 num=(cot_a*a2+cot_b*b2+cot_c*c2)/area_2d;
78 distortion=fabs(num/(Area3d*(PScalarType)4.0));
79 assert(distortion<1);
80 }
81 }
82
AreaDistorsion()83 PScalarType AreaDistorsion()
84 {
85 PScalarType sum=0;
86 ParamMesh *param_mesh=Isoparam->ParaMesh();
87 for (unsigned int i=0;i<param_mesh->face.size();i++)
88 {
89 float area3d=0;
90 PScalarType distorsion=0;
91 AreaDistorsion(¶m_mesh->face[i],distorsion,area3d);
92 sum+=distorsion*area3d;
93 }
94 return(fabs(sum));
95 }
96
AngleDistorsion()97 PScalarType AngleDistorsion()
98 {
99 PScalarType sum=0;
100 ParamMesh *param_mesh=Isoparam->ParaMesh();
101 for (unsigned int i=0;i<param_mesh->face.size();i++)
102 {
103 PScalarType distorsion=0;
104 AngleDistorsion(¶m_mesh->face[i],distorsion);
105 sum+=distorsion;
106 }
107 return ((sum)-(PScalarType)1.0);
108 }
109
AggregateDistorsion()110 PScalarType AggregateDistorsion()
111 {
112 PScalarType d0=AreaDistorsion();
113 PScalarType d1=AngleDistorsion();
114 PScalarType ret=geomAverage<PScalarType>(d0+(PScalarType)1.0,d1+(PScalarType)1.0,3,1)-(PScalarType)1;
115 return ret;
116 }
117
PrintAttributes()118 void PrintAttributes()
119 {
120 printf("\n STATISTICS \n"),
121 printf("AREA distorsion:%lf ;\n",AreaDistorsion());
122 printf("ANGLE distorsion:%lf ;\n",AngleDistorsion());
123 printf("AGGREGATE distorsion:%lf ;\n",AggregateDistorsion());
124 printf("L2 STRETCH efficiency:%lf ;\n",L2Error());
125 }
126
127
L2Error()128 PScalarType L2Error()
129 {
130 ///equilateral triagle
131 vcg::Point2f p0(-0.5,0.0);
132 vcg::Point2f p1(0.5,0.0);
133 vcg::Point2f p2(0,0.866025f);
134 ParamMesh::FaceIterator Fi;
135
136 float sum=0;
137 float A3dtot=0;
138 float A2dtot=0;
139 ParamMesh *param_mesh=Isoparam->ParaMesh();
140 for (Fi=param_mesh->face.begin();Fi!=param_mesh->face.end();Fi++)
141 {
142 if (!(*Fi).IsD())
143 {
144 ParamFace *f=&(*Fi);
145 CoordType q1=(*Fi).V(0)->P();
146 CoordType q2=(*Fi).V(1)->P();
147 CoordType q3=(*Fi).V(2)->P();
148 ///map to equilateral triangle
149 /*vcg::Point2f UV1=p0*(*Fi).V(0)->Bary.X()+p1*(*Fi).V(0)->Bary.Y()+p2*(*Fi).V(0)->Bary.Z();
150 vcg::Point2f UV2=p0*(*Fi).V(1)->Bary.X()+p1*(*Fi).V(1)->Bary.Y()+p2*(*Fi).V(1)->Bary.Z();
151 vcg::Point2f UV3=p0*(*Fi).V(2)->Bary.X()+p1*(*Fi).V(2)->Bary.Y()+p2*(*Fi).V(2)->Bary.Z();*/
152 vcg::Point2f UV1,UV2,UV3;
153 int indexDomain;
154 /*int domainType=*/Isoparam->InterpolationSpace(f,UV1,UV2,UV3,indexDomain);
155 PScalarType s1=UV1.X();
156 PScalarType t1=UV1.Y();
157 PScalarType s2=UV2.X();
158 PScalarType t2=UV2.Y();
159 PScalarType s3=UV3.X();
160 PScalarType t3=UV3.Y();
161 PScalarType A=fabs(((s2-s1)*(t3-t1)-(s3-s1)*(t2-t1))/(PScalarType)2.0);
162 if (A<(PScalarType)0.00001)
163 A=(PScalarType)0.00001;
164 PScalarType A3d=((q2 - q1) ^ (q3 - q1)).Norm()/(PScalarType)2.0;
165 A3dtot+=A3d;
166 A2dtot+=A;
167 CoordType Ss=(q1*(t2-t3)+q2*(t3-t1)+q3*(t1-t2))/((PScalarType)2.0*A);
168 CoordType St=(q1*(s3-s2)+q2*(s1-s3)+q3*(s2-s1))/((PScalarType)2.0*A);
169 PScalarType a=Ss*Ss;
170 /*PScalarType b=Ss*St;*/
171 PScalarType c=St*St;
172 PScalarType L2=sqrt((a+c)/(PScalarType)2.0);
173 sum+=L2*L2*A3d;
174 }
175 }
176 sum=sqrt(sum/A3dtot)*sqrt(A2dtot/A3dtot);
177 return sum;
178 }
179
Init(IsoParametrization * _Isoparam)180 void Init(IsoParametrization *_Isoparam)
181 {
182 float fix_num=sqrt((PScalarType)3.0)/(PScalarType)4.0;
183
184 Isoparam=_Isoparam;
185 ParamMesh *param_mesh=Isoparam->ParaMesh();
186 AbstractMesh *abstract_mesh=Isoparam->AbsMesh();
187
188 Area3d=vcg::tri::Stat<ParamMesh>::ComputeMeshArea(*param_mesh);
189 AbstractArea=(PScalarType)abstract_mesh->fn*fix_num;
190 }
191 };
192
193 template <class MeshType>
ApproxAreaDistortion(MeshType & mesh,const int & num_face)194 typename MeshType::ScalarType ApproxAreaDistortion(MeshType &mesh,const int &num_face)
195 {
196 typedef typename MeshType::FaceType FaceType;
197 typedef typename MeshType::ScalarType ScalarType;
198
199 const ScalarType maxRatio=10;
200 const ScalarType epsilon=(ScalarType)0.000001;
201
202 ScalarType sum=0,div=0;
203 ScalarType area_tot=Area<MeshType>(mesh);
204
205 for (unsigned int i=0;i<mesh.face.size();i++)
206 {
207 FaceType* f=&mesh.face[i];
208 if ((f->V(0)->father==f->V(1)->father)&&(f->V(1)->father==f->V(2)->father))
209 {
210 ScalarType area_3d=((f->P(1)-f->P(0))^(f->P(2)-f->P(0))).Norm()/area_tot;
211 ScalarType area_2d=fabs(AreaUV<FaceType>(*f)/((ScalarType)num_face));
212 if (fabs(area_2d)<epsilon)area_2d=epsilon;
213 if (fabs(area_3d)<epsilon)area_3d=epsilon;
214 ScalarType r0=area_3d/area_2d;
215 ScalarType r1=area_2d/area_3d;
216 if (r0>maxRatio)r0=maxRatio;
217 if (r1>maxRatio)r1=maxRatio;
218 sum+=area_3d*(r0+r1);
219 div+=area_3d;
220 }
221 }
222 return (sum/(div*2))-1.0;
223 }
224
225 template <class MeshType>
ApproxAngleDistortion(MeshType & mesh)226 typename MeshType::ScalarType ApproxAngleDistortion(MeshType &mesh)
227 {
228 typedef typename MeshType::FaceType FaceType;
229 typedef typename MeshType::ScalarType ScalarType;
230
231 //const ScalarType maxRatio=4;
232 const ScalarType epsilon=(ScalarType)0.000001;
233
234 ScalarType sum=0,div=0;
235 //ScalarType area_tot=Area<MeshType>(mesh);
236 vcg::Point2<ScalarType> x_axis(0.5f, (ScalarType)(sqrt(3.0)/2.0));
237 vcg::Point2<ScalarType> y_axis(1.0f,0);
238 for (unsigned int i=0;i<mesh.face.size();i++)
239 {
240 FaceType* f=&mesh.face[i];
241 if ((f->V(0)->father==f->V(1)->father)&&(f->V(1)->father==f->V(2)->father))
242 {
243 ScalarType area_3d=((f->P(1)-f->P(0))^(f->P(2)-f->P(0))).Norm();
244 vcg::Point2<ScalarType> t0=f->V(0)->Bary.X()*x_axis+f->V(0)->Bary.Y()*y_axis;
245 vcg::Point2<ScalarType> t1=f->V(1)->Bary.X()*x_axis+f->V(1)->Bary.Y()*y_axis;
246 vcg::Point2<ScalarType> t2=f->V(2)->Bary.X()*x_axis+f->V(2)->Bary.Y()*y_axis;
247 ScalarType area_2d=fabs((t1-t0)^(t2-t0));
248
249
250 ScalarType a2=(f->P(1)-f->P(0)).SquaredNorm();
251 ScalarType b2=(f->P(2)-f->P(1)).SquaredNorm();
252 ScalarType c2=(f->P(0)-f->P(2)).SquaredNorm();
253 ScalarType cot_a=((t2-t1)*(t0-t2));
254 ScalarType cot_b=((t0-t2)*(t1-t0));
255 ScalarType cot_c=((t1-t0)*(t2-t1));
256
257 ScalarType num;
258 if ((fabs(area_2d)<epsilon)||(fabs(area_3d)<epsilon))
259 num=0;
260 else
261 num=(cot_a*a2+cot_b*b2+cot_c*c2)/area_2d;
262
263 sum+=num;
264 //assert(num>=2.0*area_3d);
265 div+=area_3d;
266 }
267 }
268 return (ScalarType)(fabs(sum)/(div*2)-1.0);
269 }
270
271 template <class MeshType>
ApproxL2Error(MeshType & mesh)272 typename MeshType::ScalarType ApproxL2Error(MeshType &mesh)
273 {
274 typedef typename MeshType::ScalarType ScalarType;
275 typedef typename MeshType::CoordType CoordType;
276 ///equilateral triagle
277 vcg::Point2f p0(-0.5,0.0);
278 vcg::Point2f p1(0.5,0.0);
279 vcg::Point2f p2(0,0.866025f);
280 typename MeshType::FaceIterator Fi;
281
282 ScalarType sum=0;
283 ScalarType A3dtot=0;
284 ScalarType A2dtot=0;
285 for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
286 {
287 if (!(*Fi).IsD())
288 {
289 if (((*Fi).V(0)->father==(*Fi).V(1)->father)&&
290 ((*Fi).V(1)->father==(*Fi).V(2)->father))
291 {
292 CoordType q1=(*Fi).V(0)->RPos;
293 CoordType q2=(*Fi).V(1)->RPos;
294 CoordType q3=(*Fi).V(2)->RPos;
295 ///map to equilateral triangle
296 vcg::Point2f UV1=p0*(*Fi).V(0)->Bary.X()+p1*(*Fi).V(0)->Bary.Y()+p2*(*Fi).V(0)->Bary.Z();
297 vcg::Point2f UV2=p0*(*Fi).V(1)->Bary.X()+p1*(*Fi).V(1)->Bary.Y()+p2*(*Fi).V(1)->Bary.Z();
298 vcg::Point2f UV3=p0*(*Fi).V(2)->Bary.X()+p1*(*Fi).V(2)->Bary.Y()+p2*(*Fi).V(2)->Bary.Z();
299 ScalarType s1=UV1.X();
300 ScalarType t1=UV1.Y();
301 ScalarType s2=UV2.X();
302 ScalarType t2=UV2.Y();
303 ScalarType s3=UV3.X();
304 ScalarType t3=UV3.Y();
305 ScalarType A=fabs(((s2-s1)*(t3-t1)-(s3-s1)*(t2-t1))/2.0);
306 if (A<0.00001)
307 A=(ScalarType)0.00001;
308 ScalarType A3d=((q2 - q1) ^ (q3 - q1)).Norm()/2.0;
309 A3dtot+=A3d;
310 A2dtot+=A;
311 CoordType Ss=(q1*(t2-t3)+q2*(t3-t1)+q3*(t1-t2))/(2.0*A);
312 CoordType St=(q1*(s3-s2)+q2*(s1-s3)+q3*(s2-s1))/(2.0*A);
313 ScalarType a=Ss*Ss;
314 /*ScalarType b=Ss*St;*/
315 ScalarType c=St*St;
316 ScalarType L2=sqrt((a+c)/2.0);
317 sum+=L2*L2*A3d;
318 }
319 }
320 }
321 sum=sqrt(sum/A3dtot)*sqrt(A2dtot/A3dtot);
322 return sum;
323 }
324
325 #endif
326