Why does an `itk::PointSet` have "PointData"? Any existing use cases?

Niels_Dekker · September 30, 2024, 1:20pm

Obviously, itk::PointSet has a container of points (m_PointsContainer), but it also may have some extra data, named “PointData”, and stored in its m_PointDataContainer.

github.com

InsightSoftwareConsortium/ITK/blob/3b04ad1b377c6d5d77a2c1f454c363eab8197340/Modules/Core/Common/include/itkPointSet.h#L144


      
          
          protected:
            /** An object containing points used by the mesh.  Individual points are
             * accessed through point identifiers. */
            PointsContainerPointer m_PointsContainer{};
          
            /** An object containing data associated with the mesh's points.
             * Optionally, this can be nullptr, indicating that no data are associated with
             * the points.  The data for a point can be accessed through its point
             * identifier. */
            PointDataContainerPointer m_PointDataContainer{};
          
          public:
            /** Copy the geometric and topological structure of the given input pointSet.
             * The copying is done via reference counting.
             */
            void
            PassStructure(Self * inputPointSet);
          
            /** Restore the PointSet to its initial state. Useful for data pipeline updates
             * without memory re-allocation.

As was introduced by @willschroeder, more than 23 years ago:

Can anyone please tell me some existing use cases of these “PointData”?

dzenanz · September 30, 2024, 1:45pm

It could be per-vertex color, or some other per-vertex measurement (thickness, curvature, strain, metabolic activity, etc).

Niels_Dekker · October 1, 2024, 8:45am

Thanks @dzenanz Are there any examples like that in practice? I’m interested to know how people are actually using the PointData of a PointSet.

Your reply suggests that PointSet::m_PointDataContainer typically has one element per point, right? In code: m_PointDataContainer->size() == m_PointsContainer->size()

dzenanz · October 1, 2024, 1:12pm

Here is an example of it being used to hold cartilage thickness:

github.com

uncbiag/OAI_analysis_2/blob/728fc41a66ae3c159ed41bde8365792d0fd8b617/oai_analysis/mesh_processing.py#L508


      
          
          if mesh_type == "FC":
              vertices = np.array(mapped_mesh.GetPoints().GetData())
              vertices[:, [1, 0]] = vertices[:, [0, 1]]
              circle, z_range = get_cylinder(vertices)
              embedded, plot_xy = get_projection_from_circle_and_vertice(vertices, circle)
          
              return embedded[:, 0], embedded[:, 1], point_data
          else:
              vertice = np.array(mapped_mesh.GetPoints().GetData())
              thickness = np.array(mapped_mesh.GetPointData().GetScalars())
          
              vertice_left = vertice[vertice[:, 2] < 50]
              index_left = np.where(vertice[:, 2] < 50)[0]
          
              vertice_right = vertice[vertice[:, 2] >= 50]
              index_right = np.where(vertice[:, 2] >= 50)[0]
          
              embedded_left = do_linear_pca(vertice_left)
              embedded_right = do_linear_pca(vertice_right)

dzenanz · October 1, 2024, 1:13pm

Earlier in the same file, being set in ITK mesh:

github.com

uncbiag/OAI_analysis_2/blob/728fc41a66ae3c159ed41bde8365792d0fd8b617/oai_analysis/mesh_processing.py#L96


      
              cell_data_numpy = np.array(
                  vtk_mesh.GetCellData().GetScalars()
              )  # .astype('float64')
          
              itk_mesh.SetPoints(itk.vector_container_from_array(points_numpy))
              itk_mesh.SetCellsArray(
                  itk.vector_container_from_array(polys_numpy),
                  itk.CommonEnums.CellGeometry_TRIANGLE_CELL,
              )
              if point_data_numpy.dtype == np.float64:
                  itk_mesh.SetPointData(itk.vector_container_from_array(point_data_numpy))
                  itk_mesh.SetCellData(itk.vector_container_from_array(cell_data_numpy))
              return itk_mesh
          
          
          # Get VTK mesh using vertices and faces array
          def get_vtk_mesh(verts, faces):
              # Add Cells
              cells = vtk.vtkCellArray()
              offsets = np.arange(0, faces.shape[0]*3 + 1, 3).astype('int')
              connectivity = faces.flatten().astype('int')

ntustison · October 7, 2024, 1:42am

It’s in the B-spline fitting filter I wrote. It can be used to create B-spline objects of an arbitrary number of parametric dimensions and an arbitrary number of data dimensions. Each point represents a location in the parametric (Point) domain and the corresponding PointData represents the data associated with that point.

So, for example, the B-spline filter is used in N4 where the scalar bias field values are the 1-D PointData at their image (i.e. parametric) locations. In contrast, the B-spline SyN image method also uses the B-spline filter where the displacement vectors are the n-D PointData at the image locations.

Niels_Dekker · October 8, 2024, 3:11pm

@dzenanz @ntustison Thank you very much for your examples. So clearly the point data doesn’t always just contain pixel values. (Even though the point data container is defined as a container of pixel values: ITK/Modules/Core/Common/include/itkDefaultStaticMeshTraits.h at master · InsightSoftwareConsortium/ITK · GitHub)

I wanted to know this because we are considering to add PointSet support to the library interface of elastix But before doing so, I would like to know how it is actually being used, in practice

dzenanz · October 8, 2024, 3:15pm

Yes, in PointSet, PixelType is a placeholder for a name. It is essentially a per-vertex arbitrary data type, at least in C++.

Niels_Dekker · October 8, 2024, 6:53pm

Thanks! Do you agree that it also happens in practice that a PointSet does not need any PointData at all? I’m considering to submit another pull request that would allow instantiating a PointSet<void> (TPixelType = void), to indicate that a certain pointset never has any PointData.

But first I hope we can figure out how to fix those Windows.Python SWIG errors from PointSetBase!

dzenanz · October 8, 2024, 7:15pm

Yes, this makes sense!