mst_orient_normals Memory issue

[ahmadhasan2k8]

Issue Details

I have a multiple scanned point clouds where I calculate normals using CGAL. I run a Screened Poisson Reconstruction afterwards so I also orient the normals using mst_orient_normals. Now, because the data is scanned, sometimes there are islands of points even after some outlier removal which makes orienting the normals difficult unless I increase the nb_neighbors for mst_orient_normals to almost 500. Each point cloud have almost million points and I use MPI to run several process simultaneously. Problem is, I run out of memory beacuse I think CGAL stores all the neighborhood information in memory and process them altogether. So, even if I can use 18 or some small number for nb_neighbors and have a point cloud of 10~20 millions, and I run say 20 threads, it requires massive memory. Even my swap space becomes almost full, sometimes fails too. If I double the nb_neighbors the Mb allocated also almost doubles. Is there a way to use lower number of neighbors but still get correct orientation in case of islands of point cloud? Is it possible to not store all the data in the memory but process them as you get? I have attached a sample point cloud for testing too. nParticleShape1.2.xyz.zip

Source Code

main.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/read_xyz_points.h>
#include <CGAL/Memory_sizer.h>
#include <CGAL/Timer.h>
#include <CGAL/mst_orient_normals.h>
#include <CGAL/pca_estimate_normals.h>
#include <CGAL/property_map.h>

#include <fstream>
#include <list>
#include <utility>  // defines std::pair

// Types
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;

// Point with normal vector stored in a std::pair.
typedef std::pair<Point, Vector> PointVectorPair;

// Concurrency
#ifdef CGAL_LINKED_WITH_TBB
typedef CGAL::Parallel_tag Concurrency_tag;PLY
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif

int main(int argc, char* argv[])
{
  const char* fname = (argc > 1) ? argv[1] : "nParticleShape1.2.xyz";
  // Reads a .xyz point set file in points[].
  std::list<PointVectorPair> points;
  std::ifstream stream(fname);
  if (!stream ||
      !CGAL::read_xyz_points(
          stream,
          std::back_inserter(points),
          CGAL::parameters::point_map(
              CGAL::First_of_pair_property_map<PointVectorPair>())))
  {
    std::cerr << "Error: cannot read file " << fname << std::endl;
    return EXIT_FAILURE;
  }

  // Estimates normals direction.
  // Note: pca_estimate_normals() requiresa range of points
  // as well as property maps to access each point's position and normal.
  const int nb_neighbors = 18;  // K-nearest neighbors = 3 rings
  CGAL::pca_estimate_normals<Concurrency_tag>(
      points,
      nb_neighbors,
      CGAL::parameters::point_map(
          CGAL::First_of_pair_property_map<PointVectorPair>())
          .normal_map(CGAL::Second_of_pair_property_map<PointVectorPair>()));

  // Orients normals.
  // Note: mst_orient_normals() requires a range of points
  // as well as property maps to access each point's position and normal.
  std::list<PointVectorPair>::iterator unoriented_points_begin =
      CGAL::mst_orient_normals(
          points,
          1000,
          CGAL::parameters::point_map(
              CGAL::First_of_pair_property_map<PointVectorPair>())
              .normal_map(
                  CGAL::Second_of_pair_property_map<PointVectorPair>()));
  std::size_t memory = CGAL::Memory_sizer().virtual_size();
  std::cerr << "done: " << (memory >> 20) << " Mb allocated" << std::endl;

  // Optional: delete points with an unoriented normal
  // if you plan to call a reconstruction algorithm that expects oriented
  // normals.
  points.erase(unoriented_points_begin, points.end());

  return EXIT_SUCCESS;
}

CMakeLists.txt

cmake_minimum_required(VERSION 3.5)
project(CGAL_MST)

set(PROJECT_SRCS
${PROJECT_SOURCE_DIR}/src/main.cpp)

source_group("src" FILES ${PROJECT_SRCS})

find_package( CGAL REQUIRED )

add_executable(${PROJECT_NAME} ${PROJECT_SRCS})

target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_14)
target_link_libraries(${PROJECT_NAME} CGAL::CGAL )

Environment

• Operating system (Windows/Mac/Linux, 32/64 bits):ubunutu 16.04
• Compiler: gcc 5.4
• Release or debug mode:
• Specific flags used (if any):
• CGAL version: Latest
• Boost version: Boost 1.58
• Other libraries versions if used (Eigen, TBB, etc.):

[afabri]

The problem is that mst_orient_normals() creates a graph where each vertex has an edge to its closest 1000 vertices (that's the Riemannian Graph, the manual mentions). When you say that you have islands then this means that you need this big number so that the Riemannian graph becomes one connected component.
@palliez are you aware of any newer algorithms for orientinig normals, or variants of the Riemannian graph, say several passes with an increasing number of neighbors to consider?

[afabri]

And do you encounter a problem with islands with nParticleShape1.2.xyz.zip ??? Everything looks good for me.

[ahmadhasan2k8]

Not for nParticleShape1.2.xyz. I put this data to show memory usage. I get that forming a reimannnian graph might need more meory usage, I was more like suggesting if it was possible to form the graph iteratively somehow so that the memory usage goes down. The uploaded point cloud Des not have islands. If there are islands increasing number of neighbors create memory issues.

However, I got past the problem by using minimum number of neighbors though, 18 as suggested. Steps for that:

1. Orient normals.
2. Check if there is still unoriented normals.
3. If yes, run mst_orient_normal on the unoriented cloud.
4. Repeat till there is none left.
5. Append all the points in a single storage.

That seems to work perfectly and have no memory issues. This way I can solve the island problem and also don't have to use a large number of neighbors

[ahmadhasan2k8]

I could write an example code if needed and give a pull request if you want.

[afabri]

You still might get islands that are oriented differently from the "continent", that is the points on the continent may show upwards, and the points on the island downwards. The connectivity is needed for the global propagation. But maybe we can together come up with something.

[ahmadhasan2k8]

That's true. It does happen. If the island is too small I just remove it, but yes it's not a global propagation though. But the main purpose for me is to use it for poisson reconstruction. Even if the continent and islands might be oriented differently, poisson reconstruction works just fine. However it doesn't work well when individual connected components have different orientation. So, the process I mentioned above works really well without memory issue. But I do agree that if we could come up with something better that does not have memory issues and can propagate normals globally that would be great.

[afabri]

Are you willing to share a problematic data set?

[ahmadhasan2k8]

I probably can't share the original data but I can try to create a similar type of data which could be problematic. I'll give it a shot over the weekend, should not be too difficult

[ahmadhasan2k8]

CGAL.zip This could work as an ideal sample data.

[ahmadhasan2k8]

How about picking one point from all the islands each and one from the continent to form a graph to span a tree. And then use the continent point as seed for propagation. At this point all the islands will have one point which is globally oriented as the continent. Then, use those points as seed to propagate the orientation in their respective islands. I'm not very sure if this could work but seems to be the right way.