more properties, small change to attributes

CMakeLists.txt

@@ -13,4 +13,5 @@ else()
     target_compile_options(polymesh PRIVATE -Wall -Werror -fPIC)
     target_compile_options(polymesh PRIVATE
         $<$<COMPILE_LANGUAGE:CXX>:-std=c++11>)
+    target_link_libraries(polymesh PUBLIC -fuse-ld=gold)
 endif()

src/polymesh/algorithms/properties.hh

@@ -39,21 +39,65 @@ bool is_isolated(edge_handle v);
 int valence(vertex_handle v);
 
 /// returns the area of the (flat) polygonal face
-template <class Vec3>
-typename field_3d<Vec3>::Scalar face_area(face_handle f, vertex_attribute<Vec3> const& position);
+template <class Vec3, class Scalar = typename field_3d<Vec3>::Scalar>
+Scalar face_area(face_handle f, vertex_attribute<Vec3> const& position);
 
 /// returns the center of gravity for a given (flat) polygonal face
 template <class Vec3>
 Vec3 face_centroid(face_handle f, vertex_attribute<Vec3> const& position);
 
-/// returns the area of a given triangle
+/// returns the (CCW) oriented face normal (assumes planar polygon)
 template <class Vec3>
-typename field_3d<Vec3>::Scalar triangle_area(face_handle f, vertex_attribute<Vec3> const& position);
+Vec3 face_normal(face_handle f, vertex_attribute<Vec3> const& position);
+
+/// returns the area of a given triangle
+template <class Vec3, class Scalar = typename field_3d<Vec3>::Scalar>
+Scalar triangle_area(face_handle f, vertex_attribute<Vec3> const& position);
 
 /// returns the center of gravity for a given triangle
 template <class Vec3>
 Vec3 triangle_centroid(face_handle f, vertex_attribute<Vec3> const& position);
 
+/// returns the (CCW) oriented face normal
+template <class Vec3>
+Vec3 triangle_normal(face_handle f, vertex_attribute<Vec3> const& position);
+
+/// returns a barycentric interpolation of a triangular face
+template <class Vec3>
+Vec3 bary_interpolate(face_handle f, Vec3 bary, vertex_attribute<Vec3> const& position);
+
+/// returns the length of an edge
+template <class Vec3, class Scalar = typename field_3d<Vec3>::Scalar>
+Scalar edge_length(edge_handle e, vertex_attribute<Vec3> const& position);
+
+/// returns the length of an edge
+template <class Vec3, class Scalar = typename field_3d<Vec3>::Scalar>
+Scalar edge_length(halfedge_handle h, vertex_attribute<Vec3> const& position);
+
+/// returns the (non-normalized) vector from -> to
+template <class Vec3>
+Vec3 edge_vector(halfedge_handle h, vertex_attribute<Vec3> const& position);
+
+/// returns the (normalized) vector from -> to (0 if from == to)
+template <class Vec3>
+Vec3 edge_dir(halfedge_handle h, vertex_attribute<Vec3> const& position);
+
+/// calculates the angle between this half-edge and the next one
+template <class Vec3, class Scalar = typename field_3d<Vec3>::Scalar>
+Scalar angle_to_next(halfedge_handle h, vertex_attribute<Vec3> const& position);
+
+/// calculates the angle between this half-edge and the previous one
+template <class Vec3, class Scalar = typename field_3d<Vec3>::Scalar>
+Scalar angle_to_prev(halfedge_handle h, vertex_attribute<Vec3> const& position);
+
+/// sum of face angles at this vertex
+template <class Vec3, class Scalar = typename field_3d<Vec3>::Scalar>
+Scalar angle_sum(vertex_handle v, vertex_attribute<Vec3> const& position);
+
+/// difference between 2pi and the angle sum (positive means less than 2pi)
+template <class Vec3, class Scalar = typename field_3d<Vec3>::Scalar>
+Scalar angle_defect(vertex_handle v, vertex_attribute<Vec3> const& position);
+
 /// ======== IMPLEMENTATION ========
 
 inline bool is_boundary(vertex_handle v) { return v.is_boundary(); }
@@ -70,8 +114,8 @@ inline bool is_isolated(edge_handle v) { return v.is_isolated(); }
 
 inline int valence(vertex_handle v) { return v.adjacent_vertices().size(); }
 
-template <class Vec3>
-typename field_3d<Vec3>::Scalar triangle_area(face_handle f, vertex_attribute<Vec3> const& position)
+template <class Vec3, class Scalar>
+Scalar triangle_area(face_handle f, vertex_attribute<Vec3> const& position)
 {
     auto h = f.any_halfedge();
     auto p0 = position[h.vertex_from()];
@@ -93,7 +137,31 @@ Vec3 triangle_centroid(face_handle f, vertex_attribute<Vec3> const& position)
 }
 
 template <class Vec3>
-typename field_3d<Vec3>::Scalar face_area(face_handle f, vertex_attribute<Vec3> const& position)
+Vec3 face_normal(face_handle f, vertex_attribute<Vec3> const& position)
+{
+    auto c = face_centroid(f, position);
+    auto e = f.any_halfedge();
+    auto v0 = e.vertex_from()[position];
+    auto v1 = e.vertex_to()[position];
+    auto n = field_3d<Vec3>::cross(v0 - c, v1 - c);
+    auto l = field_3d<Vec3>::length(n);
+    return l == 0 ? field_3d<Vec3>::zero() : n / l;
+}
+
+template <class Vec3>
+Vec3 triangle_normal(face_handle f, vertex_attribute<Vec3> const& position)
+{
+    auto e = f.any_halfedge();
+    auto v0 = e.vertex_from()[position];
+    auto v1 = e.vertex_to()[position];
+    auto v2 = e.next().vertex_to()[position];
+    auto n = field_3d<Vec3>::cross(v1 - v0, v2 - v0);
+    auto l = field_3d<Vec3>::length(n);
+    return l == 0 ? field_3d<Vec3>::zero() : n / l;
+}
+
+template <class Vec3, class Scalar>
+Scalar face_area(face_handle f, vertex_attribute<Vec3> const& position)
 {
     auto varea = field_3d<Vec3>::zero();
 
@@ -150,4 +218,98 @@ Vec3 face_centroid(face_handle f, vertex_attribute<Vec3> const& position)
 
     return centroid / (3.0f * area);
 }
