[surface] getArea(point cloud segment)

[mynickmynick]

trafficstars

Is your feature request related to a problem? Please describe.

I am developing segmentation algorithms in which I need to calculate the area of segments to make selections and I think the implementation of area calculation i did works well and can be applied general purpouse, so it could be merged in PCL

Solution

Essentially there are a few steps: (1) optional voxeling of the segment (2) NormalEstimationOMP (3) GreedyProjectionTriangulation (4) Sum the area of greedy triangles Loop: area += 0.5 * P.cross(Q).norm(); where P and Q are two vectorial sides of the triangles

Describe alternatives you've considered

one might consider alternative surface mesh/triangulations but this seems to work fine on the samples I tested it there is a little are overestimate due to the fact that the mesh is not perfectly smooth. the voxeling is aimed at reducing this effect.

the area += loop is done as double

An alternative would be to extract optionally this info directly inside the greedy alg. as it comes available so to minimnize the total processing time

[mynickmynick]

i have pushed on my fork branch fork_getArea the modification which could be merged

there are both methods of GreedyProjectionTriangulation to set a flag for optional area calculation and retrieve the results (the area is updated inside addTriangle())

  /** \brief Get the calcArea flag value.
   ** true to calculate the area of the mesh*/
  inline bool
  getCalcArea() const { return (calcArea); }

  /** \brief Set the calcArea flag value.
   ** true to calculate the area of the mesh*/
  inline void
  setCalcArea(bool value) { calcArea = value; }

  /** \brief Get the Area value.*/
  inline double
  getArea() const { return (area); }

(the flag calcArea is set to false by default by the constructor) or a function that compose before optional voxeling, normal calculations etc. (one can use one or the other, but the function does not support indeces yet, while the GreedyProjectionTriangulation does)

