Add length multiplier to normals

.github/workflows/test_and_deploy.yml

@@ -19,8 +19,8 @@ jobs:
     runs-on: ${{ matrix.platform }}
     strategy:
       matrix:
-        platform: [windows-latest, macos-latest] # ubuntu-latest is failing anyway
-        python-version: ['3.8', '3.9', '3.10']
+        platform: [windows-latest, macos-latest, ubuntu-latest] # ubuntu-latest is failing anyway
+        python-version: ['3.9', '3.10']
 
     steps:
       - uses: actions/checkout@v3

docs/_toc.yml

@@ -26,6 +26,7 @@ parts:
         - file: 01_code_usage/02_example_workflows/demo_surface_tracing
         - file: 01_code_usage/02_example_workflows/demo_spherical_harmonics_code
         - file: 01_code_usage/02_example_workflows/demo_timelapse_processing
+        - file: 01_code_usage/02_example_workflows/demo_measure_intensity_on_surface
 
     - file: 02_interactive_usage/interactive_usage
       sections:
@@ -59,3 +60,5 @@ parts:
     chapters:
     - file: 06_API/approximation
     - file: 06_API/vectors
+    - file: 06_API/measurements
+    - file: 06_API/reconstruction

src/napari_stress/_measurements/__init__.py

@@ -23,3 +23,4 @@
 from .stresses import anisotropic_stress, tissue_stress_tensor, maximal_tissue_anisotropy
 from .toolbox import comprehensive_analysis
 from .measurements import distance_to_k_nearest_neighbors
+from .intensity import sample_intensity_along_vector, measure_intensity_on_surface

src/napari_stress/_measurements/intensity.py

