OMToolkit 1.0 The polygonal mesh processing tool.

- Main Page
- Related Pages
- Namespaces
- Classes
- Files

- File List
- File Members

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
OMProjector.hxx
Go to the documentation of this file.
1 //==============================================================================
13 #ifndef _OM_PROJECTOR_HXX_
14 #define _OM_PROJECTOR_HXX_
15
17 // Constructor - initializes vital variables
18 // @param mesh Pointer to a mesh to work with
20 template <class Mesh, class Matrix>
21 OMProjector<Mesh, Matrix>::OMProjector(MeshT *mesh, OpenMesh::VPropHandleT<MatrixT> propertyHandle) : tree(mesh)
22 {
23 m_mesh = mesh;
24 m_zDir = ZDIR_Z;
25 m_propertyHandle = propertyHandle;
26 }
27
29 // Method computes tangent rasters on all the vertices
30 // @param length Length of a square raster (real length in mesh space)
31 // @param resolution Number of pixels on MatrixT edge
32 // @param zDir Use ZDIR consts for choose which component will be used as Z direction
33 // @param xdirection Direction of X axis on MatrixT (for ex. Maximum curvature direction)
34 // @param lengthRelative If true, length is computed as % of median edges lengths, if not, length is static
35 // @return True, if computation ended successfully, false otherwise
37 template <class Mesh, class Matrix>
38 bool OMProjector<Mesh, Matrix>::ComputeMatrices(ScalarT length, ScalarT resolution, int zDir, bool lengthRelative, int xdirection)
39 {
40 if (!m_mesh->IsTriMesh) return false;
41
42 // init auxiliary variables and check normals
43 m_resolution = resolution;
44
45 VectorT direction;
46
47 // precomp if projecting curvature as ZDir
48 if (zDir == ZDIR_CURVATURE)
49 {
50 m_zDir = ZDIR_CURVATURE;
51 m_mesh->request_face_normals();
52 m_mesh->request_vertex_normals();
53 MeshT::FaceIter end = m_mesh->faces_end();
54 VectorT normal;
55 VectorT firstVec;
56 VectorT secondVec;
57
58 for (MeshT::FaceIter face = m_mesh->faces_begin(); face != end; ++face)
59 {
60 MeshT::FHIter first = m_mesh->fh_begin(face);
61 m_mesh->calc_edge_vector(first, firstVec);
62 m_mesh->calc_edge_vector(++first, secondVec);
63 normal = cross(firstVec, secondVec);
64 normal.normalize();
65 m_mesh->set_normal(face, normal);
66 }
67 m_mesh->update_vertex_normals();
68 }
69
70 // precomp if computing zdir or normals
71 else if (zDir == ZDIR_Z || zDir == ZDIR_NORMAL)
72 {
73 m_zDir = zDir;
74 if (!m_mesh->has_face_colors() || !m_mesh->has_vertex_normals())
75 {
76 m_mesh->request_face_normals();
77 m_mesh->request_vertex_normals();
78 m_mesh->update_normals();
79 }
80 }
81
82 // if length is set as relative, we must compute median and multiply it with length
83 std::vector<ScalarT> all;
84 MeshT::EdgeIter ende = m_mesh->edges_end();
85
86 for (MeshT::EdgeIter edge = m_mesh->edges_begin(); edge != ende; ++edge)
87 all.push_back(m_mesh->calc_edge_length(edge));
88
89 std::sort(all.begin(), all.end());
90 m_length = all[all.size()/2];
91
92 if (lengthRelative)
93 {
94 length = all[all.size()/2] * length;
95 }
96
97 // pass all vertices and compute tangent matrices
98 MeshT::VertexIter end = m_mesh->vertices_end();
99 unsigned int i = 0;
100 unsigned int total = m_mesh->n_vertices();
101
102 /*std::ofstream file;
103 file.open(filename, std::fstream::app);*/
104
105 int N = m_mesh->n_vertices();
106 switch (xdirection)
107 {
108 case XDIR_CURVATURE:
109 #pragma omp parallel// shared(m_mesh)
110 {
111 #pragma omp for schedule(dynamic)
112 for (int i = 0; i < N; ++i)
113 {
114 MeshT::VertexHandle vertex = m_mesh->vertex_handle(i);
115 if (i%5000 == 0) MDS_LOG_NOTE("Computed " << i << " vertices of " << total << ".");
116 direction = m_mesh->curvature(vertex);
117 // create a transformation MatrixT - vertex normal specifies Z direction, curvature direction specifies X direction
118 OMTransformationSolver<PointT> solver(length, resolution, m_mesh->normal(vertex), direction, m_mesh->point(vertex));
119 rasterizeVertex(solver, vertex);
120 }
121 }
122 break;
123 case XDIR_DO_NOT_CHECK:
124 #pragma omp parallel //shared(m_mesh)
125 {
126 #pragma omp for schedule(dynamic)
127 for (int i = 0; i < N; ++i)
128 //for (MeshT::VertexIter vertex = m_mesh->vertices_begin(); vertex != end; ++vertex)
129 {
130 MeshT::VertexHandle vertex = m_mesh->vertex_handle(i);
131 if (i%5000 == 0) MDS_LOG_NOTE("Computed " << i << " vertices of " << total << ".");
132 //++i;
133 // create a transformation MatrixT - vertex normal specifies Z direction, curvature direction specifies X direction
134 OMTransformationSolver<PointT> solver(length, resolution, m_mesh->normal(vertex), PointT(0.0, 0.0, 0.0), m_mesh->point(vertex));
135 rasterizeVertex(solver, vertex);
136 }
137 }
138 break;
139
140 default:
141 return false;
142 }
143 return true;
144 }
145
147 // Function rasterizes one triangle given by face handle on a MatrixT with use of transformation MatrixT
148 // @param solver OMTransformationSolver with saved transformation MatrixT
149 // @param face Face handle of a triangle
150 // @param MatrixT MatrixT on which we will rasterize a face
152 template <class Mesh, class Matrix>
153 void OMProjector<Mesh, Matrix>::rasterizeOnMatrixT(OMTransformationSolver<PointT> &solver, FaceHandleT face, MatrixT &matrix)
154 {
155 PointT vertices[3];
156 PointT current;
157 // compute face normal for face equation
158 VectorT normal = solver.transformTo2DLinear(m_mesh->normal(face));
159 normal.normalize_cond();
160
161 // find all tri vertices
162 unsigned char i = 0;
163
164 for (MeshT::FVIter it = m_mesh->fv_begin(face); it; ++it, ++i)
165 {
166 vertices[i] = solver.transformTo2D(m_mesh->point(it));
167 //std::cout << vertices[i][0] << " " << vertices[i][1] << " " << vertices[i][2] << std::endl;
168 }
169
171 VectorT realNormal= (vertices[1]-vertices[0]) % (vertices[2]-vertices[1]);
172 //std::cout << "Normal " << realNormal[0] << " " << realNormal[1] << " " << realNormal[2] << std::endl;
173 if (realNormal[2] < 0)
174 {
175 current = vertices[2];
176 vertices[2] = vertices[0];
177 vertices[0] = current;
178 }
180 if (m_zDir == ZDIR_CURVATURE)
181 {
182 int i = 0;
183 for (MeshT::FVIter it = m_mesh->fv_begin(face); it; ++it, ++i)
184 {
185 vertices[i][2] = m_mesh->curvatureMagnitude(it);
186 ++i;
187 }
188 normal = (vertices[1] - vertices[0]) % (vertices[2] - vertices[0]);
189 normal.normalize();
190 }
191
192 ScalarT equationFract = normal[2] * solver.getPixelSize();
193 if (equationFract == 0.0) equationFract = 1.0;
194 // compute D coefficient of face equation (N1 * x + N2 * y + N3 * z + d = 0)
195 ScalarT d = -(normal[0]*vertices[0][0] + normal[1]*vertices[0][1] + normal[2]*vertices[0][2]);
196
197 // Use the triangle iterator and rasterize this triangle
198 OMToolkit::OMTriIterator<PointT> triangle(vertices);
199
200
201
202 if (m_zDir == ZDIR_NORMAL)
203 {
204 VectorT Z(0.0, 0.0, 1.0);
205 while(!triangle.isEnd())
206 {
207 current[0] = triangle.getX();
208 current[1] = triangle.getY();
209
210 if (current[0] >= m_MatrixTDelim.minimum[0] && current[0] <= m_MatrixTDelim.maximum[0] &&
211 current[1] >= m_MatrixTDelim.minimum[1] && current[1] <= m_MatrixTDelim.maximum[1])// && matrix(current[1], current[0]) > (sin(acos(normal | Z))))
212 matrix(current[1], current[0]) = (sin(acos(normal | Z)));
213 ++triangle;
214 }
215 }
216 else
217 while(!triangle.isEnd())
218 {
219 current[0] = triangle.getX();
220 current[1] = triangle.getY();
221 //std::cout << current[0] << " " << current[1] << std::endl;
222 if (current[0] >= m_MatrixTDelim.minimum[0] && current[0] <= m_MatrixTDelim.maximum[0] &&
223 current[1] >= m_MatrixTDelim.minimum[1] && current[1] <= m_MatrixTDelim.maximum[1])// && abs(matrix(current[1], current[0])) > abs((normal[0] * current[0] + normal[1] * current[1] + d)/equationFract))
224 matrix(current[1], current[0]) = (normal[0] * current[0] + normal[1] * current[1] + d)/equationFract;
225 //MDS_LOG_NOTE((normal[0] * current[0] + normal[1] * current[1] + d)/equationFract);
226 //MDS_LOG_NOTE(equationFract);
227 ++triangle;
228 }
229 }
230
231
233 // Function rasterizes neighbourhood of a vertex into a raster given by transformation solver
234 // Result is saved in MatrixT property of a vertex
235 // @param solver Transformation solver with saved transformation MatrixT
236 // @param vertex Vertex on which we compute a tangent raster
238 template <class Mesh, class Matrix>
239 void OMProjector<Mesh, Matrix>::rasterizeVertex(OMTransformationSolver<PointT> &solver, VertexHandleT vertex)
240 {
241 int INITIAL;
242
243 if (m_zDir == ZDIR_Z)
244 INITIAL = m_length;
245 else
246 INITIAL = M_PI;
247
248 std::set<int> addedFaces;
249 std::set<int> searchedVertices;
250 VertexHandleT originalHandle = vertex;
251 PointT currentPoint;
252 PointT transformed;
253 int index;
254
255 // create a new instance of MatrixT to be computed
256 MatrixT matrix(m_resolution, m_resolution);
257 matrix.setConstant(INITIAL);
258
259 m_MatrixTDelim.minimum = PointT(0.0, 0.0, 0.0);
260 m_MatrixTDelim.maximum = solver.getMaxBounds();
261
262 // stack to simulate recursive seed growth
263 std::vector<VertexHandleT> stack;
264 stack.push_back(vertex);
265
266 std::vector<FaceHandleT> facesToRasterize;
268 while (!stack.empty())
269 {
270 vertex = stack.back();
271 stack.pop_back();
272
273 currentPoint = m_mesh->point(vertex);
274
275 // rasterize all faces around current vertex
276 for (MeshT::VFIter around = m_mesh->vf_begin(vertex); around; ++around)
277 {
278 // get face index and test, if it wasn't rasterized yet
279 index = around.handle().idx();
280 if (addedFaces.count(index) == 0)
281 {
282 // if not, save its index and rasterize it
283 addedFaces.insert(index);
284 facesToRasterize.push_back(around);
285 //rasterizeOnMatrixT(solver, around, matrix);
286 }
287 }
288 // then look for all neighbouring vertices and test, if they fit in tangent raster
289 for (MeshT::VVIter around = m_mesh->vv_begin(vertex); around; ++around)
290 {
291 index = around.handle().idx();
292
293 // if it was not tested yet...
294 if (searchedVertices.count(index) == 0)
295 {
296 // transform it into raster coordinates
297 transformed = solver.transformTo2D(m_mesh->point(around));
298
299 // mark as searched
300 searchedVertices.insert(index);
301
302 // if it is in raster, save it for future search
303 if (transformed[0] >= m_MatrixTDelim.minimum[0] && transformed[0] <= m_MatrixTDelim.maximum[0] &&
304 transformed[1] >= m_MatrixTDelim.minimum[1] && transformed[1] <= m_MatrixTDelim.maximum[1])
305 {
306 stack.push_back(around);
307 for (MeshT::VVIter around2 = m_mesh->vv_begin(around); around2; ++around2)
308 {
309 index = around2.handle().idx();
310 if (searchedVertices.count(index) == 0)
311 {
312 stack.push_back(around2);
313 searchedVertices.insert(index);
314 }
315 }
316 }
317
318 }
319 }
320 }
321 unsigned int total = facesToRasterize.size();
322
323 // add all tris sharing boundary edges
324 for (unsigned int i = 0; i < total; ++i)
325 {
326 for (MeshT::FFIter around = m_mesh->ff_begin(facesToRasterize[i]); around; ++around)
327 facesToRasterize.push_back(around);
328 }
329
330 total = facesToRasterize.size();
331 for (unsigned int i = 0; i < total; ++i)
332 rasterizeOnMatrixT(solver, facesToRasterize[i], matrix);
333
334 // Repair holes (find by raycasting)
335 Mesh::Normal normal = m_mesh->normal(originalHandle);
336 normal.normalize_cond();
337
338 Mesh::FaceHandle face;
339
340 for(int x = 0; x < m_resolution; ++x)
341 for(int y = 0; y < m_resolution; ++y)
342 {
343 if (matrix(y, x) == INITIAL)
344 {
345 Mesh::Point point(x, y, 0.0);
346 point = solver.transformToMesh(point);
347
348 tree.getPassingFace(point, normal, face);
349
350 if (face.is_valid())
351 rasterizeOnMatrixT(solver, face, matrix);
352 }
353 }
354
355 m_mesh->property(m_propertyHandle, originalHandle) = matrix;
356 }
357
358 #endif _OM_PROJECTOR_HXX_