-
Notifications
You must be signed in to change notification settings - Fork 10
/
metro_hausdorff.cpp
197 lines (174 loc) · 7.74 KB
/
metro_hausdorff.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
// This file is part of the implementation of
// Robust Structure Simplification for Hex Re-meshing
// Xifeng Gao, Daniele Panozzo, Wenping Wang, Zhigang Deng, Guoning Chen
// In ACM Transactions on Graphics (Proceedings of SIGGRAPH ASIA 2017)
//
// Copyright (C) 2017 Xifeng Gao<gxf.xisha@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla Public License
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
// obtain one at http://mozilla.org/MPL/2.0/.
#include "metro_hausdorff.h"
////////////////// Command line Flags and parameters
// -----------------------------------------------------------------------------------------------
// simple aux function that compute the name for the file containing the stored computations
std::string SaveFileName(const std::string &filename)
{
int pos = filename.find_last_of('.', filename.length());
std::string fileout = filename.substr(0, pos) + "_metro.ply";
return fileout;
}
int compute(Mesh &mesh0, Mesh & mesh1, double &hausdorff_ratio, double &hausdorff_ratio_threshold){
#define MSG_ERR_MESH_LOAD "error loading the input meshes.\n"
#define MSG_ERR_INVALID_OPTION "unable to parse option '%s'\n"
#define MSG_ERR_FILE_OPEN "unable to open the output file.'n"
#define MSG_ERR_UNKNOWN_FORMAT "unknown file format '%s'.\n"
// global constants
#define NO_SAMPLES_PER_FACE 10
#define N_SAMPLES_EDGE_TO_FACE_RATIO 0.1
#define BBOX_FACTOR 0.1
#define INFLATE_PERCENTAGE 0.02
#define MIN_SIZE 125 /* 125 = 5^3 */
#define N_HIST_BINS 256
#define PRINT_EVERY_N_ELEMENTS 1000
bool NumberOfSamples = false;
bool SamplesPerAreaUnit = false;
bool CleaningFlag = false;
CMesh S1, S2;
float ColorMin = 0, ColorMax = 0;
double dist1_max, dist2_max;
unsigned long n_samples_target, elapsed_time;
double n_samples_per_area_unit;
int flags;
// default parameters
flags = SamplingFlags::VERTEX_SAMPLING |
SamplingFlags::EDGE_SAMPLING |
SamplingFlags::FACE_SAMPLING |
SamplingFlags::SIMILAR_SAMPLING;
if (!(flags & SamplingFlags::USE_HASH_GRID) && !(flags & SamplingFlags::USE_AABB_TREE) && !(flags & SamplingFlags::USE_OCTREE))
flags |= SamplingFlags::USE_STATIC_GRID;
// load input meshes.
//OpenMesh("",S1);
//OpenMesh(argv[2],S2);
//string S1NewName = SaveFileName(argv[1]);
//string S2NewName = SaveFileName(argv[2]);
//tri::io::Importer<CMesh>::Open(m, filename);
S1.vert.resize(mesh0.V.cols());
for (int i = 0; i < mesh0.V.cols(); i++) {
CVertex v;
v.P()[0] = mesh0.V(0, i);
v.P()[1] = mesh0.V(1, i);
v.P()[2] = mesh0.V(2, i);
S1.vert[i] = v;
}
S1.face.resize(mesh0.Fs.size());
//Allocator<CMesh>::AddFaces(S1, mesh0.Fs.size());
for (int i = 0; i < mesh0.Fs.size(); i++) {
CFace f;
f.V(0) = &(S1.vert[mesh0.Fs[i].vs[0]]);
f.V(1) = &(S1.vert[mesh0.Fs[i].vs[1]]);
f.V(2) = &(S1.vert[mesh0.Fs[i].vs[2]]);
S1.face[i] = f;
}
S2.vert.resize(mesh1.V.cols());
for (int i = 0; i < mesh1.V.cols(); i++) {
CVertex v;
v.P()[0] = mesh1.V(0, i);
v.P()[1] = mesh1.V(1, i);
v.P()[2] = mesh1.V(2, i);
S2.vert[i] = v;
}
S2.face.resize(mesh1.Fs.size());
//Allocator<CMesh>::AddFaces(S1, mesh0.Fs.size());
for (int i = 0; i < mesh1.Fs.size(); i++) {
CFace f;
f.V(0) = &(S2.vert[mesh1.Fs[i].vs[0]]);
f.V(1) = &(S2.vert[mesh1.Fs[i].vs[1]]);
f.V(2) = &(S2.vert[mesh1.Fs[i].vs[2]]);
S2.face[i] = f;
}
S1.vn = S1.vert.size();
S1.fn = S1.face.size();
S2.vn = S2.vert.size();
S2.fn = S2.face.size();
if (!NumberOfSamples && !SamplesPerAreaUnit) {
NumberOfSamples = true;
n_samples_target = 10 * max(S1.fn, S2.fn);// take 10 samples per face
}
// compute face information
tri::UpdateComponentEP<CMesh>::Set(S1);
tri::UpdateComponentEP<CMesh>::Set(S2);
// set bounding boxes for S1 and S2
tri::UpdateBounding<CMesh>::Box(S1);
tri::UpdateBounding<CMesh>::Box(S2);
// set Bounding Box.
Box3<CMesh::ScalarType> bbox, tmp_bbox_M1 = S1.bbox, tmp_bbox_M2 = S2.bbox;
bbox.Add(S1.bbox);
bbox.Add(S2.bbox);
bbox.Offset(bbox.Diag()*0.02);
S1.bbox = bbox;
S2.bbox = bbox;
// initialize time info.
int t0 = clock();
Sampling<CMesh> ForwardSampling(S1, S2);
Sampling<CMesh> BackwardSampling(S2, S1);
// print mesh info.
//printf("Mesh info:\n");
//printf(" M1: '%s'\n\tvertices %7i\n\tfaces %7i\n\tarea %12.4f\n", argv[1], S1.vn, S1.fn, ForwardSampling.GetArea());
//printf("\tbbox (%7.4f %7.4f %7.4f)-(%7.4f %7.4f %7.4f)\n", tmp_bbox_M1.min[0], tmp_bbox_M1.min[1], tmp_bbox_M1.min[2], tmp_bbox_M1.max[0], tmp_bbox_M1.max[1], tmp_bbox_M1.max[2]);
//printf("\tbbox diagonal %f\n", (float)tmp_bbox_M1.Diag());
//printf(" M2: '%s'\n\tvertices %7i\n\tfaces %7i\n\tarea %12.4f\n", argv[2], S2.vn, S2.fn, BackwardSampling.GetArea());
//printf("\tbbox (%7.4f %7.4f %7.4f)-(%7.4f %7.4f %7.4f)\n", tmp_bbox_M2.min[0], tmp_bbox_M2.min[1], tmp_bbox_M2.min[2], tmp_bbox_M2.max[0], tmp_bbox_M2.max[1], tmp_bbox_M2.max[2]);
//printf("\tbbox diagonal %f\n", (float)tmp_bbox_M2.Diag());
// Forward distance.
//printf("\nForward distance (M1 -> M2):\n");
ForwardSampling.SetFlags(flags);
if (NumberOfSamples){
ForwardSampling.SetSamplesTarget(n_samples_target);
n_samples_per_area_unit = ForwardSampling.GetNSamplesPerAreaUnit();
}
else{
ForwardSampling.SetSamplesPerAreaUnit(n_samples_per_area_unit);
n_samples_target = ForwardSampling.GetNSamplesTarget();
}
//printf("target # samples : %lu\ntarget # samples/area : %f\n", n_samples_target, n_samples_per_area_unit);
ForwardSampling.Hausdorff();
dist1_max = ForwardSampling.GetDistMax();
//printf("\ndistances:\n max : %f (%f wrt bounding box diagonal)\n", (float)dist1_max, (float)dist1_max / bbox.Diag());
//printf(" mean : %f\n", ForwardSampling.GetDistMean());
//printf(" RMS : %f\n", ForwardSampling.GetDistRMS());
//printf("# vertex samples %9lu\n", ForwardSampling.GetNVertexSamples());
//printf("# edge samples %9lu\n", ForwardSampling.GetNEdgeSamples());
//printf("# area samples %9lu\n", ForwardSampling.GetNAreaSamples());
//printf("# total samples %9lu\n", ForwardSampling.GetNSamples());
//printf("# samples per area unit: %f\n\n", ForwardSampling.GetNSamplesPerAreaUnit());
// Backward distance.
//printf("\nBackward distance (M2 -> M1):\n");
BackwardSampling.SetFlags(flags);
if (NumberOfSamples){
BackwardSampling.SetSamplesTarget(n_samples_target);
n_samples_per_area_unit = BackwardSampling.GetNSamplesPerAreaUnit();
}
else{
BackwardSampling.SetSamplesPerAreaUnit(n_samples_per_area_unit);
n_samples_target = BackwardSampling.GetNSamplesTarget();
}
//printf("target # samples : %lu\ntarget # samples/area : %f\n", n_samples_target, n_samples_per_area_unit);
BackwardSampling.Hausdorff();
dist2_max = BackwardSampling.GetDistMax();
//printf("\ndistances:\n max : %f (%f wrt bounding box diagonal)\n", (float)dist2_max, (float)dist2_max / bbox.Diag());
//printf(" mean : %f\n", BackwardSampling.GetDistMean());
//printf(" RMS : %f\n", BackwardSampling.GetDistRMS());
//printf("# vertex samples %9lu\n", BackwardSampling.GetNVertexSamples());
//printf("# edge samples %9lu\n", BackwardSampling.GetNEdgeSamples());
//printf("# area samples %9lu\n", BackwardSampling.GetNAreaSamples());
//printf("# total samples %9lu\n", BackwardSampling.GetNSamples());
//printf("# samples per area unit: %f\n\n", BackwardSampling.GetNSamplesPerAreaUnit());
// compute time info.
elapsed_time = clock() - t0;
int n_total_sample = ForwardSampling.GetNSamples() + BackwardSampling.GetNSamples();
double mesh_dist_max = max(dist1_max, dist2_max);
hausdorff_ratio = (float)mesh_dist_max / bbox.Diag();
if (hausdorff_ratio > hausdorff_ratio_threshold) return false;
return true;
}