tests frangi time usage

brick/processing/frangi_py_test.py

+# Copyright CNRS/Inria/UNS
+# Contributor(s): Eric Debreuve (since 2019)
+#
+# eric.debreuve@cnrs.fr
+#
+# This software is governed by the CeCILL  license under French law and
+# abiding by the rules of distribution of free software.  You can  use,
+# modify and/ or redistribute the software under the terms of the CeCILL
+# license as circulated by CEA, CNRS and INRIA at the following URL
+# "http://www.cecill.info".
+#
+# As a counterpart to the access to the source code and  rights to copy,
+# modify and redistribute granted by the license, users are provided only
+# with a limited warranty  and the software's author,  the holder of the
+# economic rights,  and the successive licensors  have only  limited
+# liability.
+#
+# In this respect, the user's attention is drawn to the risks associated
+# with loading,  using,  modifying and/or developing or reproducing the
+# software by the user in light of its specific status of free software,
+# that may mean  that it is complicated to manipulate,  and  that  also
+# therefore means  that it is reserved for developers  and  experienced
+# professionals having in-depth computer knowledge. Users are therefore
+# encouraged to load and test the software's suitability as regards their
+# requirements in conditions enabling the security of their systems and/or
+# data to be ensured and,  more generally, to use and operate it in the
+# same conditions as regards security.
+#
+# The fact that you are presently reading this means that you have had
+# knowledge of the CeCILL license and that you accept its terms.
+
+import brick.processing.frangi3 as fg_
+
+import itertools as it_
+import time as tm_
+from typing import Sequence, Tuple
+
+import matplotlib.pyplot as pl_
+import numpy as np_
+import scipy.ndimage as im_
+import skimage.draw as dw_
+import skimage.measure as ms_
+import skimage.morphology as mp_
+from mpl_toolkits.mplot3d import Axes3D
+import skimage.io as io_
+
+
+array_t = np_.ndarray
+
+
+# Free parameters
+PARALLEL_CHOICES = (False, True)
+METHOD_CHOICES = ("python", "c")
+DIFFERENTIATION_CHOICES = ("indirect", "direct")
+
+SHAPE = (200, 200, 200)
+MARGIN = 4
+
+SCALE_RANGE = (0.1, 6.1)
+SCALE_STEP = 1.0
+
+# Spiral parameters
+N_SAMPLES = 500
+ANGLE_MAX = 50.0
+HEIGHT_MAX = 2.0
+ROTATION_SPEED = 0.05
+RADIAL_SPEED = 0.05
+
+
+def Spiral(
+    n_samples: int,
+    angle_max: float,
+    height_max: float,
+    rotation_speed: float,
+    radial_speed: float,
+) -> Tuple[array_t, array_t, array_t]:
+    #
+    angles = np_.linspace(0.0, angle_max, n_samples)
+    x_coords = rotation_speed * np_.exp(radial_speed * angles) * np_.cos(angles)
+    y_coords = rotation_speed * np_.exp(radial_speed * angles) * np_.sin(angles)
+    z_coords = np_.linspace(0.0, height_max, n_samples)
+
+    return x_coords, y_coords, z_coords
+
+
+def ScaledAndShiftedAndPossiblyRounded(
+    x_coords: array_t,
+    y_coords: array_t,
+    z_coords: array_t,
+    target_bbox: Sequence[int],
+    margin: int = 0,
+    rounded: bool = False,
+) -> Tuple[array_t, array_t, array_t]:
+    #
+    if target_bbox.__len__() == 3:  # It is a shape
+        target_bbox = (
+            0,
+            target_bbox[0] - 1,
+            0,
+            target_bbox[1] - 1,
+            0,
+            target_bbox[2] - 1,
+        )
+
+    corner = np_.array((np_.min(x_coords), np_.min(y_coords), np_.min(z_coords)))
+    shape = np_.diff(
+        (
+            corner[0],
+            np_.max(x_coords),
+            corner[1],
+            np_.max(y_coords),
+            corner[2],
+            np_.max(z_coords),
+        )
+    )[0::2]
+    # Even if dealing with indices, do not do + 1
+    target_shape = np_.diff(target_bbox)[0::2] - 2 * margin
+
+    scaling = target_shape / shape
+    shift = tuple(tgt_coord + margin for tgt_coord in target_bbox[0::2])
+
+    x_coords = scaling[0] * (x_coords - corner[0]) + shift[0]
+    y_coords = scaling[1] * (y_coords - corner[1]) + shift[1]
+    z_coords = scaling[2] * (z_coords - corner[2]) + shift[2]
+
+    if rounded:
+        x_coords = np_.around(x_coords)
+        y_coords = np_.around(y_coords)
+        z_coords = np_.around(z_coords)
+
+    return x_coords, y_coords, z_coords
+
+
+def WriteCurveInVolume(
+    x_coords: array_t,
+    y_coords: array_t,
+    z_coords: array_t,
+    volume: array_t,
+    target_bbox: Sequence[int] = None,
+    thickness: int = 1,
+) -> None:
+    #
+    coords = []  # curr_x, next_x, curr_y, next_y...
+    for coord in (x_coords, y_coords, z_coords):
+        iterators = it_.tee(coord)
+        next(iterators[1], None)
+        coords.extend(iterators)
+
+    for curr_x, next_x, curr_y, next_y, curr_z, next_z in zip(*coords):
+        line = dw_.line_nd(
+            (curr_x, curr_y, curr_z),
+            (next_x, next_y, next_z),
+            endpoint=True,
+            integer=True,
+        )
+
+        volume[line] = 1
+
+    if thickness > 1:
+        selem = mp_.ball(thickness // 2, dtype=np_.bool)
+        if target_bbox is None:
+            dilated = im_.binary_dilation(volume, structure=selem).astype(
+                np_.uint8, copy=False
+            )
+            volume[...] = dilated[...]
+        else:
+            x_slice = slice(target_bbox[0], target_bbox[1] + 1)
+            y_slice = slice(target_bbox[2], target_bbox[3] + 1)
+            z_slice = slice(target_bbox[4], target_bbox[5] + 1)
+            dilated = im_.binary_dilation(
+                volume[x_slice, y_slice, z_slice], structure=selem
+            ).astype(np_.uint8, copy=False)
+            volume[x_slice, y_slice, z_slice] = dilated[...]
+
+
+def PrepareCurvePlot(
+    x_coords: array_t, y_coords: array_t, z_coords: array_t, title: str = None
+) -> None:
+    #
+    fig = pl_.figure()
+    ax = fig.add_subplot(111, projection=Axes3D.name)
+    ax.plot(x_coords, y_coords, z_coords)
+    ax.set_xlabel("x")
+    ax.set_ylabel("y")
+    ax.set_zlabel("z")
+    if title is not None:
+        ax.set_title(title)
+
+
+def PrepareIsoSurfacePlot(volume: array_t, title: str = None) -> None:
+    #
+    fig = pl_.figure()
+    ax = fig.add_subplot(111, projection=Axes3D.name)
+    verts, faces = ms_.marching_cubes_classic(
+        volume, level=0.0, spacing=(0.5, 0.5, 0.5)
+    )
+    ax.plot_trisurf(verts[:, 0], verts[:, 1], faces, verts[:, 2], cmap="Spectral", lw=1)
+    ax.set_xlabel("x")
+    ax.set_ylabel("y")
+    ax.set_zlabel("z")
+    if title is not None:
+        ax.set_title(title)
+
+
+if __name__ == "__main__":
+    #
+    print(f"STARTED: {tm_.strftime('%a, %b %d %Y @ %H:%M:%S')}")
+    start_time = tm_.time()
+
+    print("--- Volume Creation")
+
+    # volume_ = np_.zeros(SHAPE, dtype=np_.uint8)
+    # x_local, y_local, z_local = Spiral(
+    #     N_SAMPLES, ANGLE_MAX, HEIGHT_MAX, ROTATION_SPEED, RADIAL_SPEED
+    # )
+    #
+    # half_shape = tuple(dim // 2 for dim in SHAPE)
+    # thickness_ = 1
+    # for x_range in ((0, half_shape[0]), (half_shape[0] + 1, SHAPE[0] - 1)):
+    #     for y_range in ((0, half_shape[1]), (half_shape[1] + 1, SHAPE[1] - 1)):
+    #         for z_range in ((0, half_shape[2]), (half_shape[2] + 1, SHAPE[2] - 1)):
+    #             print(f"    Thickness: {thickness_}")
+    #             target_bbox_ = x_range + y_range + z_range
+    #             x_global, y_global, z_global = ScaledAndShiftedAndPossiblyRounded(
+    #                 x_local,
+    #                 y_local,
+    #                 z_local,
+    #                 target_bbox_,
+    #                 margin=MARGIN,
+    #                 rounded=True,
+    #             )
+    #             WriteCurveInVolume(
+    #                 x_global,
+    #                 y_global,
+    #                 z_global,
+    #                 volume_,
+    #                 target_bbox_,
+    #                 thickness=thickness_,
+    #             )
+    #             thickness_ += 1
+
+    data_path = 'D:\\MorganeNadal\\PyCharm\\nutrimorph\\data\\DIO_6H_6_1.70bis_2.2_3.tif'
+    volume_ = io_.imread(data_path)
+    volume_ = volume_[:,:,:,1]
+    # volume_ = volume_[:, 512:, 512:]
+
+    print("--- Display Preparation")
+
+    PrepareIsoSurfacePlot(volume_, "Original Volume")
+
+    print("--- Done")
+
+    useless_choices = (
+        # (False, "c", "direct"),
+        # (False, "c", "indirect"),
+        (True, "c", "indirect"),
+        (True, "c", "direct"),
+        # (True, "python", "direct"),
+        # (True, "python", "indirect"),
+        # (False, "python", "direct"),
+        # (False, "python", "indirect"),
+    )
+    prm_choices = []
+    for parallel_prm in PARALLEL_CHOICES:
+        for method_prm in METHOD_CHOICES:
+            for differentiation_prm in DIFFERENTIATION_CHOICES:
+                choice = (parallel_prm, method_prm, differentiation_prm)
+                if choice not in useless_choices:
+                    prm_choices.append(choice)
+
+    for idx, parameters in enumerate(prm_choices, start=1):
+        prm_as_str = "/".join(elm.__str__().capitalize() for elm in parameters)
+        prm_as_str = prm_as_str.replace("False", "Sequential")
+        prm_as_str = prm_as_str.replace("True", "Parallel")
+
+        elapsed_time = tm_.gmtime(tm_.time() - start_time)
+        print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
+
+        print(f"--- Frangi Enhancement {idx} of {prm_choices.__len__()}: {prm_as_str}")
+
+        enhanced, scale_map = fg_.FrangiEnhancement(
+            volume_,
+            scale_range=SCALE_RANGE,
+            scale_step=SCALE_STEP,
+            alpha=0.5,
+            beta=0.5,
+            frangi_c=500.0,
+            in_parallel=parameters[0],
+            method=parameters[1],
+            differentiation_mode=parameters[2],
+        )
+        elapsed_time = tm_.gmtime(tm_.time() - start_time)
+        print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
+
+        print("--- Display Preparation")
+
+        PrepareIsoSurfacePlot(enhanced, f"Enhanced Volume {prm_as_str}")
+        fig_ = pl_.figure()
+        ax_ = fig_.add_subplot(111)
+        ax_.matshow(scale_map.max(axis=2))
+
+        print("--- Done")
+        elapsed_time = tm_.gmtime(tm_.time() - start_time)
+        print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
+
+    elapsed_time = tm_.gmtime(tm_.time() - start_time)
+    print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
+    print(f"DONE: {tm_.strftime('%a, %b %d %Y @ %H:%M:%S')}")
+
+    pl_.show()

nutrimorph.py

@@ -118,7 +118,7 @@ image = in_.ImageVerification(image, channel)
 # iv_.image_verification(image, channel)  # -> PySide2 user interface  # TODO: must return the modified image!
 # /!\ conflicts between some versions of PySide2 and Python3
 
-image = image[:, 800:, 800:]  # 512 # 562 # Just for development
+# image = image[:, 800:, 800:]  # 512 # 562 # Just for development
 img_shape = image.shape
 
 #
@@ -197,6 +197,9 @@ enhanced_ext, ext_scales = extension_t.EnhancedForDetection(
     method,
     in_parallel=in_parallel)
 
+elapsed_time = tm_.gmtime(tm_.time() - start_time)
+print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
+
 # Creation of the enhanced maps
 ext_nfo["coarse_map"] = extension_t.CoarseMap(enhanced_ext, ext_low_c, ext_high_c, ext_selem_pixel_c)
 ext_nfo["coarse_map"], ext_lmp = extension_t.FilteredCoarseMap(ext_nfo["coarse_map"], ext_min_area_c)

parameters.ini

@@ -80,7 +80,8 @@ max_straight_sq_dist_c     :    0.24050024 * (30 ** 2)
 max_weighted_length_c      :    0.24050024 * 20.0
 
 [Frangi]
-scale_range                :    0.024,0.745
+scale_range                :    0.024,1,466
+#0.745
 # 0.1,3.1
 ; write min,max scale range
 scale_step                 :    0.24
@@ -89,7 +90,10 @@ alpha                      :    0.19
 # 0.8
 beta                       :    0.12
 # 0.5
-frangi_c                   :    3
+frangi_c                   :    0.5
+; For an image ranging from 0 to 65280 : frangi_c btw 0 and 1.
+; For an image from 200 to 65280 : frangi_c btw 1 and 10
+; highly depends on the size of the extensions in the image and on their intensity
 # 60.12506 # 500.0
 bright_on_dark             :    True
 method                     :    c