Merge branch 'improve_sphere' into 'master'

Improve implementation of `SplineSurface`

See merge request uhlmann-group/python-spline-fitting-toolbox!5

.gitlab-ci.yml

@@ -9,9 +9,9 @@ before_script:
 tests :
   stage: test
   script:
-    - pip install -e . -v
-    - pip install pytest pytest-cov
-    - pytest --cov=sft/ tests/
+    - python -m pip install -e . -v
+    - python -m pip install pytest pytest-cov
+    - python -m pytest --cov=sft/ tests/
     - coverage xml
   artifacts:
     reports:

environment.yml

@@ -4,8 +4,15 @@ channels:
   - pytorch
 dependencies:
   - python=3.8
-  - cudatoolkit-dev
+  - cudatoolkit
   - pytorch
-  - scipy
-  - numpy
-  - matplotlib
+  - tensorflow
+  - tensorflow-gpu
+  - pip
+  - pip:
+    - scipy
+    - numpy
+    - matplotlib
+    - pytest
+    - scikit-image
+    - mpmath

sft/spline/basis.py

 import numpy as np
 import math
 import abc
+import mpmath
 
 # TODO: local refinement filters
 # TODO: quadratic prefilters
@@ -785,3 +786,112 @@ def multinomial(
             )
 
     return
+
+
+class CardinalExponential(Basis):
+    def __init__(self, alpha):
+        Basis.__init__(self, multigenerator=False, support=4.0)
+        self._alpha = alpha
+
+    def _ba3(self, t):
+        if 0 <= t <= 1:
+            return -(1 / (self._alpha ** 2)) * (np.cos(self._alpha * t) - 1)
+        if 1 < t <= 2:
+            return -(2 / (self._alpha ** 2)) * (
+                    np.cos(self._alpha) -
+                    0.5 * np.cos(self._alpha * (1 - t))
+                    - 0.5 * np.cos(self._alpha * (2 - t))
+            )
+        if 2 < t <= 3:
+            return -(1 / (self._alpha ** 2)) * (np.cos(self._alpha * (3 - t)) - 1)
+        return 0.0
+
+    def _ba4(self, t):
+        if 0 <= t <= 1:
+            return -(1 / (self._alpha ** 3)) * (-self._alpha * t +
+                                                np.sin(self._alpha * t))
+        elif 1 < t <= 2:
+            return -(1 / (self._alpha ** 3)) * (
+                    self._alpha * (t - 2)
+                    + 2 * self._alpha * np.cos(self._alpha) * (t - 1)
+                    + 2 * np.sin(self._alpha * (1 - t))
+                    + np.sin(self._alpha * (2 - t)))
+        if 2 < t <= 3:
+            return (1 / (self._alpha ** 3)) * (
+                    self._alpha * (t - 2) +
+                    2 * self._alpha * np.cos(self._alpha) * (t - 3) +
+                    np.sin(self._alpha * (2 - t)) +
+                    2 * np.sin(self._alpha * (3 - t))
+            )
+        if 3 < t <= 4:
+            return -(1 / (self._alpha ** 3)) * (
+                    self._alpha * (t - 4) + np.sin(self._alpha * (4 - t)))
+        return 0.0
+
+    def _ba3Prime(self, t):
+        if 0 <= t <= 1:
+            return (1 / self._alpha) * np.sin(self._alpha * t)
+        if 1 < t <= 2:
+            return (1 / self._alpha) * (np.sin(self._alpha * (1 - t)) +
+                                        np.sin(self._alpha * (2 - t)))
+        if 2 < t <= 3:
+            return -(1 / self._alpha) * np.sin(self._alpha * (3 - t))
+        return 0.0
+
+    def _ba4Prime(self, t):
+        if 0 <= t <= 1:
+            return (1 / (self._alpha ** 2)) * (
+                    2 * np.sin(0.5 * self._alpha * t) ** 2)
+        elif 1 < t <= 2:
+            return (1 / (self._alpha ** 2)) * (
+                    -2 * np.cos(self._alpha) +
+                    np.cos(self._alpha * (t - 2)) +
+                    2 * np.cos(self._alpha * (t - 1)) - 1)
+        elif 2 < t <= 3:
+            return (1 / (self._alpha ** 2)) * (
+                    2 * np.cos(self._alpha) -
+                    2 * np.cos(self._alpha * (t - 3)) +
+                    2 * (np.sin(0.5 * self._alpha * (t - 2)) ** 2))
+        elif 3 < t <= 4:
+            return -(1 / (self._alpha ** 2)) * (
+                    2 * np.sin(0.5 * self._alpha * (4 - t)) ** 2)
+        return 0.0
+
+    @property
+    def _lambda3(self):
+        # Watch out, calculus may be wrong for alphas
+        # that are not partitions of 2*pi
+        return ((0.5 * self._alpha) ** 2) * \
+               (mpmath.csc(0.5 * self._alpha) *
+                mpmath.csc(self._alpha) *
+                (1 - self._alpha * mpmath.csc(self._alpha))) / \
+               (mpmath.sec(0.5 * self._alpha)
+                - 0.5 * self._alpha * mpmath.csc(0.5 * self._alpha))
+
+    @property
+    def _lambda4(self):
+        # Watch out, calculus may be wrong for alphas
+        # that are not partitions of 2*pi
+        return ((0.5 * self._alpha) ** 3) * \
+               (mpmath.sec(0.5 * self._alpha) ** 2) / \
+               (np.tan(0.5 * self._alpha) - 0.5 * self._alpha)
+
+    def value(self, t):
+        return (self._lambda4 * self._ba4(t + 2) +
+                self._lambda3 * (self._ba3(t + 2) + self._ba3(t + 1)))
+
+    def firstDerivativeValue(self, t):
+        return (self._lambda4 * self._ba4Prime(t + 2) +
+                self._lambda3 * (self._ba3Prime(t + 2) + self._ba3Prime(t + 1)))
+
+    def secondDerivativeValue(self, x):
+        raise NotImplementedError(Basis._unimplemented_message)
+
+    def filterSymmetric(self, s):
+        raise NotImplementedError(Basis._unimplemented_message)
+
+    def filterPeriodic(self, s):
+        raise NotImplementedError(Basis._unimplemented_message)
+
+    def refinementMask(self):
+        raise NotImplementedError(Basis._unimplemented_message)