template <typename PointInT, typename PointN> double greedyArea( shared_ptr<pcl::PointCloud<PointInT>>& cloud, double greedySearchRadius, // suggested >= 2.0*max(EuclideanConnection,VoxelLeaf) int normal_k_search = 20, float VoxelLeaf = 0.001, bool enableVoxel = false, double greedyMu = 2.5, int greedyMaxNearestNeighbours = 100, double greedyMaximumSurfaceAngle = (M_PI / 4), // 45 degrees double greedyMinimumAngle = (M_PI / 18), // 10 degrees double greedyMaximumAngle = (2 * M_PI / 3), // 120 degrees bool greedyNormalConsistency = false, unsigned int threadNrNormals = 2);

the area is updated inside addTriangle():

    if (calcArea) {
      const Eigen::Matrix<double, 3, 1> P(
          (*input_)[(*indices_)[triangle_.vertices[0]]].x -
              (*input_)[(*indices_)[triangle_.vertices[2]]].x,
          (*input_)[(*indices_)[triangle_.vertices[0]]].y -
              (*input_)[(*indices_)[triangle_.vertices[2]]].y,
          (*input_)[(*indices_)[triangle_.vertices[0]]].z -
              (*input_)[(*indices_)[triangle_.vertices[2]]].z);
      const Eigen::Matrix<double, 3, 1> Q(
          (*input_)[(*indices_)[triangle_.vertices[1]]].x -
              (*input_)[(*indices_)[triangle_.vertices[2]]].x,
          (*input_)[(*indices_)[triangle_.vertices[1]]].y -
              (*input_)[(*indices_)[triangle_.vertices[2]]].y,
          (*input_)[(*indices_)[triangle_.vertices[1]]].z -
              (*input_)[(*indices_)[triangle_.vertices[2]]].z);
      area += 0.5 * P.cross(Q).norm();
    }

[mynickmynick]

https://github.com/mynickmynick/pcl/tree/fork_getArea

[mynickmynick]

An alternative would be to extract optionally this info directly inside the greedy alg. as it comes available so to minimnize the total processing time

I am quoting myself to tell that i did that (inside addTriangle )

but I can extract that part of code again outside of GreedyProjectionTriangulation running on the triangle mesh it outputs. I checked the elaboration time is nearly the same

[mynickmynick]

you might wanna merge only the methods for area calculation added to GreedyProjectionTriangulation and leave the function greedyArea() only as a tutorial

[mvieth]

There is already a tutorial that shows how to estimate normals and then use GreedyProjectionTriangulation: https://pcl.readthedocs.io/projects/tutorials/en/master/greedy_projection.html

I checked and did not find a method that computes the total area of a mesh. Closest I found is calculatePolygonArea.

Since computing the mesh area is not specific to GreedyProjectionTriangulation (can be done on any (triangle) mesh no matter where it is from), I would not add extra functionality to that class but instead add a free function that does this, maybe in pcl/common/common.h

[mynickmynick]

There is already a tutorial that shows how to estimate normals and then use GreedyProjectionTriangulation: https://pcl.readthedocs.io/projects/tutorials/en/master/greedy_projection.html

yes I had seen it and refered to it as a matter of facts, but it doesn't compute the area

Since computing the mesh area is not specific to GreedyProjectionTriangulation (can be done on any (triangle) mesh no matter where it is from), I would not add extra functionality to that class but instead add a free function that does this, maybe in pcl/common/common.h

I don't see clearly how to extend it to any mesh calculated, how to test if the mesh is of triangles and if the mesh is not redundant (with partially overlapping faces or similar). please name some other classes that would be involved or give frtehr hints

i see other 3 classes with a common base

class ConvexHull : public MeshConstruction<PointInT> it already has the area calculation with different methods coming from qhull and it tends to give a 2 area estimation when the cloud is flat because it meshes it with a convex double layer (basically wrapping around the cloud with a double surace). As a matter of facts I am using it to estimate flat convexity of the cloud by comparing 2greedyarea and convex hull area!!

class ConcaveHull : public MeshConstruction<PointInT> this has not implemented it in PCL but has it available in qhull an area and volume computation analog to convexhull that can be estracted (i had done some tests but using Delaunay it tends to give one dimension more calculation, not very clear...)

class OrganizedFastMesh : public MeshConstruction<PointInT>: I have no idea how this works

I suggest we give this implementation available for Greedy alg. which is simple and already available as an optional calculation (as already calc area and volume is optionally available for convex hull) and postpone extensions to other scenarios for later

[mynickmynick]

or I can just make available a function that does this

        double area = 0;
        for ()
        {
            const Eigen::Matrix<double, 3, 1>
              P((*input_)[(*indices_)[triangle_.vertices[0]]].x - (*input_)[(*indices_)[triangle_.vertices[2]]].x,
                (*input_)[(*indices_)[triangle_.vertices[0]]].y - (*input_)[(*indices_)[triangle_.vertices[2]]].y,
                (*input_)[(*indices_)[triangle_.vertices[0]]].z - (*input_)[(*indices_)[triangle_.vertices[2]]].z);
            const Eigen::Matrix<double, 3, 1>
              Q((*input_)[(*indices_)[triangle_.vertices[1]]].x - (*input_)[(*indices_)[triangle_.vertices[2]]].x,
                (*input_)[(*indices_)[triangle_.vertices[1]]].y - (*input_)[(*indices_)[triangle_.vertices[2]]].y,
                (*input_)[(*indices_)[triangle_.vertices[1]]].z - (*input_)[(*indices_)[triangle_.vertices[2]]].z);
            area += 0.5 * P.cross(Q).norm();
        }

or

        for (const auto& vv : triangles.polygons)
        {
            const Eigen::Matrix<double, 3, 1>
              P((*input_)[(*indices_)[vv.vertices[0]]].x - (*input_)[(*indices_)[v.vertices[2]]].x,
                (*input_)[(*indices_)[v.vertices[0]]].y - (*input_)[(*indices_)[v.vertices[2]]].y,
                (*input_)[(*indices_)[v.vertices[0]]].z - (*input_)[(*indices_)[v.vertices[2]]].z);
            const Eigen::Matrix<double, 3, 1>
              Q((*input_)[(*indices_)[v.vertices[1]]].x - (*input_)[(*indices_)[v.vertices[2]]].x,
                (*input_)[(*indices_)[v.vertices[1]]].y - (*input_)[(*indices_)[v.vertices[2]]].y,
                (*input_)[(*indices_)[v.vertices[1]]].z - (*input_)[(*indices_)[v.vertices[2]]].z);
            area += 0.5 * P.cross(Q).norm();

        }

[mynickmynick]

but there triangles would be pcl::PolygonMesh and I would have to waste time to check if the polygons are triangles or have more than 3 points it's inefficient

[mvieth]

I don't see clearly how to extend it to any mesh calculated, how to test if the mesh is of triangles and if the mesh is not redundant (with partially overlapping faces or similar)

and I would have to waste time to check if the polygons are triangles or have more than 3 points it's inefficient

It is not necessary that the function checks that every face is a triangle. It is sufficient if the function name and/or function documentation makes it clear that this is a requirement. Then the user is responsible.

please name some other classes that would be involved or give frtehr hints

Yes, OrganizedFastMesh could for example also produce a triangle mesh where a user wants to know the total area. Maybe the marching cubes algorithms as well, would have to check whether they guarantee triangle meshes.

as already calc area and volume is optionally available for convex hull

ConvexHull does not compute area and volume itself, it only makes available/acts as an interface to qhull. So this is not comparable in my opinion.

[mynickmynick]

please let me make some observations in points: (1) I had implemented the area calculation initially outside GreedyProjectionTriangulation methods, then I implemented it inside . I checked only with clouds not indexed and saw no difference in elaboration time. But I expect the second (area calculated inside addTriangle()) to be faster with cloud mapped by Indices. Because in the first case (outside in the function) you would have a second cycle with indirect mapping in RAM, while in thr first everything is done in just one cycle which the compiler can optimize for RAM access (and it can also be parallelized in the future). So I would recommend to have it optionally calculated inside. If your concern is to avoid to add the if (calcArea) {} inside addTriangle() we can also implement a second alternative like void addTriangleAndArea (pcl::index_t a, pcl::index_t b, pcl::index_t c, std::vectorpcl::Vertices &polygons) with the area calculation and leave the original addTriangle() without it

but as they are inline we would need two overloads or alternatives also for the calling functions so i think it's not worthwhile. it's just an if(bool) added condition with negligible elaboration time effect

(2) in industrial automation it is important to have simple reliable fast robust and efficient code which checks and handles all errors how can you do that by allowing the input to the function to be a mesh not of triangles and let the function compute anyway an area which is not correct as an output if the mesh is different? you cannot let the user do that leaving the checks to the comments/documentation that can work for an academic library, not for a library with also industrial applications

Like this it is already available to merge, very simple and efficient, few lines of code. One could test if expressing explicitely the norm(cross product) would speed up (and it probably would), leaving the //area += 0.5 * P.cross(Q).norm();
just as a comment

ConvexHull does not compute area and volume itself, it only makes available/acts as an interface to qhull. So this is not comparable in my opinion.

yes but it is just another example of area made available internally, optionally and in a specialized efficient way.

Yes, OrganizedFastMesh could for example also produce a triangle mesh where a user wants to know the total area. Maybe the marching cubes algorithms as well, would have to check whether they guarantee triangle meshes.

I don't think I have time in the following days to study that and it might take a long time. Why not postponing it to future implementations? Why pretending to generalize to different methods that will probably require a different and more complicated approach?

[mynickmynick]

and also the documentation is usually so incomplete that... I am now extending to the case of Indices and it is not clear at all how one is supposed to use the voxel on indexed clouds. you have to call getRemovedIndices after?? how am i supposed to collect the new indexing after voxeling? i find them automatically changed on the same Indices pointer?? I am using un-indexed copies of the segments that makes it simpler and with faster access.

[mynickmynick]

@mvieth what shall we do about this old but interesting topic i opened?

[mvieth]

Honestly, I still think a free function, for example called computeTriangleMeshArea, would be the most versatile solution, and applicable to not only GreedyProjectionTriangulation. Regarding your point that program correctness cannot/should not be ensured by the documentation: if a user wants to shoot himself/herself in the foot, they will probably find a way to do that with most existing PCL class. It is not possible, or in the interest of speed at least not desirable, to check for all user input whether all conventions/ preconditions stated in the documentation are fulfilled. Such checks are a speed/security trade-off and should be decided on a case-by-case basis. I think it is reasonable to expect a user to read the documentation and/or make basic inferences from the function name (also in industrial applications :wink:)

Perhaps we can ask @larshg for a third opinion.

In the end, of course, I can't force you to implement something you don't want. If you are only interested in contributing this addition to GreedyProjectionTriangulation, and we make sure that the class does not become slower if the area computation is turned off, then I am also fine with merging that :slightly_smiling_face:

[mynickmynick]

OK thanks! I will review this probably not before end of August or beginning of September because next weeks I'll have a long trip during which I won't be much available

[mynickmynick]

@mvieth i pushed the functions