WIP: better API of feature generators + eigen analysis

This commit is contained in:
Simon Giraudot
2018-07-05 09:07:32 +02:00
parent 49aea9ec26
commit 3f37fa504f
10 changed files with 423 additions and 410 deletions
@@ -48,8 +48,6 @@ public:
std::size_t idx)
: eigen (eigen), m_idx (idx)
{
std::cerr << "[" << idx << "]" << std::endl;
std::cerr << "[" << m_idx << "]" << std::endl;
std::ostringstream oss;
oss << "eigenvalue" << (idx+1);
this->set_name (oss.str());
@@ -138,8 +138,6 @@ public:
* (c[std::size_t(channel)] - mean) / (2. * sd * sd)));
}
#endif
std::cerr << "[" << channel << " / " << mean << " / " << sd << "]" << std::endl;
std::cerr << "[" << m_channel << " / " << m_mean << " / " << m_sd << "]" << std::endl;
std::ostringstream oss;
if (channel == HUE) oss << "hue";
else if (channel == SATURATION) oss << "saturation";
@@ -80,6 +80,8 @@ public:
#ifdef CGAL_LINKED_WITH_TBB
if (m_tasks != NULL)
delete m_tasks;
for (std::size_t i = 0; i < m_adders.size(); ++ i)
delete m_adders[i];
#endif
}
/// \endcond
@@ -118,6 +120,31 @@ public:
return m_features.back();
}
/// \cond SKIP_IN_MANUAL
template <typename Feature, typename ... T>
Feature_handle add_with_scale_id (std::size_t i, T&& ... t)
{
#ifdef CGAL_LINKED_WITH_TBB
if (m_tasks != NULL)
{
m_features.push_back (Feature_handle());
Parallel_feature_adder<Feature, T...>* adder
= new Parallel_feature_adder<Feature, T...>(i, m_features.back(), std::forward<T>(t)...);
m_adders.push_back (adder);
m_tasks->run (*adder);
}
else
#endif
{
m_features.push_back (Feature_handle (new Feature(std::forward<T>(t)...)));
m_features.back()->set_name (m_features.back()->name() + "_" + std::to_string(i));
}
return m_features.back();
}
/// \end
#if defined(CGAL_LINKED_WITH_TBB) || defined(DOXYGEN_RUNNING)
void begin_parallel_additions()
{
@@ -132,6 +159,7 @@ public:
for (std::size_t i = 0; i < m_adders.size(); ++ i)
delete m_adders[i];
m_adders.clear();
}
#endif
@@ -202,11 +230,18 @@ private:
template <typename Feature, typename ... T>
struct Parallel_feature_adder : Abstract_parallel_feature_adder
{
Feature_handle fh;
std::size_t scale;
mutable Feature_handle fh;
boost::shared_ptr<std::tuple<T...> > args;
Parallel_feature_adder (Feature_handle fh, T&& ... t)
: fh (fh)
: scale (std::size_t(-1)), fh (fh)
{
args = boost::make_shared<std::tuple<T...> >(std::forward<T>(t)...);
}
Parallel_feature_adder (std::size_t scale, Feature_handle fh, T&& ... t)
: scale(scale), fh (fh)
{
args = boost::make_shared<std::tuple<T...> >(std::forward<T>(t)...);
}
@@ -230,6 +265,8 @@ private:
void add_feature (Tuple& t, seq<S...>) const
{
fh.attach (new Feature (std::forward<T>(std::get<S>(t))...));
if (scale != std::size_t(-1))
fh->set_name (fh->name() + "_" + std::to_string(scale));
}
void operator()() const
@@ -195,18 +195,22 @@ private:
};
typedef CGAL::cpp11::array<float, 3> float3;
std::vector<float3> m_eigenvalues;
std::vector<float> m_sum_eigenvalues;
std::vector<float3> m_centroids;
std::vector<float3> m_smallest_eigenvectors;
struct Content
{
std::vector<float3> eigenvalues;
std::vector<float> sum_eigenvalues;
std::vector<float3> centroids;
std::vector<float3> smallest_eigenvectors;
#ifdef CGAL_CLASSIFICATION_EIGEN_FULL_STORAGE
std::vector<float3> m_middle_eigenvectors;
std::vector<float3> m_largest_eigenvectors;
std::vector<float3> middle_eigenvectors;
std::vector<float3> largest_eigenvectors;
#endif
float m_mean_range;
float mean_range;
};
boost::shared_ptr<Content> m_content; // To avoid copies with named constructors
public:
@@ -214,9 +218,12 @@ public:
Local_eigen_analysis () { }
/// \endcond
/// \name Named Constructors
/// @{
/*!
\brief Computes the local eigen analysis of an input range based
on a local neighborhood.
\brief Computes the local eigen analysis of an input point set
based on a local neighborhood.
\tparam PointRange model of `ConstRange`. Its iterator type is
`RandomAccessIterator` and its value type is the key type of
@@ -252,22 +259,24 @@ public:
#else
typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<float, 3> >
#endif
Local_eigen_analysis (const PointRange& input,
PointMap point_map,
const NeighborQuery& neighbor_query,
const ConcurrencyTag& = ConcurrencyTag(),
const DiagonalizeTraits& = DiagonalizeTraits())
static Local_eigen_analysis create_from_point_set(const PointRange& input,
PointMap point_map,
const NeighborQuery& neighbor_query,
const ConcurrencyTag& = ConcurrencyTag(),
const DiagonalizeTraits& = DiagonalizeTraits())
{
m_eigenvalues.resize (input.size());
m_sum_eigenvalues.resize (input.size());
m_centroids.resize (input.size());
m_smallest_eigenvectors.resize (input.size());
Local_eigen_analysis out;
out.m_content = boost::make_shared<Content>();
out.m_content->eigenvalues.resize (input.size());
out.m_content->sum_eigenvalues.resize (input.size());
out.m_content->centroids.resize (input.size());
out.m_content->smallest_eigenvectors.resize (input.size());
#ifdef CGAL_CLASSIFICATION_EIGEN_FULL_STORAGE
m_middle_eigenvectors.resize (input.size());
m_largest_eigenvectors.resize (input.size());
out.m_content->middle_eigenvectors.resize (input.size());
out.m_content->largest_eigenvectors.resize (input.size());
#endif
m_mean_range = 0.;
out.m_content->mean_range = 0.;
#ifndef CGAL_LINKED_WITH_TBB
CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
@@ -277,7 +286,7 @@ public:
{
tbb::mutex mutex;
Compute_eigen_values<PointRange, PointMap, NeighborQuery, DiagonalizeTraits>
f(*this, input, point_map, neighbor_query, m_mean_range, mutex);
f(out, input, point_map, neighbor_query, out.m_content->mean_range, mutex);
tbb::parallel_for(tbb::blocked_range<size_t>(0, input.size ()), f);
}
else
@@ -292,22 +301,38 @@ public:
for (std::size_t j = 0; j < neighbors.size(); ++ j)
neighbor_points.push_back (get(point_map, *(input.begin()+neighbors[j])));
m_mean_range += float(CGAL::sqrt (CGAL::squared_distance
out.m_content->mean_range += float(CGAL::sqrt (CGAL::squared_distance
(get(point_map, *(input.begin() + i)),
get(point_map, *(input.begin() + neighbors.back())))));
compute<typename PointMap::value_type, DiagonalizeTraits>
out.compute<typename PointMap::value_type, DiagonalizeTraits>
(i, get(point_map, *(input.begin()+i)), neighbor_points);
}
}
m_mean_range /= input.size();
}
out.m_content->mean_range /= input.size();
// Experimental feature, not used officially
/// \cond SKIP_IN_MANUAL
struct Input_is_clusters { };
struct Input_is_face_graph { };
return out;
}
/*!
\brief Computes the local eigen analysis of an input face graph
based on a local neighborhood.
\tparam FaceListGraph model of `FaceListGraph`.
\tparam NeighborQuery model of `NeighborQuery`
\tparam ConcurrencyTag enables sequential versus parallel
algorithm. Possible values are `Parallel_tag` (default value is %CGAL
is linked with TBB) or `Sequential_tag` (default value otherwise).
\tparam DiagonalizeTraits model of `DiagonalizeTraits` used for
matrix diagonalization. It can be omitted: if Eigen 3 (or greater)
is available and `CGAL_EIGEN3_ENABLED` is defined then an overload
using `Eigen_diagonalize_traits` is provided. Otherwise, the
internal implementation `Diagonalize_traits` is used.
\param input face graph.
\param neighbor_query object used to access neighborhoods of points.
*/
template <typename FaceListGraph,
typename NeighborQuery,
#if defined(DOXYGEN_RUNNING)
@@ -322,29 +347,31 @@ public:
#else
typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<float, 3> >
#endif
Local_eigen_analysis (const Input_is_face_graph&,
const FaceListGraph& input,
const NeighborQuery& neighbor_query,
const ConcurrencyTag& = ConcurrencyTag(),
const DiagonalizeTraits& = DiagonalizeTraits())
static Local_eigen_analysis create_from_face_graph (const FaceListGraph& input,
const NeighborQuery& neighbor_query,
const ConcurrencyTag& = ConcurrencyTag(),
const DiagonalizeTraits& = DiagonalizeTraits())
{
typedef typename boost::graph_traits<FaceListGraph>::face_descriptor face_descriptor;
typedef typename boost::graph_traits<FaceListGraph>::face_iterator face_iterator;
typedef typename CGAL::Iterator_range<face_iterator> Face_range;
typedef typename boost::property_map<FaceListGraph, CGAL::face_index_t>::type::value_type face_index;
Local_eigen_analysis out;
out.m_content = boost::make_shared<Content>();
Face_range range (faces(input));
m_eigenvalues.resize (range.size());
m_sum_eigenvalues.resize (range.size());
m_centroids.resize (range.size());
m_smallest_eigenvectors.resize (range.size());
out.m_content->eigenvalues.resize (range.size());
out.m_content->sum_eigenvalues.resize (range.size());
out.m_content->centroids.resize (range.size());
out.m_content->smallest_eigenvectors.resize (range.size());
#ifdef CGAL_CLASSIFICATION_EIGEN_FULL_STORAGE
m_middle_eigenvectors.resize (range.size());
m_largest_eigenvectors.resize (range.size());
out.m_content->middle_eigenvectors.resize (range.size());
out.m_content->largest_eigenvectors.resize (range.size());
#endif
m_mean_range = 0.;
out.m_content->mean_range = 0.;
#ifndef CGAL_LINKED_WITH_TBB
CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
@@ -354,7 +381,7 @@ public:
{
tbb::mutex mutex;
Compute_eigen_values_graph<FaceListGraph, NeighborQuery, DiagonalizeTraits>
f(*this, input, neighbor_query, m_mean_range, mutex);
f(out, input, neighbor_query, out.m_content->mean_range, mutex);
tbb::parallel_for(tbb::blocked_range<std::size_t>(0, range.size()), f);
}
@@ -366,17 +393,41 @@ public:
std::vector<face_index> neighbors;
neighbor_query (fd, std::back_inserter (neighbors));
m_mean_range += face_radius(fd, input);
out.m_content->mean_range += out.face_radius(fd, input);
compute_triangles<FaceListGraph, DiagonalizeTraits>
out.compute_triangles<FaceListGraph, DiagonalizeTraits>
(input, fd, neighbors);
}
}
m_mean_range /= range.size();
out.m_content->mean_range /= range.size();
return out;
}
/*!
\brief Computes the local eigen analysis of an input set of point
clusters based on a local neighborhood.
\tparam ClusterRange model of `ConstRange`. Its iterator type is
`RandomAccessIterator` and its value type is the key type of
`PointMap`.
\tparam PointMap model of `ReadablePropertyMap` whose key
type is the value type of the iterator of `PointRange` and value type
is `CGAL::Point_3`.
\tparam ConcurrencyTag enables sequential versus parallel
algorithm. Possible values are `Parallel_tag` (default value is %CGAL
is linked with TBB) or `Sequential_tag` (default value otherwise).
\tparam DiagonalizeTraits model of `DiagonalizeTraits` used for
matrix diagonalization. It can be omitted: if Eigen 3 (or greater)
is available and `CGAL_EIGEN3_ENABLED` is defined then an overload
using `Eigen_diagonalize_traits` is provided. Otherwise, the
internal implementation `Diagonalize_traits` is used.
\param input point range.
\param point_map property map to access the input points.
\param neighbor_query object used to access neighborhoods of points.
*/
template <typename ClusterRange,
typename PointMap,
#if defined(DOXYGEN_RUNNING)
@@ -391,25 +442,27 @@ public:
#else
typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<float, 3> >
#endif
Local_eigen_analysis (const Input_is_clusters&,
const ClusterRange& input,
PointMap point_map,
const ConcurrencyTag& = ConcurrencyTag(),
const DiagonalizeTraits& = DiagonalizeTraits())
static Local_eigen_analysis create_from_point_clusters (const ClusterRange& input,
PointMap point_map,
const ConcurrencyTag& = ConcurrencyTag(),
const DiagonalizeTraits& = DiagonalizeTraits())
{
m_eigenvalues.resize (input.size());
m_sum_eigenvalues.resize (input.size());
m_centroids.resize (input.size());
m_smallest_eigenvectors.resize (input.size());
Local_eigen_analysis out;
out.m_content = boost::make_shared<Content>();
out.m_content->eigenvalues.resize (input.size());
out.m_content->sum_eigenvalues.resize (input.size());
out.m_content->centroids.resize (input.size());
out.m_content->smallest_eigenvectors.resize (input.size());
#ifdef CGAL_CLASSIFICATION_EIGEN_FULL_STORAGE
m_middle_eigenvectors.resize (input.size());
m_largest_eigenvectors.resize (input.size());
out.m_content->middle_eigenvectors.resize (input.size());
out.m_content->largest_eigenvectors.resize (input.size());
#endif
m_mean_range = 0.;
out.m_content->mean_range = 0.;
Compute_clusters_eigen_values<ClusterRange, PointMap, DiagonalizeTraits>
f(*this, input, point_map);
f(out, input, point_map);
#ifndef CGAL_LINKED_WITH_TBB
@@ -426,8 +479,14 @@ public:
for (std::size_t i = 0; i < input.size(); ++ i)
f.apply (i);
}
return out;
}
/// \endcond
/// @}
/// \name Analysis
/// @{
/*!
\brief Returns the estimated unoriented normal vector of the point at position `index`.
@@ -436,9 +495,9 @@ public:
template <typename GeomTraits>
typename GeomTraits::Vector_3 normal_vector (std::size_t index) const
{
return typename GeomTraits::Vector_3(double(m_smallest_eigenvectors[index][0]),
double(m_smallest_eigenvectors[index][1]),
double(m_smallest_eigenvectors[index][2]));
return typename GeomTraits::Vector_3(double(m_content->smallest_eigenvectors[index][0]),
double(m_content->smallest_eigenvectors[index][1]),
double(m_content->smallest_eigenvectors[index][2]));
}
/*!
@@ -449,26 +508,28 @@ public:
typename GeomTraits::Plane_3 plane (std::size_t index) const
{
return typename GeomTraits::Plane_3
(typename GeomTraits::Point_3 (double(m_centroids[index][0]),
double(m_centroids[index][1]),
double(m_centroids[index][2])),
typename GeomTraits::Vector_3 (double(m_smallest_eigenvectors[index][0]),
double(m_smallest_eigenvectors[index][1]),
double(m_smallest_eigenvectors[index][2])));
(typename GeomTraits::Point_3 (double(m_content->centroids[index][0]),
double(m_content->centroids[index][1]),
double(m_content->centroids[index][2])),
typename GeomTraits::Vector_3 (double(m_content->smallest_eigenvectors[index][0]),
double(m_content->smallest_eigenvectors[index][1]),
double(m_content->smallest_eigenvectors[index][2])));
}
/*!
\brief Returns the normalized eigenvalues of the point at position `index`.
*/
const Eigenvalues& eigenvalue (std::size_t index) const { return m_eigenvalues[index]; }
const Eigenvalues& eigenvalue (std::size_t index) const { return m_content->eigenvalues[index]; }
/*!
\brief Returns the sum of eigenvalues of the point at position `index`.
*/
float sum_of_eigenvalues (std::size_t index) const { return m_sum_eigenvalues[index]; }
float sum_of_eigenvalues (std::size_t index) const { return m_content->sum_eigenvalues[index]; }
/// @}
/// \cond SKIP_IN_MANUAL
float mean_range() const { return m_mean_range; }
float mean_range() const { return m_content->mean_range; }
/// \endcond
private:
@@ -497,18 +558,18 @@ private:
if (neighbor_points.size() == 0)
{
Eigenvalues v = make_array( 0.f, 0.f, 0.f );
m_eigenvalues[index] = v;
m_centroids[index] = make_array(float(query.x()), float(query.y()), float(query.z()) );
m_smallest_eigenvectors[index] = make_array( 0.f, 0.f, 1.f );
m_content->eigenvalues[index] = v;
m_content->centroids[index] = make_array(float(query.x()), float(query.y()), float(query.z()) );
m_content->smallest_eigenvectors[index] = make_array( 0.f, 0.f, 1.f );
#ifdef CGAL_CLASSIFICATION_EIGEN_FULL_STORAGE
m_middle_eigenvectors[index] = make_array( 0.f, 1.f, 0.f );
m_largest_eigenvectors[index] = make_array( 1.f, 0.f, 0.f );
m_content->middle_eigenvectors[index] = make_array( 0.f, 1.f, 0.f );
m_content->largest_eigenvectors[index] = make_array( 1.f, 0.f, 0.f );
#endif
return;
}
Point centroid = CGAL::centroid (neighbor_points.begin(), neighbor_points.end());
m_centroids[index] = make_array( float(centroid.x()), float(centroid.y()), float(centroid.z()) );
m_content->centroids[index] = make_array( float(centroid.x()), float(centroid.y()), float(centroid.z()) );
CGAL::cpp11::array<float, 6> covariance = make_array( 0.f, 0.f, 0.f, 0.f, 0.f, 0.f );
@@ -536,12 +597,12 @@ private:
if (sum > 0.f)
for (std::size_t i = 0; i < 3; ++ i)
evalues[i] = evalues[i] / sum;
m_sum_eigenvalues[index] = float(sum);
m_eigenvalues[index] = make_array( float(evalues[0]), float(evalues[1]), float(evalues[2]) );
m_smallest_eigenvectors[index] = make_array( float(evectors[0]), float(evectors[1]), float(evectors[2]) );
m_content->sum_eigenvalues[index] = float(sum);
m_content->eigenvalues[index] = make_array( float(evalues[0]), float(evalues[1]), float(evalues[2]) );
m_content->smallest_eigenvectors[index] = make_array( float(evectors[0]), float(evectors[1]), float(evectors[2]) );
#ifdef CGAL_CLASSIFICATION_EIGEN_FULL_STORAGE
m_middle_eigenvectors[index] = make_array( float(evectors[3]), float(evectors[4]), float(evectors[5]) );
m_largest_eigenvectors[index] = make_array( float(evectors[6]), float(evectors[7]), float(evectors[8]) );
m_content->middle_eigenvectors[index] = make_array( float(evectors[3]), float(evectors[4]), float(evectors[5]) );
m_content->largest_eigenvectors[index] = make_array( float(evectors[6]), float(evectors[7]), float(evectors[8]) );
#endif
}
@@ -560,7 +621,7 @@ private:
if (neighbor_faces.size() == 0)
{
Eigenvalues v = {{ 0.f, 0.f, 0.f }};
m_eigenvalues[get(get(CGAL::face_index,g), query)] = v;
m_content->eigenvalues[get(get(CGAL::face_index,g), query)] = v;
CGAL::cpp11::array<Triangle,1> tr
= {{ Triangle (get(get (CGAL::vertex_point, g), target(halfedge(query, g), g)),
@@ -569,12 +630,12 @@ private:
Point c = CGAL::centroid(tr.begin(),
tr.end(), Kernel(), CGAL::Dimension_tag<2>());
m_centroids[get(get(CGAL::face_index,g), query)] = {{ float(c.x()), float(c.y()), float(c.z()) }};
m_content->centroids[get(get(CGAL::face_index,g), query)] = {{ float(c.x()), float(c.y()), float(c.z()) }};
m_smallest_eigenvectors[get(get(CGAL::face_index,g), query)] = {{ 0.f, 0.f, 1.f }};
m_content->smallest_eigenvectors[get(get(CGAL::face_index,g), query)] = {{ 0.f, 0.f, 1.f }};
#ifdef CGAL_CLASSIFICATION_EIGEN_FULL_STORAGE
m_middle_eigenvectors[get(get(CGAL::face_index,g), query)] = {{ 0.f, 1.f, 0.f }};
m_largest_eigenvectors[get(get(CGAL::face_index,g), query)] = {{ 1.f, 0.f, 0.f }};
m_content->middle_eigenvectors[get(get(CGAL::face_index,g), query)] = {{ 0.f, 1.f, 0.f }};
m_content->largest_eigenvectors[get(get(CGAL::face_index,g), query)] = {{ 1.f, 0.f, 0.f }};
#endif
return;
}
@@ -600,7 +661,7 @@ private:
c, Kernel(), (Triangle*)NULL, CGAL::Dimension_tag<2>(),
DiagonalizeTraits());
m_centroids[get(get(CGAL::face_index,g), query)] = {{ float(c.x()), float(c.y()), float(c.z()) }};
m_content->centroids[get(get(CGAL::face_index,g), query)] = {{ float(c.x()), float(c.y()), float(c.z()) }};
CGAL::cpp11::array<double, 3> evalues = {{ 0.f, 0.f, 0.f }};
CGAL::cpp11::array<double, 9> evectors = {{ 0.f, 0.f, 0.f,
@@ -615,12 +676,12 @@ private:
if (sum > 0.f)
for (std::size_t i = 0; i < 3; ++ i)
evalues[i] = evalues[i] / sum;
m_sum_eigenvalues[get(get(CGAL::face_index,g), query)] = float(sum);
m_eigenvalues[get(get(CGAL::face_index,g), query)] = {{ float(evalues[0]), float(evalues[1]), float(evalues[2]) }};
m_smallest_eigenvectors[get(get(CGAL::face_index,g), query)] = {{ float(evectors[0]), float(evectors[1]), float(evectors[2]) }};
m_content->sum_eigenvalues[get(get(CGAL::face_index,g), query)] = float(sum);
m_content->eigenvalues[get(get(CGAL::face_index,g), query)] = {{ float(evalues[0]), float(evalues[1]), float(evalues[2]) }};
m_content->smallest_eigenvectors[get(get(CGAL::face_index,g), query)] = {{ float(evectors[0]), float(evectors[1]), float(evectors[2]) }};
#ifdef CGAL_CLASSIFICATION_EIGEN_FULL_STORAGE
m_middle_eigenvectors[get(get(CGAL::face_index,g), query)] = {{ float(evectors[3]), float(evectors[4]), float(evectors[5]) }};
m_largest_eigenvectors[get(get(CGAL::face_index,g), query)] = {{ float(evectors[6]), float(evectors[7]), float(evectors[8]) }};
m_content->middle_eigenvectors[get(get(CGAL::face_index,g), query)] = {{ float(evectors[3]), float(evectors[4]), float(evectors[5]) }};
m_content->largest_eigenvectors[get(get(CGAL::face_index,g), query)] = {{ float(evectors[6]), float(evectors[7]), float(evectors[8]) }};
#endif
}
@@ -152,9 +152,10 @@ private:
t.reset();
t.start();
eigen = new Local_eigen_analysis (typename Local_eigen_analysis::Input_is_face_graph(),
input, neighborhood->n_ring_neighbor_query(nb_scale + 1),
ConcurrencyTag(), DiagonalizeTraits());
eigen = new Local_eigen_analysis
(Local_eigen_analysis::create_from_face_graph
(input, neighborhood->n_ring_neighbor_query(nb_scale + 1),
ConcurrencyTag(), DiagonalizeTraits()));
float mrange = eigen->mean_range();
if (this->voxel_size < 0)
this->voxel_size = mrange;
@@ -216,15 +217,29 @@ public:
Mesh_feature_generator(Feature_set& features,
const FaceListGraph& input,
PointMap point_map,
std::size_t nb_scales)
std::size_t nb_scales,
float voxel_size = -1.f)
: m_input (input), m_range(faces(input)), m_point_map (point_map), m_features (features)
{
m_bbox = CGAL::bounding_box
(boost::make_transform_iterator (m_range.begin(), CGAL::Property_map_to_unary_function<PointMap>(m_point_map)),
boost::make_transform_iterator (m_range.end(), CGAL::Property_map_to_unary_function<PointMap>(m_point_map)));
generate_features_impl (nb_scales);
CGAL::Real_timer t; t.start();
m_scales.reserve (nb_scales);
m_scales.push_back (new Scale (m_input, m_range, m_point_map, m_bbox, voxel_size, 0));
voxel_size = m_scales[0]->grid_resolution();
for (std::size_t i = 1; i < nb_scales; ++ i)
{
voxel_size *= 2;
m_scales.push_back (new Scale (m_input, m_range, m_point_map, m_bbox, voxel_size, i, m_scales[i-1]->grid));
}
t.stop();
CGAL_CLASSIFICATION_CERR << "Scales computed in " << t.time() << " second(s)" << std::endl;
t.reset();
}
@@ -244,6 +259,42 @@ public:
/// \endcond
void generate_point_based_features ()
{
#ifdef DO_NOT_USE_EIGEN_FEATURES
for (int j = 0; j < 3; ++ j)
{
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Eigenvalue> (i, m_range, eigen(i), std::size_t(j));
}
#else
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Anisotropy> (i, i, m_range, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Eigentropy> (i, m_range, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Linearity> (i, m_range, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Omnivariance> (i, m_range, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Planarity> (i, m_range, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Sphericity> (i, m_range, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Sum_eigen> (i, m_range, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Surface_variation> (i, m_range, eigen(i));
#endif
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Distance_to_plane> (i, m_range, m_point_map, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Dispersion> (i, m_range, m_point_map, grid(i), radius_neighbors(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Elevation> (i, m_range, m_point_map, grid(i), radius_dtm(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Verticality> (i, m_range, eigen(i));
}
/*!
\brief Returns the bounding box of the input point set.
*/
@@ -297,51 +348,6 @@ private:
m_scales.clear();
}
void add (Feature_handle fh, std::size_t i)
{
m_features_to_rename.push_back (std::make_pair (fh, i));
}
void generate_point_based_features ()
{
#ifdef DO_NOT_USE_EIGEN_FEATURES
for (int j = 0; j < 3; ++ j)
{
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Eigenvalue> (m_range, eigen(i), std::size_t(j)), i);
}
#else
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Anisotropy> (m_range, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Eigentropy> (m_range, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Linearity> (m_range, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Omnivariance> (m_range, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Planarity> (m_range, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Sphericity> (m_range, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Sum_eigen> (m_range, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Surface_variation> (m_range, eigen(i)), i);
#endif
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Distance_to_plane> (m_range, m_point_map, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Dispersion> (m_range, m_point_map, grid(i), radius_neighbors(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Elevation> (m_range, m_point_map, grid(i), radius_dtm(i)), i);
}
void generate_normal_based_features(const CGAL::Default_property_map<face_iterator, typename Geom_traits::Vector_3>&)
{
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Verticality> (m_range, eigen(i)), i);
}
template <typename T>
const T& get_parameter (const T& t)
{
@@ -355,46 +361,6 @@ private:
return Default_property_map<face_iterator, T>();
}
void generate_features_impl (std::size_t nb_scales)
{
CGAL::Real_timer t; t.start();
m_scales.reserve (nb_scales);
float voxel_size = -1;
m_scales.push_back (new Scale (m_input, m_range, m_point_map, m_bbox, voxel_size, 0));
voxel_size = m_scales[0]->grid_resolution();
for (std::size_t i = 1; i < nb_scales; ++ i)
{
voxel_size *= 2;
m_scales.push_back (new Scale (m_input, m_range, m_point_map, m_bbox, voxel_size, i, m_scales[i-1]->grid));
}
t.stop();
CGAL_CLASSIFICATION_CERR << "Scales computed in " << t.time() << " second(s)" << std::endl;
t.reset();
t.start();
#ifdef CGAL_LINKED_WITH_TBB
m_features.begin_parallel_additions();
#endif
generate_point_based_features ();
generate_normal_based_features (CGAL::Default_property_map<face_iterator, typename Geom_traits::Vector_3>());
#ifdef CGAL_LINKED_WITH_TBB
m_features.end_parallel_additions();
#endif
for (std::size_t i = 0; i < m_features_to_rename.size(); ++ i)
m_features_to_rename[i].first->set_name
(m_features_to_rename[i].first->name() + "_" + std::to_string(m_features_to_rename[i].second));
t.stop();
CGAL_CLASSIFICATION_CERR << "Features computed in " << t.time() << " second(s)" << std::endl;
}
};
@@ -171,7 +171,7 @@ private:
{
CGAL::Real_timer t;
t.start();
if (voxel_size < 0.)
if (lower_grid == NULL)
neighborhood = new Neighborhood (input, point_map);
else
neighborhood = new Neighborhood (input, point_map, voxel_size);
@@ -186,7 +186,8 @@ private:
t.start();
eigen = new Local_eigen_analysis
(input, point_map, neighborhood->k_neighbor_query(12), ConcurrencyTag(), DiagonalizeTraits());
(Local_eigen_analysis::create_from_point_set
(input, point_map, neighborhood->k_neighbor_query(12), ConcurrencyTag(), DiagonalizeTraits()));
float range = eigen->mean_range();
if (this->voxel_size < 0)
@@ -237,7 +238,6 @@ private:
const PointRange& m_input;
PointMap m_point_map;
Feature_set& m_features;
std::vector<std::pair<Feature_handle, std::size_t> > m_features_to_rename;
public:
@@ -309,37 +309,34 @@ public:
\param color_map property map to access the colors of the input points (if any).
\param echo_map property map to access the echo values of the input points (if any).
*/
template <typename VectorMap = Default,
typename ColorMap = Default,
typename EchoMap = Default>
Point_set_feature_generator(Feature_set& features,
const PointRange& input,
PointMap point_map,
std::size_t nb_scales,
VectorMap normal_map = VectorMap(),
ColorMap color_map = ColorMap(),
EchoMap echo_map = EchoMap()
#ifndef DOXYGEN_RUNNING
, float voxel_size = -1.f // Undocumented way of changing base voxel size
#endif
)
float voxel_size = -1.f)
: m_input (input), m_point_map (point_map), m_features (features)
{
m_bbox = CGAL::bounding_box
(boost::make_transform_iterator (m_input.begin(), CGAL::Property_map_to_unary_function<PointMap>(m_point_map)),
boost::make_transform_iterator (m_input.end(), CGAL::Property_map_to_unary_function<PointMap>(m_point_map)));
typedef typename Default::Get<VectorMap, typename GeomTraits::Vector_3 >::type
Vmap;
typedef typename Default::Get<ColorMap, RGB_Color >::type
Cmap;
typedef typename Default::Get<EchoMap, std::size_t >::type
Emap;
CGAL::Real_timer t; t.start();
m_scales.reserve (nb_scales);
m_scales.push_back (new Scale (m_input, m_point_map, m_bbox, voxel_size));
generate_features_impl (nb_scales, voxel_size,
get_parameter<Vmap>(normal_map),
get_parameter<Cmap>(color_map),
get_parameter<Emap>(echo_map));
if (voxel_size == -1.f)
voxel_size = m_scales[0]->grid_resolution();
for (std::size_t i = 1; i < nb_scales; ++ i)
{
voxel_size *= 2;
m_scales.push_back (new Scale (m_input, m_point_map, m_bbox, voxel_size, m_scales[i-1]->grid));
}
t.stop();
CGAL_CLASSIFICATION_CERR << "Scales computed in " << t.time() << " second(s)" << std::endl;
t.reset();
}
@@ -359,6 +356,77 @@ public:
/// \endcond
void generate_point_based_features ()
{
#ifdef DO_NOT_USE_EIGEN_FEATURES
for (int j = 0; j < 3; ++ j)
{
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Eigenvalue> (i, m_input, eigen(i), std::size_t(j));
}
#else
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Anisotropy> (i, m_input, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Eigentropy> (i, m_input, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Linearity> (i, m_input, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Omnivariance> (i, m_input, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Planarity> (i, m_input, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Sphericity> (i, m_input, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Sum_eigen> (i, m_input, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Surface_variation> (i, m_input, eigen(i));
#endif
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Distance_to_plane> (i, m_input, m_point_map, eigen(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Dispersion> (i, m_input, m_point_map, grid(i), radius_neighbors(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Elevation> (i, m_input, m_point_map, grid(i), radius_dtm(i));
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Verticality> (i, m_input, eigen(i));
}
template <typename VectorMap>
void generate_normal_based_features(const VectorMap& normal_map)
{
m_features.add<Verticality> (m_input, normal_map);
}
template <typename ColorMap>
void generate_color_based_features(const ColorMap& color_map)
{
#ifdef DO_NOT_USE_HSV_FEATURES
typedef Feature::Color_channel<GeomTraits, PointRange, ColorMap> Color_channel;
for (std::size_t i = 0; i < 3; ++ i)
m_features.add<Color_channel> (m_input, color_map, typename Color_channel::Channel(i));
#else
typedef Feature::Hsv<GeomTraits, PointRange, ColorMap> Hsv;
for (std::size_t i = 0; i <= 8; ++ i)
m_features.add<Hsv> (m_input, color_map, typename Hsv::Channel(0), 45.f * float(i), 22.5f);
for (std::size_t i = 0; i <= 4; ++ i)
m_features.add<Hsv> (m_input, color_map, typename Hsv::Channel(1), 25.f * float(i), 12.5f);
for (std::size_t i = 0; i <= 4; ++ i)
m_features.add<Hsv> (m_input, color_map, typename Hsv::Channel(2), 25.f * float(i), 12.5f);
#endif
}
template <typename EchoMap>
void generate_echo_based_features(const EchoMap& echo_map)
{
typedef Feature::Echo_scatter<GeomTraits, PointRange, PointMap, EchoMap> Echo_scatter;
for (std::size_t i = 0; i < m_scales.size(); ++ i)
m_features.add_with_scale_id<Echo_scatter> (i, m_input, echo_map, grid(i), radius_neighbors(i));
}
/*!
\brief Returns the bounding box of the input point set.
*/
@@ -412,94 +480,6 @@ private:
m_scales.clear();
}
void add (Feature_handle fh, std::size_t i)
{
m_features_to_rename.push_back (std::make_pair (fh, i));
}
void generate_point_based_features ()
{
#ifdef DO_NOT_USE_EIGEN_FEATURES
for (int j = 0; j < 3; ++ j)
{
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Eigenvalue> (m_input, eigen(i), std::size_t(j)), i);
}
#else
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Anisotropy> (m_input, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Eigentropy> (m_input, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Linearity> (m_input, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Omnivariance> (m_input, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Planarity> (m_input, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Sphericity> (m_input, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Sum_eigen> (m_input, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Surface_variation> (m_input, eigen(i)), i);
#endif
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Distance_to_plane> (m_input, m_point_map, eigen(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Dispersion> (m_input, m_point_map, grid(i), radius_neighbors(i)), i);
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Elevation> (m_input, m_point_map, grid(i), radius_dtm(i)), i);
}
template <typename VectorMap>
void generate_normal_based_features(const VectorMap& normal_map)
{
m_features.add<Verticality> (m_input, normal_map);
}
void generate_normal_based_features(const CGAL::Constant_property_map<Iterator, typename GeomTraits::Vector_3>&)
{
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Verticality> (m_input, eigen(i)), i);
}
template <typename ColorMap>
void generate_color_based_features(const ColorMap& color_map)
{
#ifdef DO_NOT_USE_HSV_FEATURES
typedef Feature::Color_channel<GeomTraits, PointRange, ColorMap> Color_channel;
for (std::size_t i = 0; i < 3; ++ i)
m_features.add<Color_channel> (m_input, color_map, typename Color_channel::Channel(i));
#else
typedef Feature::Hsv<GeomTraits, PointRange, ColorMap> Hsv;
for (std::size_t i = 0; i <= 8; ++ i)
m_features.add<Hsv> (m_input, color_map, typename Hsv::Channel(0), 45.f * float(i), 22.5f);
for (std::size_t i = 0; i <= 4; ++ i)
m_features.add<Hsv> (m_input, color_map, typename Hsv::Channel(1), 25.f * float(i), 12.5f);
for (std::size_t i = 0; i <= 4; ++ i)
m_features.add<Hsv> (m_input, color_map, typename Hsv::Channel(2), 25.f * float(i), 12.5f);
#endif
}
void generate_color_based_features(const CGAL::Constant_property_map<Iterator, RGB_Color>&)
{
}
template <typename EchoMap>
void generate_echo_based_features(const EchoMap& echo_map)
{
typedef Feature::Echo_scatter<GeomTraits, PointRange, PointMap, EchoMap> Echo_scatter;
for (std::size_t i = 0; i < m_scales.size(); ++ i)
add(m_features.add<Echo_scatter> (m_input, echo_map, grid(i), radius_neighbors(i)), i);
}
void generate_echo_based_features(const CGAL::Constant_property_map<Iterator, std::size_t>&)
{
}
void generate_gradient_features()
{
#ifdef CGAL_CLASSIFICATION_USE_GRADIENT_OF_FEATURE
@@ -535,61 +515,6 @@ private:
return Constant_property_map<Iterator, T>();
}
template<typename VectorMap, typename ColorMap, typename EchoMap>
void generate_features_impl (std::size_t nb_scales, float voxel_size,
VectorMap normal_map,
ColorMap color_map,
EchoMap echo_map)
{
CGAL::Real_timer t; t.start();
m_scales.reserve (nb_scales);
m_scales.push_back (new Scale (m_input, m_point_map, m_bbox, -1.f));
if (voxel_size == -1.f)
voxel_size = m_scales[0]->grid_resolution();
for (std::size_t i = 1; i < nb_scales; ++ i)
{
voxel_size *= 2;
m_scales.push_back (new Scale (m_input, m_point_map, m_bbox, voxel_size, m_scales[i-1]->grid));
}
t.stop();
CGAL_CLASSIFICATION_CERR << "Scales computed in " << t.time() << " second(s)" << std::endl;
t.reset();
t.start();
#ifdef CGAL_LINKED_WITH_TBB
if (boost::is_convertible<ConcurrencyTag,Parallel_tag>::value)
m_features.begin_parallel_additions();
#else
CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
"Parallel_tag is enabled but TBB is unavailable.");
#endif
generate_point_based_features ();
generate_normal_based_features (normal_map);
generate_color_based_features (color_map);
generate_echo_based_features (echo_map);
#ifdef CGAL_LINKED_WITH_TBB
if (boost::is_convertible<ConcurrencyTag,Parallel_tag>::value)
m_features.end_parallel_additions();
#else
CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
"Parallel_tag is enabled but TBB is unavailable.");
#endif
for (std::size_t i = 0; i < m_features_to_rename.size(); ++ i)
m_features_to_rename[i].first->set_name
(m_features_to_rename[i].first->name() + "_" + std::to_string(m_features_to_rename[i].second));
t.stop();
CGAL_CLASSIFICATION_CERR << "Features computed in " << t.time() << " second(s)" << std::endl;
}
};