sft/spline/curve_models.py

@@ -45,7 +45,7 @@ class SplineCurve:
         else:
             raise RuntimeError(
                 "M must be greater or equal than the "
-                "spline generator support size."
+                "spline basis support size."
             )
 
         self._basis = basis
@@ -498,7 +498,7 @@ class HermiteSplineCurve(SplineCurve):
     def __init__(self, M, basis, closed):
         if not basis.multigenerator:
             raise RuntimeError(
-                "It looks like you are trying to use a single generator to "
+                "It looks like you are trying to use a single basis to "
                 "build a multigenerator spline model."
             )
 

tests/test_curve.py

 import glob
 
-import cv2
 import matplotlib.pyplot as plt
 import numpy as np
 import pytest
@@ -59,10 +58,7 @@ def test_binary_mask(img_file, plot=False):
 
     sft_arcLength = splineCurve.arcLength(t0, tf)
 
-    length_opencv = cv2.arcLength(np.float32(splineContour), True)
-
-    np.testing.assert_almost_equal(length_direct, length_opencv, decimal=4)
-    np.testing.assert_almost_equal(sft_arcLength, length_opencv, decimal=4)
+    np.testing.assert_almost_equal(sft_arcLength, length_direct, decimal=4)
 
 
 @pytest.mark.parametrize('basis_func', [B1(),

tests/test_surface.py

 #!/usr/bin/env python
 # coding: utf-8
 
-import json
-import multiprocessing
-import time
+import tensorflow as tf
+import os.path as osp
 
-import imageio
+from skimage import io
 import numpy as np
-import os.path as osp
 import pytest
 
 from sft.spline.surface_models import SplineSphere
@@ -22,7 +20,7 @@ def test_parameter_to_world():
     Mt = 5
     Ms = 5
     model = SplineSphere((Ms, Mt))
-    model.initializeSphere(radius, center)
+    model.initialize_sphere(radius, center)
 
     # Compare tangents to theoretical tangents
     t = 0.324543643
@@ -37,175 +35,173 @@ def test_parameter_to_world():
          radius * np.cos(np.pi * s) * np.sin(2.0 * np.pi * t),
          -radius * np.sin(np.pi * s)])
 
