OMToolkit  1.0
The polygonal mesh processing tool.
 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>
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>
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_