Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Avoid phantom surfces in marching cubes hoppe #1874

Merged
merged 12 commits into from
Aug 27, 2017
Merged

Avoid phantom surfces in marching cubes hoppe #1874

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
22a4977
spatial restriction for marching cube hoppe to avoid phantom surfaces…
Mar 10, 2017
778adec
requested changes for pr
Jul 3, 2017
97969d6
space, reserve of intermediate cloud
Jul 3, 2017
10fe1eb
should have made min_p_ and max_p_ local
Jul 3, 2017
0334173
remove min_p_ and max_p_ from mc rbf. remove debug info for mc get bo…
Jul 3, 2017
e918642
move constructor to header. set boundary and voxel size to class memb…
Jul 3, 2017
6cc5a6d
remove unnecessary local bounding box
Jul 3, 2017
13b2a8f
clean undeleted comment and cut long statement
Jul 3, 2017
a5cea53
some minor changes
Jul 4, 2017
f08107a
fake commit, to make sure bb is the problem
Jul 4, 2017
96d1eb4
undo last commit and fix gtest data for marching cube. if travis test…
Jul 4, 2017
3a7390f
make getter and setter for dist_ignore inline
Jul 5, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 35 additions & 24 deletions surface/include/pcl/surface/impl/marching_cubes.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@
//////////////////////////////////////////////////////////////////////////////////////////////
template <typename PointNT>
pcl::MarchingCubes<PointNT>::MarchingCubes ()
: min_p_ (), max_p_ (), percentage_extend_grid_ (), iso_level_ (), dist_ignore_ (-1.0f)
: percentage_extend_grid_ (), iso_level_ (), dist_ignore_ (-1.0f)
{
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wasn't this supposed to be moved to the .h file?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Move this entire constructor to the .h


Expand All @@ -57,16 +57,22 @@ pcl::MarchingCubes<PointNT>::~MarchingCubes ()

//////////////////////////////////////////////////////////////////////////////////////////////
template <typename PointNT> void
pcl::MarchingCubes<PointNT>::getBoundingBox ()
pcl::MarchingCubes<PointNT>::getBoundingBox (Eigen::Vector3f &upper_boundary,
Eigen::Vector3f &lower_boundary) const
{
pcl::getMinMax3D (*input_, min_p_, max_p_);
PointNT max_pt, min_pt;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're just interested in the coordinates so this should be PointXYZ.

pcl::getMinMax3D (*input_, min_pt, max_pt);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It doesn't make much difference in this case but this looks like the overload you were actually looking for
http://docs.pointclouds.org/trunk/group__common.html#gafd9010977f5e52b35b484be7624df3f8


min_p_ -= (max_p_ - min_p_) * percentage_extend_grid_/2;
max_p_ += (max_p_ - min_p_) * percentage_extend_grid_/2;
lower_boundary = min_pt.getVector3fMap ();
upper_boundary = max_pt.getVector3fMap ();

Eigen::Vector4f bounding_box_size = max_p_ - min_p_;
Eigen::Vector3f size3_extend = 0.5f * percentage_extend_grid_ * (upper_boundary - lower_boundary);

lower_boundary -= size3_extend;
upper_boundary += size3_extend;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Everything below this line has no effect, right? Should we delete these lines?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have deleted. Is the constructor now what you want?


Eigen::Vector3f bounding_box_size = upper_boundary - lower_boundary;

bounding_box_size = max_p_ - min_p_;
PCL_DEBUG ("[pcl::MarchingCubesHoppe::getBoundingBox] Size of Bounding Box is [%f, %f, %f]\n",
bounding_box_size.x (), bounding_box_size.y (), bounding_box_size.z ());
double max_size =
Expand All @@ -76,9 +82,9 @@ pcl::MarchingCubes<PointNT>::getBoundingBox ()
// ????
// data_size_ = static_cast<uint64_t> (max_size / leaf_size_);
PCL_DEBUG ("[pcl::MarchingCubesHoppe::getBoundingBox] Lower left point is [%f, %f, %f]\n",
min_p_.x (), min_p_.y (), min_p_.z ());
lower_boundary.x (), lower_boundary.y (), lower_boundary.z ());
PCL_DEBUG ("[pcl::MarchingCubesHoppe::getBoundingBox] Upper left point is [%f, %f, %f]\n",
max_p_.x (), max_p_.y (), max_p_.z ());
upper_boundary.x (), upper_boundary.y (), upper_boundary.z ());
}