@@ -0,0 +1,237 @@
+from napari.types import LayerDataTuple
+from napari_tools_menu import register_function
+import pandas as pd
+import numpy as np
+from .._utils.frame_by_frame import frame_by_frame
+
+
+@register_function(
+    menu="Measurement > Measure intensities on surface (n-STRESS)")
+@frame_by_frame
+def _measure_intensity_on_surface(
+    surface: "napari.types.SurfaceData",
+    image: "napari.types.ImageData",
+    measurement_range: float = 1.0,
+    sampling_distance: float = 0.5,
+    center: bool = True,
+    interpolation_method: str = "linear"
+) -> LayerDataTuple:
+    """
+    Measure intensity on surface.
+
+    This function exposes the `measure_intensity_on_surface` function
+    to napari. It returns a tuple of (surface, properties, 'surface').
+
+    Parameters
+    ----------
+    surface : SurfaceData
+    image : ImageData
+    measurement_range : float, optional
+        Range of measurement, by default 1.0. This determines in which
+        range around the surface the intensity is measured.
+    sampling_distance : float, optional
+        Distance between samples, by default 0.5. This determines how many samples
+        are taken in the measurement range. Needs to be smaller than measurement_range.
+    center : bool, optional
+        Center normal vectors on surface, by default True. If set to False, the
+        normal vectors point away from the surface so that the intensity is measured
+        only on one side of the surface.
+    interpolation_method : str, optional
+        Interpolation method, by default "linear"
+
+    Returns
+    -------
+    LayerDataTuple
+        Tuple of (surface, properties, 'surface')
+    """
+    intensities = measure_intensity_on_surface(surface,
+                                               image,
+                                               measurement_range=measurement_range,
+                                               sampling_distance=sampling_distance,
+                                               center=center,
+                                               interpolation_method=interpolation_method)
+
+    intensity_mean = intensities.mean(axis=1)
+    intensity_std = intensities.std(axis=1)
+    intensity_max = intensities.max(axis=1)
+    intensity_min = intensities.min(axis=1)
+
+    intensities['intensity_mean'] = intensity_mean
+    intensities['intensity_std'] = intensity_std
+    intensities['intensity_max'] = intensity_max
+    intensities['intensity_min'] = intensity_min
+
+    properties = {
+        'metadata': {'features': intensities},
+        'colormap': 'inferno'
+    }
+
+    surface = list(surface)
+    surface[2] = intensities['intensity_max'].values
+
+    return (surface, properties, 'surface')   
+
+
+def measure_intensity_on_surface(
+    surface: "napari.types.SurfaceData",
+    image: "napari.types.ImageData",
+    measurement_range: float = 1.0,
+    sampling_distance: float = 0.5,
+    center: bool = True,
+    interpolation_method: str = "linear"
+) -> pd.DataFrame:
+    """
+    Measure intensity on surface.
+
+    Parameters
+    ----------
+    surface : SurfaceData
+    image : ImageData
+    measurement_range : float, optional
+        Range of measurement, by default 1.0. This determines in which
+        range around the surface the intensity is measured.
+    sampling_distance : float, optional
+        Distance between samples, by default 0.5. This determines how many samples
+        are taken in the measurement range. Needs to be smaller than measurement_range.
+    center : bool, optional
+        Center normal vectors on surface, by default True. If set to False, the
+        normal vectors point away from the surface so that the intensity is measured
+        only on one side of the surface.
+    interpolation_method : str, optional
+        Interpolation method, by default "linear"
+
+    Returns
+    -------
+    pd.DataFrame
+        Intensity values
+    """
+    from .._vectors import normal_vectors_on_surface
+
+    normal_vectors = normal_vectors_on_surface(surface,
+                                               length_multiplier=measurement_range,
+                                               center=True)
+    intensities = sample_intensity_along_vector(normal_vectors,
+                                                image,
+                                                sampling_distance=sampling_distance,
+                                                interpolation_method=interpolation_method)
+
+    return intensities
+
+
+@register_function(
+    menu="Measurement > Measure intensities along normals (n-STRESS)")
+@frame_by_frame
+def _sample_intensity_along_vector(
+    sample_vectors: "napari.types.VectorsData",
+    image: "napari.types.ImageData",
+    sampling_distance: float = 1.0,
+    interpolation_method: str = "linear",
+) -> LayerDataTuple:
+    """
+    Sample intensity along sample_vectors of equal length.
+
+    This function exposes the `sample_intensity_along_vector` function
+    to napari. It returns a tuple of (sample_vectors, properties, 'vectors').
+    It requires the vectors to be of equal length.
+
+    Parameters
+    ----------
+    sample_vectors : VectorsData
+        sample_vectors along which to measure intensity
+    image : ImageData
+        Image to sample
+    sampling_distance : float, optional
+        Distance between samples, by default 1.0
+    interpolation_method : str, optional
+        Interpolation method, by default "linear"
+
+    Returns
+    -------
+    LayerDataTuple
+        Tuple of (sample_vectors, properties, 'vectors')
+    """
+    intensities = sample_intensity_along_vector(
+        sample_vectors,
+        image,
+        sampling_distance=sampling_distance,
+        interpolation_method=interpolation_method,
+    )
+    intensity_mean = intensities.mean(axis=1)
+    intensity_std = intensities.std(axis=1)
+    intensity_max = intensities.max(axis=1)
+    intensity_min = intensities.min(axis=1)
+
+    intensities['intensity_mean'] = intensity_mean
+    intensities['intensity_std'] = intensity_std
+    intensities['intensity_max'] = intensity_max
+    intensities['intensity_min'] = intensity_min
+
+    properties = {
+        'features': intensities,
+        'edge_color': 'intensity_max',
+        'edge_width': 1,
+        'edge_colormap': 'inferno'
+    }
+    return (sample_vectors, properties, 'vectors')
+
+
+def sample_intensity_along_vector(
+    sample_vectors: "napari.types.sample_vectorsData",
+    image: "napari.types.ImageData",
+    sampling_distance: float = 1.0,
+    interpolation_method: str = "linear",
+) -> pd.DataFrame:
+    """
+    Sample intensity along sample_vectors of equal length.
+
+    Args:
+        sample_vectors (napari.types.sample_vectorsData): sample_vectors
+            along which to measure intensity
+        image (napari.types.ImageData): Image to sample
+        sampling_distance (float, optional): Distance between samples.
+            Defaults to 1.0.
+        interpolation_method (str, optional): Interpolation method.
+            Defaults to 'linear'.
+    """
+    from scipy.interpolate import RegularGridInterpolator
+
+    # Create coords for interpolator
+    X1 = np.arange(0, image.shape[0], 1)
+    X2 = np.arange(0, image.shape[1], 1)
+    X3 = np.arange(0, image.shape[2], 1)
+    interpolator = RegularGridInterpolator(
+        (X1, X2, X3),
+        image,
+        bounds_error=False,
+        fill_value=np.nan,
+        method=interpolation_method,
+    )
+
+    # get start point and unit normal vector
+    start_points = sample_vectors[:, 0]
+    lengths = np.linalg.norm(sample_vectors[:, 1], axis=1)[:, np.newaxis]
+    unit_vector = sample_vectors[:, 1] / lengths
+
+    # calculate number of steps
+    steps = np.round(lengths / sampling_distance).mean().astype(int)
+
+    # sample intensity
+    intensity = np.zeros((len(start_points), steps))
+
+    coordinates_along_sample_vectors = []
+    for idx in range(len(start_points)):
+        _directions = [
+            sampling_distance * unit_vector[idx] * i for i in np.arange(steps)
+        ]
+        _directions = np.stack(_directions)
+
+        _sample_vectors = start_points[idx] + _directions
+        coordinates_along_sample_vectors.append(_sample_vectors)
+
+    for idx in range(len(start_points)):
+        intensity[idx, :] = interpolator(coordinates_along_sample_vectors[idx])
+
+    intensity = pd.DataFrame(np.stack(intensity))
+    intensity.columns = [f"step_{idx}" for idx in range(steps)]
+
+    return intensity