-    dds_comp = model.parametersToWorld((s * (Ms - 1), t * Mt), d=(True, False))
-    ddt_comp = model.parametersToWorld((s * (Ms - 1), t * Mt), d=(False, True))
+    dds_comp = model.parametersToWorld((s * (Ms - 1), t * Mt), ds=True, dt=False)
+    ddt_comp = model.parametersToWorld((s * (Ms - 1), t * Mt), ds=False, dt=True)
 
     np.testing.assert_almost_equal(dds_comp, dds_true)
     np.testing.assert_almost_equal(ddt_comp, ddt_true)
 
 
-def test_ray_trace():
-    M = (5, 4)
+@pytest.mark.parametrize('use_tf', [True, False])
+@pytest.mark.parametrize('Ms', np.arange(3, 5))
+@pytest.mark.parametrize('Mt', np.arange(3, 5))
+@pytest.mark.parametrize('s_rate_s', [5, 10, 20])
+@pytest.mark.parametrize('s_rate_t', [5, 10, 20])
+@pytest.mark.parametrize('radius', np.logspace(-6, 6, 7))
+@pytest.mark.parametrize('center', [(0, 0, 0), (-12.3, 37, 1337)])
+def test_initialize_from_points(use_tf, Ms, Mt, s_rate_s, s_rate_t, radius, center):
 
-    # Create a spline sphere
-    splineSurface = SplineSphere(M)
+    # We sample points to recreate the sphere
+    J = ((Ms - 1) * s_rate_s) + 1
+    I = Mt * s_rate_t
 
-    # Initialize it with a radius and a center
-    radius = 7
-    center = (10, 10, 10)
-    splineSurface.initializeSphere(radius, center)
+    sample_points = np.zeros((I * J, 3))
+    for j in range(0, J):
+        for i in range(0, I):
+            ip = i + I * j
+            point = np.array([
+                np.cos(2.0 * np.pi * i / I) * np.sin(np.pi * j / (J - 1)),
+                np.sin(2.0 * np.pi * i / I) * np.sin(np.pi * j / (J - 1)),
+                np.cos(np.pi * j / (J - 1))])
+            sample_points[ip] = center + radius * point
 
-    # Testing sample function
-    t0 = time.time()
-    res2 = splineSurface.sample((20, 20), cpuCount=multiprocessing.cpu_count())
-    print("Elapsed time: " + str(time.time() - t0))
+    if use_tf:
+        sample_points = tf.convert_to_tensor(sample_points)
 
-    surface = np.zeros((21, 21, 21), dtype=np.int8)
-    for p in res2:
-        surface[int(p[2]), int(p[1]), int(p[0])] = 255
+    surface = SplineSphere(M=(Ms, Mt))
+    surface.initialize_from_points(sample_points, sampling_rates=(s_rate_s, s_rate_t))
 
-    imageio.volwrite(osp.join(TEST_DATA, "surface.tif"), surface)
+    # Initialize sphere for the given parameters
+    sphere = SplineSphere(M=(Ms, Mt))
+    sphere.initialize_sphere(radius=radius, center=center)
 
-    # Testing initializeFromPoints function (perfect samples)
-    samplingRate = (20, 20)
+    # Check for identical internals
+    np.testing.assert_almost_equal(sphere.coefs, surface.coefs)
+    np.testing.assert_almost_equal(sphere.tans, surface.tans)
+    np.testing.assert_almost_equal(sphere.north_pole, surface.north_pole)
+    np.testing.assert_almost_equal(sphere.south_pole, surface.south_pole)
 
-    radius = 7
-    center = (15, 15, 15)
+    # Check for same sampled points
+    sample_points_from_sphere = sphere.sample(sampling_rates=(s_rate_s, s_rate_t))
+    np.testing.assert_almost_equal(sample_points, sample_points_from_sphere)
 
-    J = ((M[0] - 1) * samplingRate[0]) + 1
-    I = M[1] * samplingRate[1]
+    # Check for the internals proper types
+    assert use_tf == tf.is_tensor(surface.coefs)
+    assert use_tf == tf.is_tensor(surface.control_points)
 