Expand All @@ -99,6 +105,8 @@ pcl::MarchingCubes<PointNT>::interpolateEdge (Eigen::Vector3f &p1,
template <typename PointNT> void
pcl::MarchingCubes<PointNT>::createSurface (const std::vector<float> &leaf_node,
const Eigen::Vector3i &index_3d,
const Eigen::Vector3f &upper_boundary,
const Eigen::Vector3f &lower_boundary,
pcl::PointCloud<PointNT> &cloud)
{
int cubeindex = 0;
Expand All @@ -116,25 +124,24 @@ pcl::MarchingCubes<PointNT>::createSurface (const std::vector<float> &leaf_node,
if (edgeTable[cubeindex] == 0)
return;

//Eigen::Vector4f index_3df (index_3d[0], index_3d[1], index_3d[2], 0.0f);
Eigen::Vector3f center;// TODO coeff wise product = min_p_ + Eigen::Vector4f (1.0f/res_x_, 1.0f/res_y_, 1.0f/res_z_) * index_3df * (max_p_ - min_p_);
center[0] = min_p_[0] + (max_p_[0] - min_p_[0]) * float (index_3d[0]) / float (res_x_);
center[1] = min_p_[1] + (max_p_[1] - min_p_[1]) * float (index_3d[1]) / float (res_y_);
center[2] = min_p_[2] + (max_p_[2] - min_p_[2]) * float (index_3d[2]) / float (res_z_);
const Eigen::Vector3f size_voxel = (upper_boundary - lower_boundary).cwiseQuotient (
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is similar to delta right? Increases the time you compute this to three time.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

excuse me, what do you mean?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

https://github.com/t-lou/pcl/blob/10fe1eb3e7e720a4a80f8ae3c304a493eeb77785/surface/include/pcl/surface/impl/marching_cubes.hpp#L127-L128

is similar to this

https://github.com/t-lou/pcl/blob/10fe1eb3e7e720a4a80f8ae3c304a493eeb77785/surface/include/pcl/surface/impl/marching_cubes_hoppe.hpp#L66-L67

and this

https://github.com/t-lou/pcl/blob/10fe1eb3e7e720a4a80f8ae3c304a493eeb77785/surface/include/pcl/surface/impl/marching_cubes_rbf.hpp#L103-L104

You're computing the same property three times. This normally hints me that it should become a member of the base class. What's your opinion?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would prefer to calculate once and use them as parameters, as they are data-related, not algorithm related. In this way the functions would be more stateless.
Do you think I should put delta (voxel_size), min_p_ and max_p_ all to class definition?

Copy link
Member

@SergioRAgostinho SergioRAgostinho Jul 3, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. I'm just trying to understand what are the advantages of your approach and the avantages of mine.

The class has data members (i.e. input_ and grid_) so it's only natural to have other members related to this data.

Currently your protected methods have a parameters in the signature that never changes, as long as the input cloud (which is a member) doesn't change, and they're only invoked following a public call to performReconstruction.

If we add the members to the class definition and simply don't expose setters and getter functions to them we have:

  • More member which are data related but not accessible to the outside
  • You don't have them in the signature because you don't need them, since you're anyway you're alway invoking the methods with the same bound data.
  • You prevent have to recompute delta, because it persists valid throughout a call to performReconstruction

Opinion @taketwo

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

PCL classes are by design stateful. We use separate function calls to set parameters, input point cloud, input indices. After that we call a compute method, passing a single argument: storage for the output. Therefore I agree with @SergioRAgostinho, there is no reason to pretend there is no state and pass around variables that are computed from the input point cloud, which in fact is the state.

Eigen::Vector3f (res_x_, res_y_, res_z_));

const Eigen::Vector3f center = lower_boundary + size_voxel.cwiseProduct (index_3d.cast<float> ());

std::vector<Eigen::Vector3f, Eigen::aligned_allocator<Eigen::Vector3f> > p;
p.resize (8);
for (int i = 0; i < 8; ++i)
{
Eigen::Vector3f point = center;
if(i & 0x4)
point[1] = static_cast<float> (center[1] + (max_p_[1] - min_p_[1]) / float (res_y_));
point[1] = static_cast<float> (center[1] + size_voxel[1]);

if(i & 0x2)
point[2] = static_cast<float> (center[2] + (max_p_[2] - min_p_[2]) / float (res_z_));
point[2] = static_cast<float> (center[2] + size_voxel[2]);

if((i & 0x1) ^ ((i >> 1) & 0x1))
point[0] = static_cast<float> (center[0] + (max_p_[0] - min_p_[0]) / float (res_x_));
point[0] = static_cast<float> (center[0] + size_voxel[0]);

p[i] = point;
}
Expand Down Expand Up @@ -237,7 +244,7 @@ pcl::MarchingCubes<PointNT>::performReconstruction (pcl::PolygonMesh &output)

performReconstruction (points, output.polygons);

pcl::toPCLPointCloud2(points, output.cloud);
pcl::toPCLPointCloud2 (points, output.cloud);
}


Expand All @@ -259,24 +266,28 @@ pcl::MarchingCubes<PointNT>::performReconstruction (pcl::PointCloud<PointNT> &po
// the point cloud really generated from Marching Cubes, prev intermediate_cloud_
pcl::PointCloud<PointNT> intermediate_cloud;

Eigen::Vector3f upper_boundary, lower_boundary;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I noticed that every operation on these guys is coefficient wise. You should replace their type to Eigen::Array. Have a look at https://eigen.tuxfamily.org/dox/group__TutorialArrayClass.html


// Create grid
grid_ = std::vector<float> (res_x_*res_y_*res_z_, NAN);

// Populate tree
tree_->setInputCloud (input_);

getBoundingBox ();
getBoundingBox (upper_boundary, lower_boundary);

// Transform the point cloud into a voxel grid
// This needs to be implemented in a child class
voxelizeData ();
voxelizeData (upper_boundary, lower_boundary);

// Run the actual marching cubes algorithm, store it into a point cloud,
// and copy the point cloud + connectivity into output
intermediate_cloud.clear();
intermediate_cloud.clear ();
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unnecessary. Intermediate_cloud is always initialized empty.


// preallocate memory assuming a hull. suppose 6 point per voxel
intermediate_cloud.reserve((size_t)(res_y_*res_z_ + res_x_*res_z_ + res_x_*res_y_) * 2 * 6);
double size_reserve = std::min((double) intermediate_cloud.points.max_size (),
2.0 * 6.0 * (double) (res_y_*res_z_ + res_x_*res_z_ + res_x_*res_y_));
intermediate_cloud.reserve ((size_t) size_reserve);

for (int x = 1; x < res_x_-1; ++x)
for (int y = 1; y < res_y_-1; ++y)
Expand All @@ -286,10 +297,10 @@ pcl::MarchingCubes<PointNT>::performReconstruction (pcl::PointCloud<PointNT> &po
std::vector<float> leaf_node;
getNeighborList1D (leaf_node, index_3d);
if(!leaf_node.empty())
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Space after if and empty

createSurface (leaf_node, index_3d, intermediate_cloud);
createSurface (leaf_node, index_3d, upper_boundary, lower_boundary, intermediate_cloud);
}

points.swap(intermediate_cloud);
points.swap (intermediate_cloud);

polygons.resize (points.size () / 3);
for (size_t i = 0; i < polygons.size (); ++i)
Expand Down
16 changes: 8 additions & 8 deletions surface/include/pcl/surface/impl/marching_cubes_hoppe.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,13 +58,13 @@ pcl::MarchingCubesHoppe<PointNT>::~MarchingCubesHoppe ()

//////////////////////////////////////////////////////////////////////////////////////////////
template <typename PointNT> void
pcl::MarchingCubesHoppe<PointNT>::voxelizeData ()
pcl::MarchingCubesHoppe<PointNT>::voxelizeData (const Eigen::Vector3f &upper_boundary,
const Eigen::Vector3f &lower_boundary)
{
const bool is_far_ignored = dist_ignore_ > 0.0f;

const Eigen::Vector3f min_p = min_p_.segment(0, 3);
const Eigen::Vector3f delta = ( max_p_ - min_p_ ).segment(0, 3).cwiseQuotient(
Eigen::Vector3f((float)res_x_, (float)res_y_, (float)res_z_) );
const Eigen::Vector3f delta = (upper_boundary - lower_boundary).cwiseQuotient (
Eigen::Vector3f (res_x_, res_y_, res_z_));

for (int x = 0; x < res_x_; ++x)
{
Expand All @@ -76,9 +76,9 @@ pcl::MarchingCubesHoppe<PointNT>::voxelizeData ()

for (int z = 0; z < res_z_; ++z)
{
std::vector<int> nn_indices(1, 0);
std::vector<float> nn_sqr_dists(1, 0.0f);
const Eigen::Vector3f point = min_p + delta.cwiseProduct(Eigen::Vector3f(x, y, z));
std::vector<int> nn_indices (1, 0);
std::vector<float> nn_sqr_dists (1, 0.0f);
const Eigen::Vector3f point = lower_boundary + delta.cwiseProduct (Eigen::Vector3f (x, y, z));
PointNT p;

p.getVector3fMap () = point;
Expand All @@ -89,7 +89,7 @@ pcl::MarchingCubesHoppe<PointNT>::voxelizeData ()
{
const Eigen::Vector3f normal = input_->points[nn_indices[0]].getNormalVector3fMap ();

if(normal.norm() > 0.5f)
if(!std::isnan (normal(0)) && normal.norm () > 0.5f)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Space after if and normal

grid_[z_start + z] = normal.dot (
point - input_->points[nn_indices[0]].getVector3fMap ());
}
Expand Down
12 changes: 7 additions & 5 deletions surface/include/pcl/surface/impl/marching_cubes_rbf.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,8 @@ pcl::MarchingCubesRBF<PointNT>::~MarchingCubesRBF ()

//////////////////////////////////////////////////////////////////////////////////////////////
template <typename PointNT> void
pcl::MarchingCubesRBF<PointNT>::voxelizeData ()
pcl::MarchingCubesRBF<PointNT>::voxelizeData (const Eigen::Vector3f &upper_boundary,
const Eigen::Vector3f &lower_boundary)
{
// Initialize data structures
unsigned int N = static_cast<unsigned int> (input_->size ());
Expand Down Expand Up @@ -99,14 +100,15 @@ pcl::MarchingCubesRBF<PointNT>::voxelizeData ()
weights[i + N] = w (i + N, 0);
}

const Eigen::Vector3d delta = (upper_boundary - lower_boundary).cast<double> ().cwiseQuotient (
Eigen::Vector3d (res_x_, res_y_, res_z_));
const Eigen::Vector3d base = lower_boundary.cast<double> ();

for (int x = 0; x < res_x_; ++x)
for (int y = 0; y < res_y_; ++y)
for (int z = 0; z < res_z_; ++z)
{
Eigen::Vector3d point;
point[0] = min_p_[0] + (max_p_[0] - min_p_[0]) * float (x) / float (res_x_);
point[1] = min_p_[1] + (max_p_[1] - min_p_[1]) * float (y) / float (res_y_);
point[2] = min_p_[2] + (max_p_[2] - min_p_[2]) * float (z) / float (res_z_);
Eigen::Vector3d point = base + delta.cwiseProduct (Eigen::Vector3d (x, y, z));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Const qualifier


double f = 0.0;
std::vector<double>::const_iterator w_it (weights.begin());
Expand Down
25 changes: 17 additions & 8 deletions surface/include/pcl/surface/marching_cubes.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -446,9 +446,6 @@ namespace pcl
/** \brief The grid resolution */
int res_x_, res_y_, res_z_;

/** \brief Min and max data points. */
Eigen::Vector4f min_p_, max_p_;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Was this (or the bound named equivalent) worth removing? I saw you compute it only once and then have to pass it around to a number of methods. Its lifetime is equivalent to input_, grid_, res_x_, etc... I really don't understand why you felt they were "stateful", or to be more clear, that they shouldn't belong to the state of the Class. Can you explain a little bit the reasoning?

I also question if delta should not follow the same.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I just have another compelling argument: There's not a single time you invoke the methods (the ones which now have the bounds on the signature) and that these bounds actually change. You don't invoke the method changing the bounds ever. As the long as the input cloud is the same, the bounds will always be the same.

Copy link
Contributor Author

@t-lou t-lou Jul 3, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because I think it's not the property of the class, but of the data. And using them as parameter doesn't make the program too complicated or slow.

Sorry, what do you mean by "if delta should not follow the same"?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you mean to do marching cube meshing on the same input cloud for several times with different parameters? But the old code would also invoke "getBoundingBox ()" each time we execute "performReconstruction"


/** \brief Parameter that defines how much free space should be left inside the grid between
* the bounding box of the point cloud and the grid limits, as a percentage of the bounding box.*/
float percentage_extend_grid_;
Expand All @@ -461,9 +458,13 @@ namespace pcl
* or distance between voxel centroid and point are larger that it. */
float dist_ignore_;

/** \brief Convert the point cloud into voxel data. */
/** \brief Convert the point cloud into voxel data.
* \param[in] upper_boundary The upper boundary of point cloud (after extension)
* \param[in] lower_boundary The upper boundary of point cloud (after extension)
*/
virtual void
voxelizeData () = 0;
voxelizeData (const Eigen::Vector3f &upper_boundary,
const Eigen::Vector3f &lower_boundary) = 0;

/** \brief Interpolate along the voxel edge.
* \param[in] p1 The first point on the edge
Expand All @@ -479,16 +480,24 @@ namespace pcl
/** \brief Calculate out the corresponding polygons in the leaf node
* \param leaf_node the leaf node to be checked
* \param index_3d the 3d index of the leaf node to be checked
* \param upper_boundary the upper boundary of point cloud
* \param lower_boundary the lower boundary of point cloud
* \param cloud point cloud to store the vertices of the polygon
*/
*/
void
createSurface (const std::vector<float> &leaf_node,
const Eigen::Vector3i &index_3d,
const Eigen::Vector3f &upper_boundary,
const Eigen::Vector3f &lower_boundary,
pcl::PointCloud<PointNT> &cloud);

/** \brief Get the bounding box for the input data points. */
/** \brief Get the bounding box for the input data points.
* \param[in] upper_boundary The upper boundary of point cloud (after extension)
* \param[in] lower_boundary The upper boundary of point cloud (after extension)
*/
void
getBoundingBox ();
getBoundingBox (Eigen::Vector3f &upper_boundary,
Eigen::Vector3f &lower_boundary) const;


/** \brief Method that returns the scalar value at the given grid position.
Expand Down
11 changes: 6 additions & 5 deletions surface/include/pcl/surface/marching_cubes_hoppe.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,8 +61,6 @@ namespace pcl
using MarchingCubes<PointNT>::res_x_;
using MarchingCubes<PointNT>::res_y_;
using MarchingCubes<PointNT>::res_z_;
using MarchingCubes<PointNT>::min_p_;
using MarchingCubes<PointNT>::max_p_;
using MarchingCubes<PointNT>::dist_ignore_;

typedef typename pcl::PointCloud<PointNT>::Ptr PointCloudPtr;
Expand All @@ -77,10 +75,13 @@ namespace pcl
/** \brief Destructor. */
~MarchingCubesHoppe ();

/** \brief Convert the point cloud into voxel data. */
/** \brief Convert the point cloud into voxel data.
* \param[in] upper_boundary The upper boundary of point cloud (after extension)
* \param[in] lower_boundary The lower boundary of point cloud (after extension)
*/
void
voxelizeData ();

voxelizeData (const Eigen::Vector3f &upper_boundary,
const Eigen::Vector3f &lower_boundary);

public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
Expand Down
8 changes: 6 additions & 2 deletions surface/include/pcl/surface/marching_cubes_rbf.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,9 +78,13 @@ namespace pcl
/** \brief Destructor. */
~MarchingCubesRBF ();

/** \brief Convert the point cloud into voxel data. */
/** \brief Convert the point cloud into voxel data.
* \param[in] upper_boundary The upper boundary of point cloud (after extension)
* \param[in] lower_boundary The lower boundary of point cloud (after extension)
*/
void
voxelizeData ();
voxelizeData (const Eigen::Vector3f &upper_boundary,
const Eigen::Vector3f &lower_boundary);


/** \brief Set the off-surface points displacement value.
Expand Down