+
+template <class Vec3, class Scalar>
+Scalar edge_length(edge_handle e, vertex_attribute<Vec3> const& position)
+{
+    return field_3d<Vec3>::length(position[e.vertexA()] - position[e.vertexB()]);
+}
+
+template <class Vec3, class Scalar>
+Scalar edge_length(halfedge_handle h, vertex_attribute<Vec3> const& position)
+{
+    return field_3d<Vec3>::length(position[h.vertex_from()] - position[h.vertex_to()]);
+}
+
+template <class Vec3>
+Vec3 edge_vector(halfedge_handle h, vertex_attribute<Vec3> const& position)
+{
+    return position[h.vertex_to()] - position[h.vertex_from()];
+}
+
+template <class Vec3>
+Vec3 edge_dir(halfedge_handle h, vertex_attribute<Vec3> const& position)
+{
+    auto d = position[h.vertex_to()] - position[h.vertex_from()];
+    auto l = field_3d<Vec3>::length(d);
+    if (l == 0)
+        return field_3d<Vec3>::zero();
+    return d / l;
+}
+
+template <class Vec3, class Scalar>
+Scalar angle_to_next(halfedge_handle h, vertex_attribute<Vec3> const& position)
+{
+    auto v0 = h.vertex_from()[position];
+    auto v1 = h.vertex_to()[position];
+    auto v2 = h.next().vertex_to()[position];
+
+    auto v01 = v0 - v1;
+    auto v21 = v2 - v1;
+
+    auto l01 = field_3d<Vec3>::length(v01);
+    auto l21 = field_3d<Vec3>::length(v21);
+
+    if (l01 == 0 || l21 == 0)
+        return 0;
+
+    auto ca = field_3d<Vec3>::dot(v01, v21) / (l01 * l21);
+    return std::acos(ca);
+}
+
+template <class Vec3, class Scalar>
+Scalar angle_to_prev(halfedge_handle h, vertex_attribute<Vec3> const& position)
+{
+    auto v0 = h.vertex_to()[position];
+    auto v1 = h.vertex_from()[position];
+    auto v2 = h.prev().vertex_from()[position];
+
+    auto v01 = v0 - v1;
+    auto v21 = v2 - v1;
+
+    auto l01 = field_3d<Vec3>::length(v01);
+    auto l21 = field_3d<Vec3>::length(v21);
+
+    if (l01 == 0 || l21 == 0)
+        return 0;
+
+    auto ca = field_3d<Vec3>::dot(v01, v21) / (l01 * l21);
+    return std::acos(ca);
+}
+
+template <class Vec3, class Scalar>
+Scalar angle_sum(vertex_handle v, vertex_attribute<Vec3> const& position)
+{
+    Scalar sum = 0;
+    for (auto h : v.outgoing_halfedges())
+        if (!h.is_boundary())
+            sum += angle_to_prev(h, position);
+    return sum;
+}
+
+template <class Vec3, class Scalar>
+Scalar angle_defect(vertex_handle v, vertex_attribute<Vec3> const& position)
+{
+    return 2 * M_PI - angle_sum(v, position);
+}
+
+template <class Vec3>
+Vec3 bary_interpolate(face_handle f, Vec3 bary, vertex_attribute<Vec3> const& position)
+{
+    auto h = f.any_halfedge();
+    auto v0 = h.vertex_to()[position];
+    auto v1 = h.next().vertex_to()[position];
+    auto v2 = h.next().next().vertex_to()[position];
+    return v0 * bary[0] + v1 * bary[1] + v2 * bary[2];
+}
 }

src/polymesh/attributes.hh

@@ -55,7 +55,7 @@ public:
     /// returns a new attribute where the given function was applied to each entry
     template <class FuncT>
     auto map(FuncT&& f) const -> attribute<tmp::decayed_result_type_of<FuncT, AttrT>>;
-    /// applies to given function to each attribute entry
+    /// applies to given function to each attribute entry (calls f(e))
     template <class FuncT>
     void apply(FuncT&& f);
 

src/polymesh/formats/ply.cc

+#include "ply.hh"
+
+// TODO

src/polymesh/formats/ply.hh

+#pragma once
+
+#include <glm/glm.hpp>
+
+#include <iostream>
+#include <string>
+
+#include "../Mesh.hh"
+
+namespace polymesh
+{
+// TODO
+// void write_ply(std::string const& filename, Mesh const& mesh, vertex_attribute<glm::vec3> const& position);
+// void write_ply(std::ostream& out, Mesh const& mesh, vertex_attribute<glm::vec3> const& position);
+// bool read_off(std::string const& filename, Mesh& mesh, vertex_attribute<glm::vec3>& position);
+// bool read_off(std::istream& input, Mesh& mesh, vertex_attribute<glm::vec3>& position);
+}

src/polymesh/impl/impl_attributes.hh

@@ -253,6 +253,6 @@ void primitive_attribute<tag, AttrT>::apply(FuncT &&f)
     auto s = size();
     auto d = data();
     for (auto i = 0; i < s; ++i)
-        d[i] = f(d[i]);
+        f(d[i]);
 }
 }