-    samplePoints = np.zeros((I * J, 3))
-    for j in range(0, J):
-        for i in range(0, I):
-            ip = i + I * j
-            point =  np.array([
-                np.cos(2.0 * np.pi * i / I) * np.sin(np.pi * j / (J - 1)),
-                np.sin(2.0 * np.pi * i / I) * np.sin(np.pi * j / (J - 1)),
-                np.cos(np.pi * j / (J - 1))])
-            samplePoints[ip] = center + radius * point
-
-    surface = SplineSphere(M)
-    surface.initializeFromPoints(samplePoints, samplingRate)
-
-    # Export as JSON for interactive VTK visualisation
-    results = {
-        'center': surface.centroid().tolist(),
-        'L': np.linalg.norm(surface.coefs[0] -
-                            surface.coefs[M[1] * (M[0] - 2) + 1]),
-        'R': [[1, 0, 0], [0, 1, 0], [0, 0, 1]],
-        'longitude_count': M[0],
-        'latitude_count': M[1],
-        'controlPoints': surface.coefs.tolist(),
-        'tangentVectors': surface.tans.tolist()
-    }
-
-    data = {}
-    data[0] = results
-    with open(osp.join(TEST_DATA, "surface.json"), 'w') as fw:
-        json.dump(data, fw, indent=4, separators=(',', ': '))
-
-    def rayTrace(volume, p0, p1):
-        i0 = np.asarray(p0)
-        i1 = np.asarray(p1)
-
-        if np.sum((i0 - i1) ** 2) < 0.25:
-            return i0
-
-        if np.any(i1 < 0) or np.any(i1 >= volume.shape[::-1]):
-            i1 = edgeTrace(volume, i0, i1)
-
-        im = np.asarray((i0 + i1) / 2)
-        label0 = volume[tuple(i0.astype(int)[::-1])]
-        if label0 != 1:
-            return i0
-
-        labelm = volume[tuple(im.astype(int)[::-1])]
-        if labelm == label0:
-            return rayTrace(volume, im, i1)
-        else:
-            return rayTrace(volume, i0, im)
-
-    def edgeTrace(volume, i0, i1):
-        inf = 0
-        sup = 1
-        step = 1.0 / max(volume.shape)
-        while (sup - inf) > step:
-            factor = (inf + sup) / 2.0
-            ii = i0 + factor * (i1 - i0)
-            ii[np.logical_or(ii < 0, ii >= volume.shape[::-1])] = np.nan
-            try:
-                _ = volume[tuple(ii.astype(int)[::-1])]
-                inf = factor
-            except:
-                sup = factor
-        ii = i0 + inf * (i1 - i0)
-        return ii
-
-    # Testing initializeFromPoints function (image samples)
-    samplingRate = (20, 20)
-
-    data = imageio.volread(osp.join(TEST_DATA, "surface.tif"))
-    indices = np.where(data == 1)[::-1]
-    center = np.asarray(indices).mean(axis=1)
-
-    dS = np.array([0, 0, 1])
-    dT = np.array([[1, 0, 0], [0, 1, 0]])
-
-    poles = np.stack((rayTrace(data, center, center - (max(data.shape) * dS)),
-                      rayTrace(data, center, center + (max(data.shape) * dS))))
-    distances = np.sqrt(np.sum((poles - center) ** 2, axis=1))
-    L = distances[0]
-    Lsum = np.sum(distances)
-
-    s0 = center - L * dS
-
-    J = ((M[0] - 1) * samplingRate[0]) + 1
-    I = M[1] * samplingRate[1]
-
-    alphas = np.asarray(range(I)) * 2.0 * np.pi / I
-    if np.any(np.abs(np.cross(dT[0], dT[1]) - dS) > 1e-8):
-        alphas = -alphas
-
-    samplePoints = np.zeros((I * J, 3))
-    for j in range(0, J):
-        beta = j * np.pi / (J - 1)
-        c = s0 + (dS * ((Lsum / 2.0) - ((Lsum / 2.0) * np.cos(beta))))
 
