Surface vertex classification#

In this notebook we demonstrate how to retrieve surface/vertex measurements to classify vertices on the surface. The used example data is a derived from of AV Luque and JV Veenvliet (2023) licensed CC-BY. See the creating_surfaces for how to create the surface from raw imaging data.

See also

You need to additionally install the napari-accelerated-pixel-and-object-classification plugin, e.g. by calling mamba install napari-accelerated-pixel-and-object-classification -c conda-forge from the terminal.

import napari_process_points_and_surfaces as nppas

import napari
import apoc
import numpy as np

We will be working with a simple geometry (i.e., an ellipsoid) to demonstrate the curvature property and the different settings.

surface = nppas.gastruloid()
surface

The nppas gastruloid example is derived from AV Luque and JV Veenvliet (2023) which is licensed CC-BY (https://creativecommons.org/licenses/by/4.0/legalcode) and can be downloaded from here: https://zenodo.org/record/7603081
nppas.SurfaceTuple
origin (z/y/x)[0. 0. 0.]
center of mass(z/y/x)57.561,308.175,440.144
scale(z/y/x)1.000,1.000,1.000
bounds (z/y/x)13.983...112.865
110.951...460.844
168.949...807.707
average size171.204
number of vertices3324
number of faces6643

Quantification#

We can create a table (pandas Dataframe) like this.

requested_measurements = [nppas.Quality.AREA, 
                          nppas.Quality.ASPECT_RATIO,
                          nppas.Quality.GAUSS_CURVATURE, 
                          nppas.Quality.MEAN_CURVATURE,
                          nppas.Quality.SPHERE_FITTED_CURVATURE_HECTA_VOXEL,
                          nppas.Quality.SPHERE_FITTED_CURVATURE_KILO_VOXEL,
                         ]

df = nppas.surface_quality_table(surface, requested_measurements)
df

c:\structure\code\vedo\vedo\pointcloud.py:526: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  C, residue, rank, _ = np.linalg.lstsq(A, f)  # solve AC=f
vertex_index Quality.AREA Quality.ASPECT_RATIO Quality.GAUSS_CURVATURE Quality.MEAN_CURVATURE Quality.SPHERE_FITTED_CURVATURE_HECTA_VOXEL Quality.SPHERE_FITTED_CURVATURE_KILO_VOXEL
0 0 29.997389 1.600400 0.030287 0.000490 0.000257 0.000019
1 1 46.087046 1.602183 0.011136 0.000018 0.000250 0.000019
2 2 35.886338 1.400599 0.012633 0.000142 0.000253 0.000019
3 3 22.887296 1.751932 0.036979 0.000548 0.000379 0.000019
4 4 29.952347 1.220882 0.010277 0.000047 0.000391 0.000019
... ... ... ... ... ... ... ...
3319 3319 25.079661 1.340802 0.031878 0.000606 0.000168 0.000019
3320 3320 47.213916 1.254924 0.004615 0.000003 0.000169 0.000019
3321 3321 35.964707 1.140267 0.015661 0.000198 0.000163 0.000019
3322 3322 45.673529 1.189562 0.011380 0.000100 0.000152 0.000019
3323 3323 30.105530 1.151230 0.014396 0.000190 0.000163 0.000019

3324 rows × 7 columns

To get an overview about measurements, we can summarize them:

df.describe().T
count mean std min 25% 50% 75% max
vertex_index 3324.0 1661.500000 9.597005e+02 0.000000 830.750000 1661.500000 2492.250000 3323.000000
Quality.AREA 3324.0 33.753233 1.079078e+01 5.677486 26.694735 32.956835 39.255080 125.564101
Quality.ASPECT_RATIO 3324.0 7.126810 8.990960e+01 1.038034 1.292444 1.437911 1.648299 3421.965459
Quality.GAUSS_CURVATURE 3324.0 0.016958 3.527488e-02 -1.031106 0.005509 0.013645 0.024739 0.348243
Quality.MEAN_CURVATURE 3324.0 0.000383 7.652709e-03 -0.028035 -0.000135 0.000010 0.000270 0.426018
Quality.SPHERE_FITTED_CURVATURE_HECTA_VOXEL 3324.0 0.000258 6.857515e-05 0.000152 0.000214 0.000241 0.000275 0.000545
Quality.SPHERE_FITTED_CURVATURE_KILO_VOXEL 3324.0 0.000019 3.388642e-21 0.000019 0.000019 0.000019 0.000019 0.000019

From that table, we can extract a single column as list.

sp_curvature = list(df['Quality.SPHERE_FITTED_CURVATURE_HECTA_VOXEL'])
sp_curvature[:5]

[0.0002572409483117654,
 0.00025042866718866794,
 0.0002531992227477513,
 0.00037887302875519685,
 0.00039058363169657424]

Visualizing measurements#

To visualize the measurements, we need to attach them to the surface:

quantified_surface = nppas.set_vertex_values(surface, sp_curvature)

The visualization can be customized as well, e.g. by changing the view angle and the colormap.

quantified_surface.azimuth = -135
quantified_surface
nppas.SurfaceTuple
origin (z/y/x)[0. 0. 0.]
center of mass(z/y/x)57.561,308.175,440.144
scale(z/y/x)1.000,1.000,1.000
bounds (z/y/x)13.983...112.865
110.951...460.844
168.949...807.707
average size171.204
number of vertices3324
number of faces6643
min0.00015164454247100238
max0.0005448765616528308
 
quantified_surface.cmap = 'jet'
quantified_surface
nppas.SurfaceTuple
origin (z/y/x)[0. 0. 0.]
center of mass(z/y/x)57.561,308.175,440.144
scale(z/y/x)1.000,1.000,1.000
bounds (z/y/x)13.983...112.865
110.951...460.844
168.949...807.707
average size171.204
number of vertices3324
number of faces6643
min0.00015164454247100238
max0.0005448765616528308

Interacting with surface data in Napari#

We now open Napari to interact with the data

viewer = napari.Viewer(ndisplay=3)
viewer.camera.angles = (40, -30, 55)

surface_layer = viewer.add_surface(surface, colormap='hsv')

surface_layer.properties = df.to_dict(orient='list')
surface_layer.features = df

from napari_skimage_regionprops import add_table
add_table(surface_layer, viewer)

Napari status bar display of label properties disabled because https://github.com/napari/napari/issues/5417 and https://github.com/napari/napari/issues/4342

<napari_skimage_regionprops._table.TableWidget at 0x2c7b7ebd5e0>

viewer.window.add_dock_widget(nppas.SurfaceAnnotationWidget(viewer))

<napari._qt.widgets.qt_viewer_dock_widget.QtViewerDockWidget at 0x2c7b92c8700>

Manual annotation#

Use the napari window and the Surface Annotation Widget to draw with value 2 and 3 in concave and convex region on the surface.

napari.utils.nbscreenshot(viewer)

from napari_accelerated_pixel_and_object_classification._surface_vertex_classifier import SurfaceVertexClassifierWidget
viewer.window.add_dock_widget(SurfaceVertexClassifierWidget(viewer))

<napari._qt.widgets.qt_viewer_dock_widget.QtViewerDockWidget at 0x2c7b9330d30>

Selected surface layer: surface
Selected measurements: ['Quality.AREA', 'Quality.ASPECT_RATIO', 'Quality.GAUSS_CURVATURE', 'Quality.MEAN_CURVATURE', 'Quality.SPHERE_FITTED_CURVATURE_HECTA_VOXEL', 'Quality.SPHERE_FITTED_CURVATURE_KILO_VOXEL']
selected properties       Quality.AREA  Quality.ASPECT_RATIO  Quality.GAUSS_CURVATURE  \
0        29.997389              1.600400                 0.030287   
1        46.087046              1.602183                 0.011136   
2        35.886338              1.400599                 0.012633   
3        22.887296              1.751932                 0.036979   
4        29.952347              1.220882                 0.010277   
...            ...                   ...                      ...   
3319     25.079661              1.340802                 0.031878   
3320     47.213916              1.254924                 0.004615   
3321     35.964707              1.140267                 0.015661   
3322     45.673529              1.189562                 0.011380   
3323     30.105530              1.151230                 0.014396   

      Quality.MEAN_CURVATURE  Quality.SPHERE_FITTED_CURVATURE_HECTA_VOXEL  \
0                   0.000490                                     0.000257   
1                   0.000018                                     0.000250   
2                   0.000142                                     0.000253   
3                   0.000548                                     0.000379   
4                   0.000047                                     0.000391   
...                      ...                                          ...   
3319                0.000606                                     0.000168   
3320                0.000003                                     0.000169   
3321                0.000198                                     0.000163   
3322                0.000100                                     0.000152   
3323                0.000190                                     0.000163   

      Quality.SPHERE_FITTED_CURVATURE_KILO_VOXEL  
0                                       0.000019  
1                                       0.000019  
2                                       0.000019  
3                                       0.000019  
4                                       0.000019  
...                                          ...  
3319                                    0.000019  
3320                                    0.000019  
3321                                    0.000019  
3322                                    0.000019  
3323                                    0.000019  

[3324 rows x 6 columns]
As there are no 0 in the annotated classes, the minimum is subtracted
annotated classes [0. 0. 0. ... 0. 0. 0.]
RFC predictions finished. [2. 2. 3. ... 3. 2. 3.]

Train classifier#

Next, we can use the Surface Vertex Classification widget to train a classifier.

napari.utils.nbscreenshot(viewer)

Apply classifier#

The classifier is an apoc TableRowClassifier, which can be applied to the table/DataFrame above.

table_row_classifier = apoc.TableRowClassifier(opencl_filename='table_row_classifier.cl')

result = table_row_classifier.predict(df)

result

array([1, 1, 2, ..., 2, 1, 2], dtype=uint32)

df['classification'] = result
df
vertex_index Quality.AREA Quality.ASPECT_RATIO Quality.GAUSS_CURVATURE Quality.MEAN_CURVATURE Quality.SPHERE_FITTED_CURVATURE_HECTA_VOXEL Quality.SPHERE_FITTED_CURVATURE_KILO_VOXEL _CLUSTER_ID classification
0 0 29.997389 1.600400 0.030287 0.000490 0.000257 0.000019 2.0 1
1 1 46.087046 1.602183 0.011136 0.000018 0.000250 0.000019 2.0 1
2 2 35.886338 1.400599 0.012633 0.000142 0.000253 0.000019 3.0 2
3 3 22.887296 1.751932 0.036979 0.000548 0.000379 0.000019 3.0 2
4 4 29.952347 1.220882 0.010277 0.000047 0.000391 0.000019 3.0 2
... ... ... ... ... ... ... ... ... ...
3319 3319 25.079661 1.340802 0.031878 0.000606 0.000168 0.000019 3.0 2
3320 3320 47.213916 1.254924 0.004615 0.000003 0.000169 0.000019 2.0 1
3321 3321 35.964707 1.140267 0.015661 0.000198 0.000163 0.000019 3.0 2
3322 3322 45.673529 1.189562 0.011380 0.000100 0.000152 0.000019 2.0 1
3323 3323 30.105530 1.151230 0.014396 0.000190 0.000163 0.000019 3.0 2

3324 rows × 9 columns

annotated_surface = nppas.set_vertex_values(surface, result)
annotated_surface.azimuth = -135
annotated_surface
nppas.SurfaceTuple
origin (z/y/x)[0. 0. 0.]
center of mass(z/y/x)57.561,308.175,440.144
scale(z/y/x)1.000,1.000,1.000
bounds (z/y/x)13.983...112.865
110.951...460.844
168.949...807.707
average size171.204
number of vertices3324
number of faces6643
min1
max2