Merge pull request #18 from masadcv/add-tiovanen-method

Add  Toivenan's method, including samples, docs and updated requirements

FastGeodis/__init__.py

@@ -161,6 +161,153 @@ def signed_generalised_geodesic3d(
     )
 
 
+def generalised_geodesic2d_toivanen(
+    image: torch.Tensor, 
+    softmask: torch.Tensor, 
+    v: float, 
+    lamb: float, 
+    iter: int = 2
+):
+    r"""Computes Generalised Geodesic Distance using Toivanen's raster scanning method from:
+
+    Toivanen, Pekka J. 
+    "New geodosic distance transforms for gray-scale images." 
+    Pattern Recognition Letters 17.5 (1996): 437-450.
+    
+    For more details on generalised geodesic distance, check the following reference:
+
+    Criminisi, Antonio, Toby Sharp, and Andrew Blake.
+    "Geos: Geodesic image segmentation."
+    European Conference on Computer Vision, Berlin, Heidelberg, 2008.
+
+    The function expects input as torch.Tensor, which can be run on CPU or GPU depending on Tensor's device location
+
+    Args:
+        image: input image, can be grayscale or multiple channels.
+        softmask: softmask in range [0, 1] with seed information.
+        v: weighting factor for establishing relationship between unary and spatial distances.
+        lamb: weighting factor between 0.0 and 1.0. 0.0 returns euclidean distance, whereas 1.0 returns geodesic distance
+        iter: number of passes of the iterative distance transform method
+
+    Returns:
+        torch.Tensor with distance transform
+    """
+    return FastGeodisCpp.generalised_geodesic2d_toivanen(
+        image, softmask, v, lamb, 1 - lamb, iter
+    )
+
+
+def generalised_geodesic3d_toivanen(
+    image: torch.Tensor,
+    softmask: torch.Tensor,
+    spacing: List,
+    v: float,
+    lamb: float,
+    iter: int = 4,
+):
+    r"""Computes Generalised Geodesic Distance using Toivanen's raster scanning method from:
+
+    Toivanen, Pekka J. 
+    "New geodosic distance transforms for gray-scale images." 
+    Pattern Recognition Letters 17.5 (1996): 437-450.
+    
+    For more details on generalised geodesic distance, check the following reference:
+
+
+    Criminisi, Antonio, Toby Sharp, and Andrew Blake.
+    "Geos: Geodesic image segmentation."
+    European Conference on Computer Vision, Berlin, Heidelberg, 2008.
+
+    The function expects input as torch.Tensor, which can be run on CPU or GPU depending on Tensor's device location
+
+    Args:
+        image: input image, can be grayscale or multiple channels.
+        softmask: softmask in range [0, 1] with seed information.
+        spacing: spacing for 3D data
+        v: weighting factor for establishing relationship between unary and spatial distances.
+        lamb: weighting factor between 0.0 and 1.0. 0.0 returns euclidean distance, whereas 1.0 returns geodesic distance
+        iter: number of passes of the iterative distance transform method
+
+    Returns:
+        torch.Tensor with distance transform
+    """
+    return FastGeodisCpp.generalised_geodesic3d_toivanen(
+        image, softmask, spacing, v, lamb, 1 - lamb, iter
+    )
+
+def signed_generalised_geodesic2d_toivanen(
+    image: torch.Tensor, 
+    softmask: torch.Tensor, 
+    v: float, 
+    lamb: float, 
+    iter: int = 2
+):
+    r"""Computes Signed Generalised Geodesic Distance using Toivanen's raster scanning method from:
+
+    Toivanen, Pekka J. 
+    "New geodosic distance transforms for gray-scale images." 
+    Pattern Recognition Letters 17.5 (1996): 437-450.
+
+    For more details on generalised geodesic distance, check the following reference:
+
+    Criminisi, Antonio, Toby Sharp, and Andrew Blake.
+    "Geos: Geodesic image segmentation."
+    European Conference on Computer Vision, Berlin, Heidelberg, 2008.
+
+    The function expects input as torch.Tensor, which can be run on CPU or GPU depending on Tensor's device location
+
+    Args:
+        image: input image, can be grayscale or multiple channels.
+        softmask: softmask in range [0, 1] with seed information.
+        v: weighting factor for establishing relationship between unary and spatial distances.
+        lamb: weighting factor between 0.0 and 1.0. 0.0 returns euclidean distance, whereas 1.0 returns geodesic distance
+        iter: number of passes of the iterative distance transform method
+
+    Returns:
+        torch.Tensor with distance transform
+    """
+    return FastGeodisCpp.signed_generalised_geodesic2d_toivanen(
+        image, softmask, v, lamb, 1 - lamb, iter
+    )
+
+
+def signed_generalised_geodesic3d_toivanen(
+    image: torch.Tensor,
+    softmask: torch.Tensor,
+    spacing: List,
+    v: float,
+    lamb: float,
+    iter: int = 4,
+):
+    r"""Computes Signed Generalised Geodesic Distance using Toivanen's raster scanning method from:
+
+    Toivanen, Pekka J. 
+    "New geodosic distance transforms for gray-scale images." 
+    Pattern Recognition Letters 17.5 (1996): 437-450.
+
+    For more details on generalised geodesic distance, check the following reference:
+
+    Criminisi, Antonio, Toby Sharp, and Andrew Blake.
+    "Geos: Geodesic image segmentation."
+    European Conference on Computer Vision, Berlin, Heidelberg, 2008.
+
+    The function expects input as torch.Tensor, which can be run on CPU or GPU depending on Tensor's device location
+
+    Args:
+        image: input image, can be grayscale or multiple channels.
+        softmask: softmask in range [0, 1] with seed information.
+        spacing: spacing for 3D data
+        v: weighting factor for establishing relationship between unary and spatial distances.
+        lamb: weighting factor between 0.0 and 1.0. 0.0 returns euclidean distance, whereas 1.0 returns geodesic distance
+        iter: number of passes of the iterative distance transform method
+
+    Returns:
+        torch.Tensor with distance transform
+    """
+    return FastGeodisCpp.signed_generalised_geodesic3d_toivanen(
+        image, softmask, spacing, v, lamb, 1 - lamb, iter
+    )
+
 def GSF2d(
     image: torch.Tensor,
     softmask: torch.Tensor,
@@ -221,3 +368,74 @@ def GSF3d(
         torch.Tensor with distance transform
     """
     return FastGeodisCpp.GSF3d(image, softmask, theta, spacing, v, lamb, iter)
+
+def GSF2d_toivanen(
+    image: torch.Tensor,
+    softmask: torch.Tensor,
+    theta: float,
+    v: float,
+    lamb: float,
+    iter: int,
+):
+    r"""Computes Geodesic Symmetric Filtering (GSF) using Toivanen's raster scanning method from:
+
+    Toivanen, Pekka J. 
+    "New geodosic distance transforms for gray-scale images." 
+    Pattern Recognition Letters 17.5 (1996): 437-450.
+
+    For more details on GSF, check the following reference:
+
+    Criminisi, Antonio, Toby Sharp, and Andrew Blake.
+    "Geos: Geodesic image segmentation."
+    European Conference on Computer Vision, Berlin, Heidelberg, 2008.
+
+    The function expects input as torch.Tensor, which can be run on CPU or GPU depending on Tensor's device location
+
+    Args:
+        image: input image, can be grayscale or multiple channels.
+        softmask: softmask in range [0, 1] with seed information.
+        v: weighting factor for establishing relationship between unary and spatial distances.
+        lamb: weighting factor between 0.0 and 1.0. 0.0 returns euclidean distance, whereas 1.0 returns geodesic distance
+        iter: number of passes of the iterative distance transform method
+
+    Returns:
+        torch.Tensor with distance transform
+    """
+    return FastGeodisCpp.GSF2d_toivanen(image, softmask, theta, v, lamb, iter)
+
+
+def GSF3d_toivanen(
+    image: torch.Tensor,
+    softmask: torch.Tensor,
+    theta: float,
+    spacing: List,
+    v: float,
+    lamb: float,
+    iter: int,
+):
+    r"""Computes Geodesic Symmetric Filtering (GSF) using Toivanen's raster scanning method from:
+
+    Toivanen, Pekka J. 
+    "New geodosic distance transforms for gray-scale images." 
+    Pattern Recognition Letters 17.5 (1996): 437-450.
+    
+    For more details on GSF, check the following reference:
+
+    Criminisi, Antonio, Toby Sharp, and Andrew Blake.
+    "Geos: Geodesic image segmentation."
+    European Conference on Computer Vision, Berlin, Heidelberg, 2008.
+
+    The function expects input as torch.Tensor, which can be run on CPU or GPU depending on Tensor's device location
+
+    Args:
+        image: input image, can be grayscale or multiple channels.
+        softmask: softmask in range [0, 1] with seed information.
+        spacing: spacing for 3D data
+        v: weighting factor for establishing relationship between unary and spatial distances.
+        lamb: weighting factor between 0.0 and 1.0. 0.0 returns euclidean distance, whereas 1.0 returns geodesic distance
+        iter: number of passes of the iterative distance transform method
+
+    Returns:
+        torch.Tensor with distance transform
+    """
+    return FastGeodisCpp.GSF3d_toivanen(image, softmask, theta, spacing, v, lamb, iter)

FastGeodis/fastgeodis.cpp

@@ -144,6 +144,38 @@ torch::Tensor generalised_geodesic3d(torch::Tensor &image, const torch::Tensor &
     return generalised_geodesic3d_cpu(image, mask, spacing, v, l_grad, l_eucl, iterations);
 }
 
+torch::Tensor generalised_geodesic2d_toivanen(torch::Tensor &image, const torch::Tensor &mask, const float &v, const float &l_grad, const float &l_eucl, const int &iterations)
+{
+
+    // check input dimensions
+    check_input_dimensions(image, mask, 4);
+
+    // toivanen method is only implementable on cpu
+    check_cpu(image);    
+    check_cpu(mask);
+
+    return generalised_geodesic2d_toivanen_cpu(image, mask, v, l_grad, l_eucl, iterations);
+}
+
+torch::Tensor generalised_geodesic3d_toivanen(torch::Tensor &image, const torch::Tensor &mask, const std::vector<float> &spacing, const float &v, const float &l_grad, const float &l_eucl, const int &iterations)
+{
+    // check input dimensions
+    check_input_dimensions(image, mask, 5);
+
+    // toivanen method is only implementable on cpu
+    check_cpu(image);    
+    check_cpu(mask);
+
+    if (spacing.size() != 3)
+    {
+        throw std::invalid_argument(
+            "function only supports 3D spacing inputs, received " + std::to_string(spacing.size()));
+    }
+
+    return generalised_geodesic3d_toivanen_cpu(image, mask, spacing, v, l_grad, l_eucl, iterations);
+}
+
+
 torch::Tensor getDs2d(torch::Tensor &image, const torch::Tensor &mask, const float &v, const float &l_grad, const float &l_eucl, const int &iterations)
 {
     torch::Tensor D_M = generalised_geodesic2d(image, mask, v, l_grad, l_eucl, iterations);
@@ -186,13 +218,63 @@ torch::Tensor GSF3d(torch::Tensor &image, const torch::Tensor &mask, const float
     return Dd_Md + De_Me;
 }
 
+torch::Tensor getDs2d_toivanen(torch::Tensor &image, const torch::Tensor &mask, const float &v, const float &l_grad, const float &l_eucl, const int &iterations)
+{
+    torch::Tensor D_M = generalised_geodesic2d_toivanen(image, mask, v, l_grad, l_eucl, iterations);
+    torch::Tensor D_Mb = generalised_geodesic2d_toivanen(image, 1 - mask, v, l_grad, l_eucl, iterations);
+
+    return D_M - D_Mb;
+}
+
+torch::Tensor GSF2d_toivanen(torch::Tensor &image, const torch::Tensor &mask, const float &theta, const float &v, const float &lambda, const int &iterations)
+{
+    torch::Tensor Ds_M = getDs2d_toivanen(image, mask, v, lambda, 1 - lambda, iterations);
+
+    torch::Tensor Md = (Ds_M > theta).type_as(Ds_M);
+    torch::Tensor Me = (Ds_M > -theta).type_as(Ds_M);
+
+    torch::Tensor Dd_Md = -getDs2d_toivanen(image, 1 - Md, v, lambda, 1 - lambda, iterations);
+    torch::Tensor De_Me = getDs2d_toivanen(image, Me, v, lambda, 1 - lambda, iterations);
+
+    return Dd_Md + De_Me;
+}
+
+torch::Tensor getDs3d_toivanen(torch::Tensor &image, const torch::Tensor &mask, const std::vector<float> &spacing, const float &v, const float &l_grad, const float &l_eucl, const int &iterations)
+{
+    torch::Tensor D_M = generalised_geodesic3d_toivanen(image, mask, spacing, v, l_grad, l_eucl, iterations);
+    torch::Tensor D_Mb = generalised_geodesic3d_toivanen(image, 1 - mask, spacing, v, l_grad, l_eucl, iterations);
+
+    return D_M - D_Mb;
+}
+
+torch::Tensor GSF3d_toivanen(torch::Tensor &image, const torch::Tensor &mask, const float &theta, const std::vector<float> &spacing, const float &v, const float &lambda, const int &iterations)
+{
+    torch::Tensor Ds_M = getDs3d_toivanen(image, mask, spacing, v, lambda, 1 - lambda, iterations);
+
+    torch::Tensor Md = (Ds_M > theta).type_as(Ds_M);
+    torch::Tensor Me = (Ds_M > -theta).type_as(Ds_M);
+
+    torch::Tensor Dd_Md = -getDs3d_toivanen(image, 1 - Md, spacing, v, lambda, 1 - lambda, iterations);
+    torch::Tensor De_Me = getDs3d_toivanen(image, Me, spacing, v, lambda, 1 - lambda, iterations);
+
+    return Dd_Md + De_Me;
+}
 
 PYBIND11_MODULE(TORCH_EXTENSION_NAME, m)
 {
     m.def("generalised_geodesic2d", &generalised_geodesic2d, "Generalised Geodesic distance 2d");
     m.def("GSF2d", &GSF2d, "Geodesic Symmetric Filtering 2d");
     m.def("signed_generalised_geodesic2d", &getDs2d, "Signed Generalised Geodesic distance 2d");
+
+    m.def("generalised_geodesic2d_toivanen", &generalised_geodesic2d_toivanen, "Generalised Geodesic distance 2d using Toivanen's method");
+    m.def("GSF2d_toivanen", &GSF2d_toivanen, "Geodesic Symmetric Filtering 2d using Toivanen's method");
+    m.def("signed_generalised_geodesic2d_toivanen", &getDs2d_toivanen, "Signed Generalised Geodesic distance 2d using Toivanen's method");
+
     m.def("generalised_geodesic3d", &generalised_geodesic3d, "Generalised Geodesic distance 3d");
     m.def("GSF3d", &GSF3d, "Geodesic Symmetric Filtering 3d");
     m.def("signed_generalised_geodesic3d", &getDs3d, "Signed Generalised Geodesic distance 3d");
+
+    m.def("generalised_geodesic3d_toivanen", &generalised_geodesic3d_toivanen, "Generalised Geodesic distance 3d using Toivanen's method");
+    m.def("GSF3d_toivanen", &GSF3d_toivanen, "Geodesic Symmetric Filtering 3d using Toivanen's method");
+    m.def("signed_generalised_geodesic3d_toivanen", &getDs3d_toivanen, "Signed Generalised Geodesic distance 3d using Toivanen's method");
 }

FastGeodis/fastgeodis.h

@@ -70,6 +70,24 @@ torch::Tensor generalised_geodesic3d_cpu(
     const float &l_eucl,
     const int &iterations);
 
+torch::Tensor generalised_geodesic2d_toivanen_cpu(
+    torch::Tensor &image,
+    const torch::Tensor &mask,
+    const float &v,
+    const float &l_grad,
+    const float &l_eucl,
+    const int &iterations);
+
+torch::Tensor generalised_geodesic3d_toivanen_cpu(
+    torch::Tensor &image,
+    const torch::Tensor &mask,
+    const std::vector<float> &spacing,
+    const float &v,
+    const float &l_grad,
+    const float &l_eucl,
+    const int &iterations);
+
+
 torch::Tensor generalised_geodesic2d(
     torch::Tensor &image,
     const torch::Tensor &mask,
@@ -86,3 +104,20 @@ torch::Tensor generalised_geodesic3d(
     const float &l_grad,
     const float &l_eucl,
     const int &iterations);
+
+torch::Tensor generalised_geodesic2d_toivanen(
+    torch::Tensor &image,
+    const torch::Tensor &mask,
+    const float &v,
+    const float &l_grad,
+    const float &l_eucl,
+    const int &iterations);
+
+torch::Tensor generalised_geodesic3d_toivanen(
+    torch::Tensor &image,
+    const torch::Tensor &mask,
+    const std::vector<float> &spacing,
+    const float &v,
+    const float &l_grad,
+    const float &l_eucl,
+    const int &iterations);