-        for i in range(0, I):
-            ip = i + I * j
-            if (j == 0) or (j == J - 1):
-                samplePoints[ip] = c
-            else:
-                v = (np.cos(alphas[i]) * dT[0]) + (np.sin(alphas[i]) * dT[1])
-                samplePoints[ip] = rayTrace(data, c, c + (max(data.shape) * v))
-
-    surface = SplineSphere(M)
-    surface.initializeFromPoints(samplePoints, samplingRate)
-
-    # Export as JSON for interactive VTK visualisation
-    results = {
-        'center': surface.centroid().tolist(),
-        'L': np.linalg.norm(surface.coefs[0] -
-                            surface.coefs[M[1] * (M[0] - 2) + 1]),
-        'R': [[1, 0, 0], [0, 1, 0], [0, 0, 1]],
-        'longitude_count': M[0],
-        'latitude_count': M[1],
-        'controlPoints': surface.coefs.tolist(),
-        'tangentVectors': surface.tans.tolist()
-    }
-
-    data = {}
-    data[0] = results
-    with open(osp.join(TEST_DATA, "surface.json"), 'w') as fw:
-        json.dump(data, fw, indent=4, separators=(',', ': '))
+@pytest.mark.parametrize('Ms', np.arange(3, 5))
+@pytest.mark.parametrize('Mt', np.arange(3, 5))
+@pytest.mark.parametrize('N', np.logspace(1, 4, dtype=int))
+def test_sampling_N(Ms, Mt, N, radius=1337):
+
+    # Initialize sphere for the given parameters
+    sphere = SplineSphere(M=(Ms, Mt))
+    sphere.initialize_sphere(radius=radius, center=(0.0, 0.0, 0.0))
 
+    # Check for same sampled points
+    sampled_points = sphere.sample(N=N)
+    assert sampled_points.shape == (N, 3)
 
+    # We assert that point are on the sphere
+    np.testing.assert_almost_equal(
+        np.linalg.norm(sampled_points, axis=1),
+        radius * np.ones((N,)),
+        decimal=6
+    )
+
+
+@pytest.mark.parametrize('use_tf', [True, False])
+@pytest.mark.parametrize('Ms', np.arange(3, 5))
+@pytest.mark.parametrize('Mt', np.arange(3, 5))
+@pytest.mark.parametrize('radius', np.logspace(-6, 6, 7))
+@pytest.mark.parametrize('center', [(0, 0, 0), (-12.3, 37, 1337)])
+def test_initialize_from_control_points(use_tf, Ms, Mt, radius, center):
+    # We sample points to recreate the sphere
+    control_points = np.zeros(((Ms - 2) * Mt + 2, 3))
+
+    PIMs = np.pi / (Ms - 1)
+    PIMt = np.pi / Mt
+
+    scaleMs = 2.0 * (1.0 - np.cos(PIMs)) / (np.cos(PIMs / 2) - np.cos(3.0 * PIMs / 2))
+    scaleMt = 2.0 * (1.0 - np.cos(2.0 * PIMt)) / (np.cos(PIMt) - np.cos(3.0 * PIMt))
+    for l in range(1, Ms - 1):
+        for k in range(0, Mt):
+            ip = k + 1 + Mt * (l - 1)
+            theta = l * PIMs
+            phi = 2.0 * k * PIMt
+            point = np.array([
+                scaleMs * scaleMt * np.cos(phi) * np.sin(theta),
+                scaleMs * scaleMt * np.sin(phi) * np.sin(theta),
+                scaleMs * np.cos(theta)])
+            control_points[ip] = center + radius * point
+
+    # North and south poles
+    control_points[0] = center + np.array([0, 0, radius])
+    control_points[-1] = center + np.array([0, 0, -radius])
+
+    surface = SplineSphere(M=(Ms, Mt))
+
+    if use_tf:
+        control_points = tf.convert_to_tensor(control_points)
+
+    surface.control_points = control_points
+
+    tans = np.array([
+        surface.north_pole + np.array([radius * np.pi, 0, 0]),
+        surface.north_pole + np.array([0, radius * np.pi, 0]),
+        surface.south_pole + np.array([radius * np.pi, 0, 0]),
+        surface.south_pole + np.array([0, radius * np.pi, 0]),
+    ])
+    surface.tans = tans
+
+    # Initialize sphere for the given parameters
+    sphere = SplineSphere(M=(Ms, Mt))
+    sphere.initialize_sphere(radius=radius, center=center)
+
+    # Check for identical internals
+    np.testing.assert_almost_equal(sphere.north_pole, surface.north_pole)
+    np.testing.assert_almost_equal(sphere.south_pole, surface.south_pole)
+    np.testing.assert_almost_equal(sphere.coefs, surface.coefs)
+    np.testing.assert_almost_equal(sphere.tans, surface.tans)
+
+    # Check for same sampled points
+    sample_points_from_sphere = sphere.sample(sampling_rates=(10, 10))
+    sample_points_from_surface = surface.sample(sampling_rates=(10, 10))
+    np.testing.assert_almost_equal(sample_points_from_sphere,
+                                   sample_points_from_surface)
+
+    # Check for the internals proper types
+    assert use_tf == tf.is_tensor(surface.coefs)
+    assert use_tf == tf.is_tensor(surface.control_points)
+
+
+@pytest.mark.parametrize('use_tf', [True, False])
+@pytest.mark.parametrize('Ms', np.arange(3, 5))
+@pytest.mark.parametrize('Mt', np.arange(3, 5))
+@pytest.mark.parametrize('s_rate_s', [5, 10, 20])
+@pytest.mark.parametrize('s_rate_t', [5, 10, 20])
+@pytest.mark.skip(reason="To be sorted out")
+def test_initialise_from_binary_mask(use_tf, Ms, Mt, s_rate_s, s_rate_t):
+
+    data = io.imread(osp.join(TEST_DATA, "volume.tif"))
+    if use_tf:
+        data = tf.convert_to_tensor(data)
+
+    surface = SplineSphere((Ms, Mt))
+    surface.initialise_from_binary_mask(data, sampling_rate=(s_rate_s, s_rate_t))
+
+    # Check for the internals proper types
+    assert use_tf == tf.is_tensor(surface.coefs)
+    assert use_tf == tf.is_tensor(surface.control_points)
+
+
+@pytest.mark.skip(reason="Computing the winding number just takes too long.")
 @pytest.mark.parametrize('Ms', np.arange(3, 5))
 @pytest.mark.parametrize('Mt', np.arange(3, 5))
 def test_winding_number_sphere(Ms, Mt):
     unitSphere = SplineSphere((Ms, Mt))
