1 /****************************************************************************
2 * VCGLib o o *
3 * Visual and Computer Graphics Library o o *
4 * _ O _ *
5 * Copyright(C) 2004-2016 \/)\/ *
6 * Visual Computing Lab /\/| *
7 * ISTI - Italian National Research Council | *
8 * \ *
9 * All rights reserved. *
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 * This program is distributed in the hope that it will be useful, *
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
20 * for more details. *
21 * *
22 ****************************************************************************/
23
24 #ifndef __VCGLIB__SAMPLING
25 #define __VCGLIB__SAMPLING
26
27 #include <time.h>
28 #include <vcg/complex/algorithms/closest.h>
29 #include <vcg/space/box3.h>
30 #include <vcg/math/histogram.h>
31 #include <vcg/space/color4.h>
32 #include <vcg/simplex/face/distance.h>
33 #include <vcg/complex/algorithms/update/color.h>
34 #include <vcg/space/index/grid_static_ptr.h>
35 #include <vcg/space/index/aabb_binary_tree/aabb_binary_tree.h>
36 #include <vcg/space/index/octree.h>
37 #include <vcg/space/index/spatial_hashing.h>
38 namespace vcg
39 {
40
41 struct SamplingFlags{
42 enum{
43 HIST = 0x0001,
44 VERTEX_SAMPLING = 0x0002,
45 EDGE_SAMPLING = 0x0004,
46 FACE_SAMPLING = 0x0008,
47 MONTECARLO_SAMPLING = 0x0010,
48 SUBDIVISION_SAMPLING = 0x0020,
49 SIMILAR_SAMPLING = 0x0040,
50 NO_SAMPLING = 0x0070,
51 SAVE_ERROR = 0x0100,
52 INCLUDE_UNREFERENCED_VERTICES = 0x0200,
53 USE_STATIC_GRID = 0x0400,
54 USE_HASH_GRID = 0x0800,
55 USE_AABB_TREE = 0x1000,
56 USE_OCTREE = 0x2000
57 };
58 };
59 // -----------------------------------------------------------------------------------------------
60 template <class MetroMesh>
61 class Sampling
62 {
63 public:
64
65 private:
66 typedef typename MetroMesh::CoordType CoordType;
67 typedef typename MetroMesh::ScalarType ScalarType;
68 typedef typename MetroMesh::VertexType VertexType;
69 typedef typename MetroMesh::VertexPointer VertexPointer;
70 typedef typename MetroMesh::VertexIterator VertexIterator;
71 typedef typename MetroMesh::FaceIterator FaceIterator;
72 typedef typename MetroMesh::FaceType FaceType;
73 typedef typename MetroMesh::FaceContainer FaceContainer;
74
75 typedef GridStaticPtr <FaceType, typename MetroMesh::ScalarType > MetroMeshGrid;
76 typedef SpatialHashTable <FaceType, typename MetroMesh::ScalarType > MetroMeshHash;
77 typedef AABBBinaryTreeIndex <FaceType, typename MetroMesh::ScalarType, vcg::EmptyClass> MetroMeshAABB;
78 typedef Octree <FaceType, typename MetroMesh::ScalarType > MetroMeshOctree;
79
80 typedef Point3<typename MetroMesh::ScalarType> Point3x;
81
82
83
84
85 // data structures
86 MetroMesh &S1;
87 MetroMesh &S2;
88 MetroMeshGrid gS2;
89 MetroMeshHash hS2;
90 MetroMeshAABB tS2;
91 MetroMeshOctree oS2;
92
93
94 unsigned int n_samples_per_face ;
95 float n_samples_edge_to_face_ratio ;
96 float bbox_factor ;
97 float inflate_percentage ;
98 unsigned int min_size ;
99 int n_hist_bins ;
100 int print_every_n_elements ;
101 int referredBit ;
102 // parameters
103 double dist_upper_bound;
104 double n_samples_per_area_unit;
105 unsigned long n_samples_target;
106 int Flags;
107
108 // results
109 Histogram<double> hist;
110 unsigned long n_total_samples;
111 unsigned long n_total_area_samples;
112 unsigned long n_total_edge_samples;
113 unsigned long n_total_vertex_samples;
114 double max_dist;
115 double mean_dist;
116 double RMS_dist;
117 double volume;
118 double area_S1;
119
120 // globals
121 int n_samples;
122
123 // private methods
124 inline double ComputeMeshArea(MetroMesh & mesh);
125 float AddSample(const Point3x &p);
126 inline void AddRandomSample(FaceIterator &T);
127 inline void SampleEdge(const Point3x & v0, const Point3x & v1, int n_samples_per_edge);
128 void VertexSampling();
129 void EdgeSampling();
130 void FaceSubdiv(const Point3x & v0, const Point3x &v1, const Point3x & v2, int maxdepth);
131 void SimilarTriangles(const Point3x &v0, const Point3x &v1, const Point3x &v2, int n_samples_per_edge);
132 void MontecarloFaceSampling();
133 void SubdivFaceSampling();
134 void SimilarFaceSampling();
135
136 public :
137 // public methods
138 Sampling(MetroMesh &_s1, MetroMesh &_s2);
139 ~Sampling();
140 void Hausdorff();
GetArea()141 double GetArea() {return area_S1;}
GetDistMax()142 double GetDistMax() {return max_dist;}
GetDistMean()143 double GetDistMean() {return mean_dist;}
GetDistRMS()144 double GetDistRMS() {return RMS_dist;}
GetDistVolume()145 double GetDistVolume() {return volume;}
GetNSamples()146 unsigned long GetNSamples() {return n_total_samples;}
GetNAreaSamples()147 unsigned long GetNAreaSamples() {return n_total_area_samples;}
GetNEdgeSamples()148 unsigned long GetNEdgeSamples() {return n_total_edge_samples;}
GetNVertexSamples()149 unsigned long GetNVertexSamples() {return n_total_vertex_samples;}
GetNSamplesPerAreaUnit()150 double GetNSamplesPerAreaUnit() {return n_samples_per_area_unit;}
GetNSamplesTarget()151 unsigned long GetNSamplesTarget() {return n_samples_target;}
GetHist()152 Histogram<double> &GetHist() {return hist;}
SetFlags(int flags)153 void SetFlags(int flags) {Flags = flags;}
ClearFlag(int flag)154 void ClearFlag(int flag) {Flags &= (flag ^ -1);}
SetParam(double _n_samp)155 void SetParam(double _n_samp) {n_samples_target = _n_samp;}
156 void SetSamplesTarget(unsigned long _n_samp);
157 void SetSamplesPerAreaUnit(double _n_samp);
158 };
159
160 // -----------------------------------------------------------------------------------------------
161
162 // constructor
163 template <class MetroMesh>
Sampling(MetroMesh & _s1,MetroMesh & _s2)164 Sampling<MetroMesh>::Sampling(MetroMesh &_s1, MetroMesh &_s2):S1(_s1),S2(_s2)
165 {
166 Flags = 0;
167 area_S1 = ComputeMeshArea(_s1);
168 // set default numbers
169 n_samples_per_face = 10;
170 n_samples_edge_to_face_ratio = 0.1f;
171 bbox_factor = 0.1f;
172 inflate_percentage = 0.02f;
173 min_size = 125; /* 125 = 5^3 */
174 n_hist_bins = 256;
175 print_every_n_elements = S1.fn/100;
176
177 if(print_every_n_elements <= 1)
178 print_every_n_elements = 2;
179
180 referredBit = VertexType::NewBitFlag();
181 // store the unreferred vertices
182 FaceIterator fi; VertexIterator vi; int i;
183 for(fi = _s1.face.begin(); fi!= _s1.face.end(); ++fi)
184 for(i=0;i<3;++i) (*fi).V(i)->SetUserBit(referredBit);
185 }
186
187 template <class MetroMesh>
~Sampling()188 Sampling<MetroMesh>::~Sampling()
189 {
190 VertexType::DeleteBitFlag(referredBit);
191 }
192
193
194 // set sampling parameters
195 template <class MetroMesh>
SetSamplesTarget(unsigned long _n_samp)196 void Sampling<MetroMesh>::SetSamplesTarget(unsigned long _n_samp)
197 {
198 n_samples_target = _n_samp;
199 n_samples_per_area_unit = n_samples_target / (double)area_S1;
200 }
201
202 template <class MetroMesh>
SetSamplesPerAreaUnit(double _n_samp)203 void Sampling<MetroMesh>::SetSamplesPerAreaUnit(double _n_samp)
204 {
205 n_samples_per_area_unit = _n_samp;
206 n_samples_target = (unsigned long)((double) n_samples_per_area_unit * area_S1);
207 }
208
209
210 // auxiliary functions
211 template <class MetroMesh>
ComputeMeshArea(MetroMesh & mesh)212 inline double Sampling<MetroMesh>::ComputeMeshArea(MetroMesh & mesh)
213 {
214 FaceIterator face;
215 double area = 0.0;
216
217 for(face=mesh.face.begin(); face != mesh.face.end(); face++)
218 if(!(*face).IsD())
219 area += DoubleArea(*face);
220
221 return area/2.0;
222 }
223
224 template <class MetroMesh>
AddSample(const Point3x & p)225 float Sampling<MetroMesh>::AddSample(const Point3x &p )
226 {
227 FaceType *f=0;
228 Point3x normf, bestq, ip;
229 ScalarType dist;
230
231 dist = dist_upper_bound;
232
233 // compute distance between p_i and the mesh S2
234 if(Flags & SamplingFlags::USE_AABB_TREE)
235 f=tri::GetClosestFaceEP<MetroMesh,MetroMeshAABB>(S2, tS2, p, dist_upper_bound, dist, normf, bestq, ip);
236 if(Flags & SamplingFlags::USE_HASH_GRID)
237 f=tri::GetClosestFaceEP<MetroMesh,MetroMeshHash>(S2, hS2, p, dist_upper_bound, dist, normf, bestq, ip);
238 if(Flags & SamplingFlags::USE_STATIC_GRID)
239 f=tri::GetClosestFaceEP<MetroMesh,MetroMeshGrid>(S2, gS2, p, dist_upper_bound, dist, normf, bestq, ip);
240 if (Flags & SamplingFlags::USE_OCTREE)
241 f=tri::GetClosestFaceEP<MetroMesh,MetroMeshOctree>(S2, oS2, p, dist_upper_bound, dist, normf, bestq, ip);
242
243 // update distance measures
244 if(dist == dist_upper_bound)
245 return -1.0;
246
247 if(dist > max_dist)
248 max_dist = dist; // L_inf
249 mean_dist += dist; // L_1
250 RMS_dist += dist*dist; // L_2
251 n_total_samples++;
252
253 if(Flags & SamplingFlags::HIST)
254 hist.Add((float)fabs(dist));
255
256 return (float)dist;
257 }
258
259
260 // -----------------------------------------------------------------------------------------------
261 // --- Vertex Sampling ---------------------------------------------------------------------------
262
263 template <class MetroMesh>
VertexSampling()264 void Sampling<MetroMesh>::VertexSampling()
265 {
266 // Vertex sampling.
267 int cnt = 0;
268 float error;
269
270 printf("Vertex sampling\n");
271 VertexIterator vi;
272 typename std::vector<VertexPointer>::iterator vif;
273 for(vi=S1.vert.begin();vi!=S1.vert.end();++vi)
274 if( (*vi).IsUserBit(referredBit) || // it is referred
275 ((Flags&SamplingFlags::INCLUDE_UNREFERENCED_VERTICES) != 0) ) //include also unreferred
276 {
277 error = AddSample((*vi).cP());
278
279 n_total_vertex_samples++;
280
281 // save vertex quality
282 if(Flags & SamplingFlags::SAVE_ERROR) (*vi).Q() = error;
283
284 // print progress information
285 if(!(++cnt % print_every_n_elements))
286 printf("Sampling vertices %d%%\r", (100 * cnt/S1.vn));
287 }
288 printf(" \r");
289 }
290
291
292 // -----------------------------------------------------------------------------------------------
293 // --- Edge Sampling -----------------------------------------------------------------------------
294
295 template <class MetroMesh>
SampleEdge(const Point3x & v0,const Point3x & v1,int n_samples_per_edge)296 inline void Sampling<MetroMesh>::SampleEdge(const Point3x & v0, const Point3x & v1, int n_samples_per_edge)
297 {
298 // uniform sampling of the segment v0v1.
299 Point3x e((v1-v0)/(double)(n_samples_per_edge+1));
300 int i;
301
302 for(i=1; i <= n_samples_per_edge; i++)
303 {
304 AddSample(v0 + e*i);
305 n_total_edge_samples++;
306 }
307 }
308
309
310 template <class MetroMesh>
EdgeSampling()311 void Sampling<MetroMesh>::EdgeSampling()
312 {
313 // Edge sampling.
314 typedef std::pair<VertexPointer, VertexPointer> pvv;
315 std::vector< pvv > Edges;
316
317 printf("Edge sampling\n");
318
319 // compute edge list.
320 FaceIterator fi;
321 for(fi=S1.face.begin(); fi != S1.face.end(); fi++)
322 for(int i=0; i<3; ++i)
323 {
324 Edges.push_back(std::make_pair((*fi).V0(i),(*fi).V1(i)));
325 if(Edges.back().first > Edges.back().second)
326 std::swap(Edges.back().first, Edges.back().second);
327 }
328 sort(Edges.begin(), Edges.end());
329 typename std::vector< pvv>::iterator edgeend = unique(Edges.begin(), Edges.end());
330 Edges.resize(edgeend-Edges.begin());
331
332 // sample edges.
333 typename std::vector<pvv>::iterator ei;
334 double n_samples_per_length_unit;
335 double n_samples_decimal = 0.0;
336 int cnt=0;
337 if(Flags & SamplingFlags::FACE_SAMPLING)
338 n_samples_per_length_unit = sqrt((double)n_samples_per_area_unit);
339 else
340 n_samples_per_length_unit = n_samples_per_area_unit;
341 for(ei=Edges.begin(); ei!=Edges.end(); ++ei)
342 {
343 n_samples_decimal += Distance((*ei).first->cP(),(*ei).second->cP()) * n_samples_per_length_unit;
344 n_samples = (int) n_samples_decimal;
345 SampleEdge((*ei).first->cP(), (*ei).second->cP(), (int) n_samples);
346 n_samples_decimal -= (double) n_samples;
347
348 // print progress information
349 if(!(++cnt % print_every_n_elements))
350 printf("Sampling edge %lu%%\r", (100 * cnt/Edges.size()));
351 }
352 printf(" \r");
353 }
354
355
356 // -----------------------------------------------------------------------------------------------
357 // --- Face Sampling -----------------------------------------------------------------------------
358
359 // Montecarlo sampling.
360 template <class MetroMesh>
AddRandomSample(FaceIterator & T)361 inline void Sampling<MetroMesh>::AddRandomSample(FaceIterator &T)
362 {
363 // random sampling over the input face.
364 double rnd_1, rnd_2;
365
366 // vertices of the face T.
367 Point3x p0(T->V(0)->cP());
368 Point3x p1(T->V(1)->cP());
369 Point3x p2(T->V(2)->cP());
370 // calculate two edges of T.
371 Point3x v1(p1 - p0);
372 Point3x v2(p2 - p0);
373
374 // choose two random numbers.
375 rnd_1 = (double)rand() / (double)RAND_MAX;
376 rnd_2 = (double)rand() / (double)RAND_MAX;
377 if(rnd_1 + rnd_2 > 1.0)
378 {
379 rnd_1 = 1.0 - rnd_1;
380 rnd_2 = 1.0 - rnd_2;
381 }
382
383 // add a random point on the face T.
384 AddSample (p0 + (v1 * rnd_1 + v2 * rnd_2));
385 n_total_area_samples++;
386 }
387
388 template <class MetroMesh>
MontecarloFaceSampling()389 void Sampling<MetroMesh>::MontecarloFaceSampling()
390 {
391 // Montecarlo sampling.
392 double n_samples_decimal = 0.0;
393 FaceIterator fi;
394
395 srand(clock());
396 // printf("Montecarlo face sampling\n");
397 for(fi=S1.face.begin(); fi != S1.face.end(); fi++)
398 if(!(*fi).IsD())
399 {
400 // compute # samples in the current face.
401 n_samples_decimal += 0.5*DoubleArea(*fi) * n_samples_per_area_unit;
402 n_samples = (int) n_samples_decimal;
403
404 // for every sample p_i in T...
405 for(int i=0; i < n_samples; i++)
406 AddRandomSample(fi);
407
408 n_samples_decimal -= (double) n_samples;
409
410 // print progress information
411 // if(!(++cnt % print_every_n_elements))
412 // printf("Sampling face %d%%\r", (100 * cnt/S1.fn));
413 }
414 // printf(" \r");
415 }
416
417
418 // Subdivision sampling.
419 template <class MetroMesh>
FaceSubdiv(const Point3x & v0,const Point3x & v1,const Point3x & v2,int maxdepth)420 void Sampling<MetroMesh>::FaceSubdiv(const Point3x & v0, const Point3x & v1, const Point3x & v2, int maxdepth)
421 {
422 // recursive face subdivision.
423 if(maxdepth == 0)
424 {
425 // ground case.
426 AddSample((v0+v1+v2)/3.0f);
427 n_total_area_samples++;
428 n_samples++;
429 return;
430 }
431
432 // compute the longest edge.
433 double maxd01 = SquaredDistance(v0,v1);
434 double maxd12 = SquaredDistance(v1,v2);
435 double maxd20 = SquaredDistance(v2,v0);
436 int res;
437 if(maxd01 > maxd12)
438 if(maxd01 > maxd20) res = 0;
439 else res = 2;
440 else
441 if(maxd12 > maxd20) res = 1;
442 else res = 2;
443
444 // break the input triangle along the median to the the longest edge.
445 Point3x pp;
446 switch(res)
447 {
448 case 0 : pp = (v0+v1)/2;
449 FaceSubdiv(v0,pp,v2,maxdepth-1);
450 FaceSubdiv(pp,v1,v2,maxdepth-1);
451 break;
452 case 1 : pp = (v1+v2)/2;
453 FaceSubdiv(v0,v1,pp,maxdepth-1);
454 FaceSubdiv(v0,pp,v2,maxdepth-1);
455 break;
456 case 2 : pp = (v2+v0)/2;
457 FaceSubdiv(v0,v1,pp,maxdepth-1);
458 FaceSubdiv(pp,v1,v2,maxdepth-1);
459 break;
460 }
461 }
462
463 template <class MetroMesh>
SubdivFaceSampling()464 void Sampling<MetroMesh>::SubdivFaceSampling()
465 {
466 // Subdivision sampling.
467 int cnt = 0, maxdepth;
468 double n_samples_decimal = 0.0;
469 typename MetroMesh::FaceIterator fi;
470
471 printf("Subdivision face sampling\n");
472 for(fi=S1.face.begin(); fi != S1.face.end(); fi++)
473 {
474 // compute # samples in the current face.
475 n_samples_decimal += 0.5*DoubleArea(*fi) * n_samples_per_area_unit;
476 n_samples = (int) n_samples_decimal;
477 if(n_samples)
478 {
479 // face sampling.
480 maxdepth = ((int)(log((double)n_samples)/log(2.0)));
481 n_samples = 0;
482 FaceSubdiv((*fi).V(0)->cP(), (*fi).V(1)->cP(), (*fi).V(2)->cP(), maxdepth);
483 }
484 n_samples_decimal -= (double) n_samples;
485
486 // print progress information
487 if(!(++cnt % print_every_n_elements))
488 printf("Sampling face %d%%\r", (100 * cnt/S1.fn));
489 }
490 printf(" \r");
491 }
492
493
494 // Similar Triangles sampling.
495 template <class MetroMesh>
SimilarTriangles(const Point3x & v0,const Point3x & v1,const Point3x & v2,int n_samples_per_edge)496 void Sampling<MetroMesh>::SimilarTriangles(const Point3x & v0, const Point3x & v1, const Point3x & v2, int n_samples_per_edge)
497 {
498 Point3x V1((v1-v0)/(double)(n_samples_per_edge-1));
499 Point3x V2((v2-v0)/(double)(n_samples_per_edge-1));
500 int i, j;
501
502 // face sampling.
503 for(i=1; i < n_samples_per_edge-1; i++)
504 for(j=1; j < n_samples_per_edge-1-i; j++)
505 {
506 AddSample( v0 + (V1*(double)i + V2*(double)j) );
507 n_total_area_samples++;
508 n_samples++;
509 }
510 }
511
512 template <class MetroMesh>
SimilarFaceSampling()513 void Sampling<MetroMesh>::SimilarFaceSampling()
514 {
515 // Similar Triangles sampling.
516 int cnt = 0, n_samples_per_edge;
517 double n_samples_decimal = 0.0;
518 FaceIterator fi;
519
520 printf("Similar Triangles face sampling\n");
521 for(fi=S1.face.begin(); fi != S1.face.end(); fi++)
522 {
523 // compute # samples in the current face.
524 n_samples_decimal += 0.5*DoubleArea(*fi) * n_samples_per_area_unit;
525 n_samples = (int) n_samples_decimal;
526 if(n_samples)
527 {
528 // face sampling.
529 n_samples_per_edge = (int)((sqrt(1.0+8.0*(double)n_samples) +5.0)/2.0);
530 n_samples = 0;
531 SimilarTriangles((*fi).V(0)->cP(), (*fi).V(1)->cP(), (*fi).V(2)->cP(), n_samples_per_edge);
532 }
533 n_samples_decimal -= (double) n_samples;
534
535 // print progress information
536 if(!(++cnt % print_every_n_elements))
537 printf("Sampling face %d%%\r", (100 * cnt/S1.fn));
538 }
539 printf(" \r");
540 }
541
542
543 // -----------------------------------------------------------------------------------------------
544 // --- Distance ----------------------------------------------------------------------------------
545
546 template <class MetroMesh>
Hausdorff()547 void Sampling<MetroMesh>::Hausdorff()
548 {
549 Box3< ScalarType> bbox;
550
551 // set grid meshes.
552 if(Flags & SamplingFlags::USE_HASH_GRID) hS2.Set(S2.face.begin(),S2.face.end());
553 if(Flags & SamplingFlags::USE_AABB_TREE) tS2.Set(S2.face.begin(),S2.face.end());
554 if(Flags & SamplingFlags::USE_STATIC_GRID) gS2.Set(S2.face.begin(),S2.face.end());
555 if(Flags & SamplingFlags::USE_OCTREE) oS2.Set(S2.face.begin(),S2.face.end());
556
557 // set bounding box
558 bbox = S2.bbox;
559 dist_upper_bound = /*bbox_factor * */bbox.Diag();
560 if(Flags & SamplingFlags::HIST)
561 hist.SetRange(0.0, dist_upper_bound/100.0, n_hist_bins);
562
563 // initialize sampling statistics.
564 n_total_area_samples = n_total_edge_samples = n_total_vertex_samples = n_total_samples = n_samples = 0;
565 max_dist = -HUGE_VAL;
566 mean_dist = RMS_dist = 0;
567
568 // Vertex sampling.
569 if(Flags & SamplingFlags::VERTEX_SAMPLING)
570 VertexSampling();
571 // Edge sampling.
572 if(n_samples_target > n_total_samples)
573 {
574 n_samples_target -= (int) n_total_samples;
575 n_samples_per_area_unit = n_samples_target / area_S1;
576 if(Flags & SamplingFlags::EDGE_SAMPLING)
577 {
578 EdgeSampling();
579 if(n_samples_target > n_total_samples) n_samples_target -= (int) n_total_samples;
580 else n_samples_target=0;
581 }
582 // Face sampling.
583 if((Flags & SamplingFlags::FACE_SAMPLING) && (n_samples_target > 0))
584 {
585 n_samples_per_area_unit = n_samples_target / area_S1;
586 if(Flags & SamplingFlags::MONTECARLO_SAMPLING) MontecarloFaceSampling();
587 if(Flags & SamplingFlags::SUBDIVISION_SAMPLING) SubdivFaceSampling();
588 if(Flags & SamplingFlags::SIMILAR_SAMPLING) SimilarFaceSampling();
589 }
590 }
591
592 // compute vertex colour
593 if(Flags & SamplingFlags::SAVE_ERROR)
594 vcg::tri::UpdateColor<MetroMesh>::PerVertexQualityRamp(S1);
595
596 // compute statistics
597 n_samples_per_area_unit = (double) n_total_samples / area_S1;
598 volume = mean_dist / n_samples_per_area_unit / 2.0;
599 mean_dist /= n_total_samples;
600 RMS_dist = sqrt(RMS_dist / n_total_samples);
601 }
602 }
603 #endif
604