14 #ifndef _OM_CURVATURE_HXX_
15 #define _OM_CURVATURE_HXX_
27 template <
class Mesh,
class Scalar>
31 m_mesh->request_face_normals();
32 m_mesh->request_vertex_normals();
33 m_mesh->update_normals();
36 m_mesh->add_property(m_areas,
"<VoronoiA>");
37 Mesh::FIter end = m_mesh->faces_end();
38 for (Mesh::FIter face = m_mesh->faces_begin(); face!= end; ++face)
39 m_mesh->property(m_areas, face) =
AreasVector(0.0, 0.0, 0.0);
41 m_curvVectorHandle = curvVectorHandle;
42 m_curvMagHandle = curvMagHandle;
52 template <
class Mesh,
class Scalar>
56 Normal normal = m_mesh->normal(vertex);
72 normalVector(0, 0) = normal[0];
73 normalVector(1, 0) = normal[1];
74 normalVector(2, 0) = normal[2];
75 normalVectorT = normalVector.transpose();
81 identity.setIdentity();
87 normalMult = normalVector * normalVectorT;
90 Scalar sumArea = getVoronoiArea(vertex) * 2;
92 for (Mesh::VOHIter outgoing = m_mesh->voh_begin(vertex); outgoing; ++outgoing)
95 m_mesh->calc_edge_vector( outgoing, outVector);
97 edgeVectorIn(0, 0) = -outVector[0];
98 edgeVectorIn(1, 0) = -outVector[1];
99 edgeVectorIn(2, 0) = -outVector[2];
101 edgeVectorOut(0, 0) = outVector[0];
102 edgeVectorOut(1, 0) = outVector[1];
103 edgeVectorOut(2, 0) = outVector[2];
107 Tij = (identity - normalMult) * edgeVectorIn;
109 TijT = Tij.transpose();
112 vecLength = outVector.length();
113 Kij = ((normalVectorT * edgeVectorOut)(0) * 2) / (vecLength * vecLength);
116 face1 = m_mesh->face_handle(outgoing);
117 face2 = m_mesh->opposite_face_handle(outgoing);
120 if (m_mesh->is_valid_handle(face1))
122 if (isObtuse(face1, he))
123 w += getVoronoiAreaTriO(face1, vertex);
124 else w += getVoronoiAreaTri(face1, vertex);
126 if (m_mesh->is_valid_handle(face2))
128 if (isObtuse(face2, he))
129 w += getVoronoiAreaTriO(face2, vertex);
130 else w += getVoronoiAreaTri(face2, vertex);
134 M += w * Kij * Tij * TijT;
145 if ((E - normalVector).norm() > (E + normalVector).norm())
146 W = E - normalVector;
148 W = E + normalVector;
150 vecLength = W.norm();
157 result = Q.transpose() * M * Q;
162 subM(0, 0) = result(1, 1);
163 subM(0, 1) = result(1, 2);
164 subM(1, 0) = result(2, 1);
165 subM(1, 1) = result(2, 2);
168 Scalar r = sqrt(subM(0, 0) * subM(0, 0) + subM(1, 0)*subM(1, 0));
169 Scalar c = subM(0, 0)/r;
170 Scalar s = subM(1, 0)/r;
179 max = T1.operator*(c) - T2.operator*(s);
180 min = T1.operator*(s) + T2.operator*(c);
188 template <
class Mesh,
class Scalar>
191 if (!m_mesh->IsTriMesh)
return false;
196 Mesh::VertexIter end = m_mesh->vertices_end();
199 for (Mesh::VertexIter it = m_mesh->vertices_begin(); it != m_mesh->vertices_end(); ++it)
202 if (!m_mesh->is_boundary(it))
204 aux = getMeanOperator(it).length()/2;
205 delta = aux * aux - getGaussOperator(it);
206 m_mesh->property(m_curvMagHandle, it) = aux - sqrt(delta);
207 getCurvatureDirections(it, max, min);
209 m_mesh->property(m_curvVectorHandle, it) = min;
214 m_mesh->property(m_curvMagHandle, it) = 0.0;
215 m_mesh->property(m_curvVectorHandle, it) =
Normal(0.0, 0.0, 0.0);
226 template <
class Mesh,
class Scalar>
229 if (!m_mesh->IsTriMesh)
return false;
235 Mesh::VertexIter end = m_mesh->vertices_end();
237 for (Mesh::VertexIter it = m_mesh->vertices_begin(); it != m_mesh->vertices_end(); ++it)
240 if (!m_mesh->is_boundary(it))
242 aux = getMeanOperator(it).length()/2;
243 delta = aux * aux - getGaussOperator(it);
244 m_mesh->property(m_curvMagHandle, it) = aux + sqrt(delta);
245 getCurvatureDirections(it, max, min);
247 m_mesh->property(m_curvVectorHandle, it) = max;
252 m_mesh->property(m_curvMagHandle, it) = 0.0;
253 m_mesh->property(m_curvVectorHandle, it) =
Normal(0.0, 0.0, 0.0);
263 template <
class Mesh,
class Scalar>
266 if (!m_mesh->IsTriMesh)
return false;
268 Mesh::VertexIter end = m_mesh->vertices_end();
270 for (Mesh::VertexIter it = m_mesh->vertices_begin(); it != m_mesh->vertices_end(); ++it)
273 if (!m_mesh->is_boundary(it))
274 m_mesh->property(m_curvMagHandle, it) = getGaussOperator(it);
277 m_mesh->property(m_curvMagHandle, it) = 0.0;
280 m_mesh->property(m_curvVectorHandle, it) =
Normal(0.0, 0.0, 0.0);
289 template <
class Mesh,
class Scalar>
292 if (!m_mesh->IsTriMesh)
return false;
294 Mesh::VertexIter end = m_mesh->vertices_end();
296 for (Mesh::VertexIter it = m_mesh->vertices_begin(); it != m_mesh->vertices_end(); ++it)
299 if (!m_mesh->is_boundary(it))
301 Normal vector = getMeanOperator(it);
302 m_mesh->property(m_curvMagHandle, it) = vector.length() / 2.0;
304 m_mesh->property(m_curvVectorHandle, it) = vector;
309 m_mesh->property(m_curvMagHandle, it) = 0.0;
310 m_mesh->property(m_curvVectorHandle, it) =
Normal(0.0, 0.0, 0.0);
322 template <
class Mesh,
class Scalar>
327 for(it = m_mesh->fe_begin(face); it; ++it)
329 if (m_mesh->calc_sector_angle(it.current_halfedge_handle()) >
M_PIHALF)
331 half_edge = it.current_halfedge_handle();
344 template <
class Mesh,
class Scalar>
347 Mesh::FHIter it = m_mesh->fh_begin(face);
348 m_mesh->fv_begin(face);
353 if (m_mesh->to_vertex_handle(it) == vertex)
358 if (m_mesh->to_vertex_handle(it) != vertex)
return 0.0;
361 Scalar result = m_mesh->property(m_areas, face)[i];
362 if (result != 0)
return result;
369 Scalar angleIt = m_mesh->calc_sector_angle(it);
370 Scalar angleNext = m_mesh->calc_sector_angle(next);
395 Scalar lenght = m_mesh->calc_edge_length(it);
396 Scalar result = (lenght * lenght * tan(angleIt))/8.0;
397 m_mesh->property(m_areas, face)[i] = result;
403 Scalar lenght = m_mesh->calc_edge_length(next);
404 Scalar result = (lenght * lenght * tan(angleIt))/8.0;
405 m_mesh->property(m_areas, face)[i] = result;
412 Scalar result = m_mesh->calc_sector_area(it);
413 Scalar lenghtB = m_mesh->calc_edge_length(next);
414 result -= (lenghtB * lenghtB * tan(angleNext))/8.0;
416 Scalar lenghtC = m_mesh->calc_edge_length(it);
417 result -= ((lenghtC * lenghtC * tan(m_mesh->calc_sector_angle(prev)))/8.0);
418 m_mesh->property(m_areas, face)[i] = result;
496 template <
class Mesh,
class Scalar>
499 Mesh::FHIter it = m_mesh->fh_begin(face);
504 if (m_mesh->to_vertex_handle(it) == vertex)
508 if (m_mesh->to_vertex_handle(it) != vertex)
return 0.0;
511 Scalar result = m_mesh->property(m_areas, face)[i];
512 if (result != 0)
return result;
515 Mesh::FHIter hedge = it;
535 Scalar auxLenght = m_mesh->calc_edge_length(hedge);
536 Scalar angle = m_mesh->calc_sector_angle(next);
539 result = auxLenght * auxLenght * 1/tan(angle);
541 angle= m_mesh->calc_sector_angle(prev);
542 auxLenght = m_mesh->calc_edge_length(next);
544 result += auxLenght * auxLenght * 1/tan(angle);
547 m_mesh->property(m_areas, face)[i] = result;
556 template <
class Mesh,
class Scalar>
560 Normal result(0.0, 0.0, 0.0);
567 for (Mesh::VertexOHalfedgeIter it = m_mesh->voh_begin(vertex); it; ++it)
570 opp = m_mesh->opposite_halfedge_handle(it);
571 tang = tan(m_mesh->calc_sector_angle(m_mesh->next_halfedge_handle(it)));
574 tang = tan(m_mesh->calc_sector_angle(m_mesh->next_halfedge_handle(opp)));
578 m_mesh->calc_edge_vector(opp, vector);
584 aux = getVoronoiArea(vertex);
590 else return Normal(0.0, 0.0, 0.0);
598 template<
class Mesh,
class Scalar>
602 Normal result(0.0, 0.0, 0.0);
609 for (Mesh::VertexIHalfedgeIter it = m_mesh->vih_begin(vertex); it; ++it)
610 aux += m_mesh->calc_sector_angle(it);
613 area = getVoronoiArea(vertex);
627 template<
class Mesh,
class Scalar>
630 Mesh::VFIter it, end;
631 end = m_mesh->vf_end(vertex);
635 for (it = m_mesh->vf_begin(vertex); it; ++it)
637 if (!isObtuse(it, hedge))
639 result += getVoronoiAreaTri(it, vertex);
643 result += getVoronoiAreaTriO(it, vertex);
649 #endif // _OM_CURVATURE_HXX_