-    unitSphere.initializeSphere(1.0, np.array([0, 0, 0]))
+    unitSphere.initialize_sphere(1.0, np.array([0, 0, 0]))
 
     assert unitSphere.isInside(np.array([0.0, 0.0, 0.0])) == 1
     assert unitSphere.isInside(np.array([0.0, 0.5, 0.0])) == 1
@@ -219,24 +215,30 @@ def test_winding_number_sphere(Ms, Mt):
 
 @pytest.mark.parametrize('Ms', np.arange(3, 5))
 @pytest.mark.parametrize('Mt', np.arange(3, 5))
-def test_value_parametrisation(Ms, Mt):
+@pytest.mark.parametrize('radius', np.logspace(-6, 6, 3))
+@pytest.mark.parametrize('center', [(0, 0, 0), (-12.3, 37, 1337)])
+def test_value_parametrisation(Ms, Mt, radius, center):
     unitSphere = SplineSphere((Ms, Mt))
-    unitSphere.initializeSphere(radius=1.0, center=np.array([0, 0, 0]))
+    unitSphere.initialize_sphere(radius=radius, center=center)
 
-    St = 20
-    Ss = 20
+    s_rate_t = 20
+    s_rate_s = 20
 
     vals_comp = []
     vals_th = []
 
-    for j in range(0, (Ss * (Ms - 1)) + 1):
-        for i in range(0, St * Mt):
-            val = unitSphere.parametersToWorld((j / Ss, i / St))
-            val_th = np.array((np.cos(i * np.pi * 2.0 / (St * Mt)) *
-                              np.sin(j * np.pi / (Ss * (Ms - 1))),
-                              np.sin(i * np.pi * 2.0 / (St * Mt)) *
-                              np.sin(j * np.pi / (Ss * (Ms - 1))),
-                              np.cos(j * np.pi / (Ss * (Ms - 1)))))
+    for j in range(0, (s_rate_s * (Ms - 1)) + 1):
+        for i in range(0, s_rate_t * Mt):
+            s = j / s_rate_s
+            t = i / s_rate_t
+            phi = t * np.pi * 2.0 / Mt
+            theta = s * np.pi / (Ms - 1)
+
+            val = unitSphere.parametersToWorld((s, t))
+            val_th = center + radius * np.array(
+                (np.cos(phi) * np.sin(theta),
+                 np.sin(phi) * np.sin(theta),
+                 np