src/napari_stress/_reconstruction/refine_surfaces.py

@@ -97,6 +97,7 @@ def trace_refinement_of_surface(
 
     """
     from .. import _vectors as vectors
+    from .._measurements.intensity import sample_intensity_along_vector
     from .._measurements.measurements import distance_to_k_nearest_neighbors
 
     if isinstance(selected_fit_type, str):
@@ -125,7 +126,7 @@ def trace_refinement_of_surface(
     vector_step = trace_vectors / n_samples
 
     # measure intensity along the vectors
-    intensity_along_vector = vectors.sample_intensity_along_vector(
+    intensity_along_vector = sample_intensity_along_vector(
         np.stack([start_points, trace_vectors]).transpose((1, 0, 2)),
         intensity_image,
         sampling_distance=sampling_distance,

src/napari_stress/_tests/test_measurements.py

@@ -2,6 +2,33 @@
 import numpy as np
 
 
+def test_intensity_measurement_on_normals():
+    from napari_stress import vectors, measurements, sample_data, frame_by_frame
+    import napari_segment_blobs_and_things_with_membranes as nsbatwm
+    import napari_process_points_and_surfaces as nppas
+
+    droplet = sample_data.get_droplet_4d()[0][0]
+    droplet_rescaled = frame_by_frame(nsbatwm.rescale)(droplet, scale_x=1, scale_y=1, scale_z=2)
+    droplet_binary = frame_by_frame(nsbatwm.threshold_otsu)(droplet_rescaled)
+
+    surface = frame_by_frame(nppas.label_to_surface)(droplet_binary, 1)
+    surface_smooth = frame_by_frame(nppas.smooth_surface)(surface)
+    surface_decimated = frame_by_frame(nppas.decimate_quadric)(surface_smooth,
+                                                               number_of_vertices = 1000)
+
+    # measure intensity on normal vectors
+    normals = vectors.normal_vectors_on_surface(surface_decimated, length_multiplier=5, center=True)
+    vectors_LDtuple = measurements.intensity._sample_intensity_along_vector(normals, droplet_rescaled)
+    assert len(vectors_LDtuple[0]) == len(normals)
+    assert 'intensity_mean' in vectors_LDtuple[1]['features'].keys()
+
+    # measure intensity directly on surface
+    intensities = measurements.intensity._measure_intensity_on_surface(surface_decimated,
+                                                                       droplet_rescaled)
+    assert len(intensities[0][2]) == len(surface_decimated[0])
+    assert 'intensity_mean' in intensities[1]['metadata']['features'][0].keys()
+
+
 def test_mean_curvature_on_ellipsoid():
     """Test that the mean curvature is computed correctly."""
     from napari_stress import reconstruction, sample_data, measurements

src/napari_stress/_tests/test_reconstruction.py

@@ -56,26 +56,6 @@ def test_quadrature_point_reconstuction(make_napari_viewer):
     reconstruction.reconstruct_surface_from_quadrature_points(lebedev_points)
 
 
-def test_vector_tools():
-    from napari_stress import vectors
-    import vedo
-    import numpy as np
-
-    sphere = vedo.Sphere(r=10, pos=(50, 50, 50))
-    image = np.random.rand(100, 100, 100)
-    sampling_distance = 0.25
-
-    vectors.normal_vectors_on_pointcloud(sphere.points())
-    normal_vectors = vectors.normal_vectors_on_surface(
-        (sphere.points(),
-         np.asarray(sphere.faces())))
-    df_intensity = vectors.sample_intensity_along_vector(
-        normal_vectors, image, sampling_distance=sampling_distance)
-
-    assert df_intensity.shape[0] == normal_vectors.shape[0]
-    assert df_intensity.shape[1] == 1/sampling_distance
-
-
 def test_surface_tracing():
     from napari_stress import reconstruction
     from skimage import filters, morphology

src/napari_stress/_tests/test_vectors.py

@@ -1,3 +1,6 @@
+import numpy as np
+
+
 def test_normal_vectors_on_pointcloud():
     import vedo
     from napari_stress import vectors
@@ -6,6 +9,19 @@ def test_normal_vectors_on_pointcloud():
     normal_vectors = vectors.normal_vectors_on_pointcloud(sphere)
     assert len(normal_vectors) == len(sphere)
 
+    # test length multiplier
+    multiplier = 2
+    normal_vectors = vectors.normal_vectors_on_pointcloud(
+        sphere, length_multiplier=multiplier)
+    vector_legth = np.linalg.norm(normal_vectors[:, 1], axis=1)
+    assert np.allclose(vector_legth, multiplier)
+
+    # if centering is False, base points should be the same as input points
+    center = False
+    normal_vectors = vectors.normal_vectors_on_pointcloud(
+        sphere, center=center)
+    assert np.allclose(normal_vectors[:, 0], sphere)
+
 
 def test_normal_vectors_on_surface():
     import vedo
@@ -17,6 +33,19 @@ def test_normal_vectors_on_surface():
     normal_vectors = vectors.normal_vectors_on_surface(sphere)
     assert len(normal_vectors) == len(sphere[0])
 
+    # test length multiplier
+    multiplier = 2
+    normal_vectors = vectors.normal_vectors_on_surface(
+        sphere, length_multiplier=multiplier)
+    vector_legth = np.linalg.norm(normal_vectors[:, 1], axis=1)
+    assert np.allclose(vector_legth, multiplier)
+
+    # Test centering
+    center = False
+    normal_vectors = vectors.normal_vectors_on_surface(
+        sphere, center=center)
+    assert np.allclose(normal_vectors[:, 0], sphere[0])
+
 
 def test_pairwise_point_distances():
     import vedo