Repository: rapidsai/cucim (language: python)

Issue:
Adding scikit-image' Blob detection
<!--

Thank you for contributing to cuCIM :)

Here are some guidelines to help the review process go smoothly.

1. Please write a description in this text box of the changes that are being
   made.

2. Please ensure that you have written units tests for the changes made/features
   added.

3. If you are closing an issue please use one of the automatic closing words as
   noted here: https://help.github.com/articles/closing-issues-using-keywords/

4. If your pull request is not ready for review but you want to make use of the
   continuous integration testing facilities please label it with `[WIP]`.

5. If your pull request is ready to be reviewed without requiring additional
   work on top of it, then remove the `[WIP]` label (if present) and replace
   it with `[REVIEW]`. If assistance is required to complete the functionality,
   for example when the C/C++ code of a feature is complete but Python bindings
   are still required, then add the label `[HELP-REQ]` so that others can triage
   and assist. The additional changes then can be implemented on top of the
   same PR. If the assistance is done by members of the rapidsAI team, then no
   additional actions are required by the creator of the original PR for this,
   otherwise the original author of the PR needs to give permission to the
   person(s) assisting to commit to their personal fork of the project. If that
   doesn't happen then a new PR based on the code of the original PR can be
   opened by the person assisting, which then will be the PR that will be
   merged.

6. Once all work has been done and review has taken place please do not add
   features or make changes out of the scope of those requested by the reviewer
   (doing this just add delays as already reviewed code ends up having to be
   re-reviewed/it is hard to tell what is new etc!). Further, please do not
   rebase your branch on main/force push/rewrite history, doing any of these
   causes the context of any comments made by reviewers to be lost. If
   conflicts occur against main they should be resolved by merging main
   into the branch used for making the pull request.

Many thanks in advance for your cooperation!

-->
Description of the pull request can be find in #398, here a short summary:
This pull request will implement the scikit-image' blob detection.
The following helper functions moved to into blob.cu:
    - _compute_disk_overlap
    - _compute_sphere_overlap
    - _blob_overlap
    - _prune_blobs
Also, some functions to calculation the approximated hessian determinant matrix moved into _hessian_det_appx.cu (in skimage this is implemented with the help of pythran):
    - _clip
    - _integ
    - _hessian_matrix_det

Credits also to @annendominik, who helped with this implementation.

This PR will close #398 

Unit tests have been added; but because I am on a windows system; I did not test it. However; I have tested all blob detection with this script and basically like it works:

```python
from math import sqrt
import time
from skimage import data
from skimage.feature import blob_dog, blob_log, blob_doh
from skimage.color import rgb2gray

import matplotlib.pyplot as plt

print("------- skimage implementation -------")
image = data.hubble_deep_field()  # [0:500, 0:500]
image_gray = rgb2gray(image)

_start_time = time.time()
blobs_log = blob_log(image_gray, max_sigma=30, num_sigma=10, threshold=.1)
print("elapsed time blobs_log: %3.3f seconds" % (time.time() - _start_time))

# Compute radii in the 3rd column.
blobs_log[:, 2] = blobs_log[:, 2] * sqrt(2)

_start_time = time.time()
blobs_dog = blob_dog(image_gray, max_sigma=30, threshold=.1)
print("elapsed time blob_dog: %3.3f seconds" % (time.time() - _start_time))
blobs_dog[:, 2] = blobs_dog[:, 2] * sqrt(2)

_start_time = time.time()
blobs_doh = blob_doh(image_gray, max_sigma=30, threshold=.01)
print("elapsed time blobs_doh: %3.3f seconds" % (time.time() - _start_time))

blobs_list = [blobs_log, blobs_dog, blobs_doh]
colors = ['yellow', 'lime', 'red']
titles = ['Laplacian of Gaussian', 'Difference of Gaussian',
          'Determinant of Hessian']
sequence = zip(blobs_list, colors, titles)

fig, axes = plt.subplots(1, 3, figsize=(9, 3), sharex=True, sharey=True)
ax = axes.ravel()

for idx, (blobs, color, title) in enumerate(sequence):
    ax[idx].set_title(title)
    ax[idx].imshow(image)
    for blob in blobs:
        y, x, r = blob
        c = plt.Circle((x, y), r, color=color, linewidth=2, fill=False)
        ax[idx].add_patch(c)
    ax[idx].set_axis_off()
plt.tight_layout()
print()


print("------- cucim.skimage implementation -------")
from cucim.skimage.feature import blob_dog, blob_log, blob_doh
from cucim.skimage.color import rgb2gray
import cupy as cp

gpu_sync = cp.cuda.stream.get_current_stream().synchronize

image = data.hubble_deep_field()  # [0:500, 0:500]
image_gray = rgb2gray(cp.asarray(image))
print("dtype: %s" % str(image_gray.dtype))

blob_log(image_gray, max_sigma=30, num_sigma=10, threshold=.1)  # cache function
gpu_sync()
_start_time = time.time()
blobs_log = blob_log(image_gray, max_sigma=30, num_sigma=10, threshold=.1)
gpu_sync()
print("elapsed time blobs_log: %3.3f seconds" % (time.time() - _start_time))

# Compute radii in the 3rd column.
blobs_log[:, 2] = blobs_log[:, 2] * sqrt(2)

blob_dog(image_gray, max_sigma=30, threshold=.1)
gpu_sync()
_start_time = time.time()
blobs_dog = blob_dog(image_gray, max_sigma=30, threshold=.1)
gpu_sync()
print("elapsed time blob_dog: %3.3f seconds" % (time.time() - _start_time))
blobs_dog[:, 2] = blobs_dog[:, 2] * sqrt(2)

blob_doh(image_gray, max_sigma=30, threshold=.01)
gpu_sync()
_start_time = time.time()
blobs_doh = blob_doh(image_gray, max_sigma=30, threshold=.01)
gpu_sync()
print("elapsed time blobs_doh: %3.3f seconds" % (time.time() - _start_time))

image_gray = image_gray.astype(cp.float64)
print("dtype: %s" % str(image_gray.dtype))

blob_log(image_gray, max_sigma=30, num_sigma=10, threshold=.1)  # cache function
gpu_sync()
_start_time = time.time()
blobs_log = blob_log(image_gray, max_sigma=30, num_sigma=10, threshold=.1)
gpu_sync()
print("elapsed time blobs_log: %3.3f seconds" % (time.time() - _start_time))

# Compute radii in the 3rd column.
blobs_log[:, 2] = blobs_log[:, 2] * sqrt(2)

try:
    blob_dog(image_gray, max_sigma=30, threshold=.1)
    gpu_sync()
    _start_time = time.time()
    blobs_dog = blob_dog(image_gray, max_sigma=30, threshold=.1)
    gpu_sync()
    print("elapsed time blob_dog: %3.3f seconds" % (time.time() - _start_time))
    blobs_dog[:, 2] = blobs_dog[:, 2] * sqrt(2)
except Exception as e:
    import sys
    print(e, file=sys.stderr)

blob_doh(image_gray, max_sigma=30, threshold=.01)
gpu_sync()
_start_time = time.time()
blobs_doh = blob_doh(image_gray, max_sigma=30, threshold=.01)
gpu_sync()
print("elapsed time blobs_doh: %3.3f seconds" % (time.time() - _start_time))

blobs_list = [cp.asnumpy(blobs_log), cp.asnumpy(blobs_dog), cp.asnumpy(blobs_doh)]
colors = ['yellow', 'lime', 'red']
titles = ['Laplacian of Gaussian', 'Difference of Gaussian',
          'Determinant of Hessian']
sequence = zip(blobs_list, colors, titles)

fig_cucim, axes_cucim = plt.subplots(1, 3, figsize=(9, 3), sharex=True, sharey=True)
ax = axes_cucim.ravel()

for idx, (blobs, color, title) in enumerate(sequence):
    ax[idx].set_title(title)
    ax[idx].imshow(image)
    for blob in blobs:
        y, x, r = blob
        c = plt.Circle((x, y), r, color=color, linewidth=2, fill=False)
        ax[idx].add_patch(c)
    ax[idx].set_axis_off()
plt.tight_layout()

plt.show()
```

However, for dtype=`cp.float64` I got the following error:
```
Traceback (most recent call last):
  File "C:/Users/RuneMonzel/autolucid/01_python_src/dev/annen/test.py", line 100, in <module>
    blob_dog(image_gray, max_sigma=30, threshold=.1)
  File "c:\users\runemonzel\cucim\python\cucim\src\cucim\skimage\feature\blob.py", line 320, in blob_dog
    gaussian_images = [gaussian(image, s, mode='reflect') for s in sigma_list]
  File "c:\users\runemonzel\cucim\python\cucim\src\cucim\skimage\feature\blob.py", line 320, in <listcomp>
    gaussian_images = [gaussian(image, s, mode='reflect') for s in sigma_list]
  File "c:\users\runemonzel\cucim\python\cucim\src\cucim\skimage\_shared\utils.py", line 294, in fixed_func
    return func(*args, **kwargs)
  File "c:\users\runemonzel\cucim\python\cucim\src\cucim\skimage\_shared\filters.py", line 141, in gaussian
    ndi.gaussian_filter(image, sigma, output=output, mode=mode, cval=cval,
  File "c:\users\runemonzel\cucim\python\cucim\src\cucim\skimage\_vendored\_ndimage_filters.py", line 419, in gaussian_filter
    return _run_1d_correlates(input, list(zip(sigmas, orders)), get, output,
  File "c:\users\runemonzel\cucim\python\cucim\src\cucim\skimage\_vendored\_ndimage_filters.py", line 248, in _run_1d_correlates
    return _filters_core._run_1d_filters(
  File "c:\users\runemonzel\cucim\python\cucim\src\cucim\skimage\_vendored\_ndimage_filters_core.py", line 100, in _run_1d_filters
    fltr(input, arg, axis, temp, mode, cval, origin, **filter_kwargs)
  File "c:\users\runemonzel\cucim\python\cucim\src\cucim\skimage\_vendored\_ndimage_filters.py", line 112, in correlate1d
    return _correlate_or_convolve1d(
  File "c:\users\runemonzel\cucim\python\cucim\src\cucim\skimage\_vendored\_ndimage_filters.py", line 204, in _correlate_or_convolve1d
    out = _shmem_convolve1d(input, weights, axis=axis, output=output,
  File "c:\users\runemonzel\cucim\python\cucim\src\cucim\skimage\filters\_separable_filtering.py", line 967, in _shmem_convolve1d
    conv_axis_kernel, block, patch_per_block = _get_separable_conv_kernel(
  File "cupy\_util.pyx", line 67, in cupy._util.memoize.decorator.ret
  File "c:\users\runemonzel\cucim\python\cucim\src\cucim\skimage\filters\_separable_filtering.py", line 858, in _get_separable_conv_kernel
    return m.get_function(func_name), block, patch_per_block
  File "cupy\_core\raw.pyx", line 485, in cupy._core.raw.RawModule.get_function
  File "cupy\_core\raw.pyx", line 96, in cupy._core.raw.RawKernel.kernel.__get__
  File "cupy\_core\raw.pyx", line 113, in cupy._core.raw.RawKernel._kernel
  File "cupy\_util.pyx", line 67, in cupy._util.memoize.decorator.ret
  File "cupy\_core\raw.pyx", line 547, in cupy._core.raw._get_raw_module
  File "cupy\_core\core.pyx", line 2081, in cupy._core.core.compile_with_cache
  File "cupy\_core\core.pyx", line 2141, in cupy._core.core.compile_with_cache
  File "C:\Users\RuneMonzel\python38_envs\cucim\lib\site-packages\cupy\cuda\compiler.py", line 492, in _compile_module_with_cache
    return _compile_with_cache_cuda(
  File "C:\Users\RuneMonzel\python38_envs\cucim\lib\site-packages\cupy\cuda\compiler.py", line 561, in _compile_with_cache_cuda
    mod.load(cubin)
  File "cupy\cuda\function.pyx", line 264, in cupy.cuda.function.Module.load
  File "cupy\cuda\function.pyx", line 266, in cupy.cuda.function.Module.load
  File "cupy_backends\cuda\api\driver.pyx", line 210, in cupy_backends.cuda.api.driver.moduleLoadData
  File "cupy_backends\cuda\api\driver.pyx", line 60, in cupy_backends.cuda.api.driver.check_status
cupy_backends.cuda.api.driver.CUDADriverError: CUDA_ERROR_INVALID_PTX: a PTX JIT compilation failed

Process finished with exit code 1
```

@grlee77 
Any thoughts on that??



[FEA] Implementing Scikit-Image's Blob Detection
Hello to the cucim dev team,

@annendominik and I are currently implementing Scikit's image blob detection and want to add it to this beautiful lib `cucim`.

**Is your feature request related to a problem? Please describe.**
Implementing [blob detection](https://github.com/scikit-image/scikit-image/blob/main/skimage/feature/blob.py ) from [scikit-image](https://github.com/scikit-image/scikit-image) with [cupy](https://github.com/cupy/cupy) backend

**Describe the solution you'd like**
A lot of numpy/scipy functions can be one-by-one replaced with cupy functions, however, the scikit-image is using the ahead of time compiler [pythran](https://github.com/serge-sans-paille/pythran) in its background. So the `scikit-image` code has a lot of for loops which need to be replaced in either cupy arrays or using cuda kernels. An example is the function [`_blob_overlap`](https://github.com/scikit-image/scikit-image/blob/f0a3fd10edc7626baa220cd9977cdee25233ac5f/skimage/feature/blob.py#L85) which is using the function `spatial.cKDTree` which is currently not implemented in `cupyx.scipy`. Also, `_blob_overlap` uses large for loops which should be avoided. (AFAIK it is not possible to cythonized/pythranized, but `numba` would be a possibility?)

I would like to write `ElementwiseKernel` instead of `RawKernel` if that's ok. Because they are easier to handle. Those kernel are going to be embedded in the *.py file.

**Describe alternatives you've considered**
Alternative to Elementwise Kernel would be RawKernel and using *.cu files like proposed in https://github.com/rapidsai/cucim/issues/375. However, this needs more time and understanding of the GPU.

**Additional context**
If someone would review the code for its fastness I would appreciate it. Currently, we just avoid for-loops as speed improvements.



Discussion / hints:

Hi @monzelr, it would be great if you want to work on this. 

Yes, scikit-image 0.19 did start using Pythran, however there is still much more Cython-based code than Pythran overall. Either way, this sort of compiled code is less trivially converted to CuPy-based solutions and typically requires use of `ElementwiseKernel` or `RawModule` as you mention.

There are currently a couple uses of `scipy.spatial.cKDTree` in cuCIM, but this is unfortunately done on the host at the moment due to the lack of a GPU equivalent in CuPy. (It is used in the private `ensure_spacing` helper used within `feature.peak_local_max` as well as within `feature.corner_peaks`).


I don't think Pythran is being used on the specific blob detection functions you asked about, though? 

I also don't see any large for loops in `_blob_overlap`, but think perhaps you mean in `_prune_blobs` which does have such a loop over points returned by `cKDTree`.

One other tip: A few of the helper functions have many simple math operations (e.g. `_compute_sphere_overlap`). In that case, adding a simple `@cupy.fuse()` decorator to the function can fuse all of these into a single CUDA kernel and should make it several times more efficient than if each individual multiplication, division, etc. was a separate CUDA kernel. Basically `fuse` can only be used with a subset of functions, but all basic math operations and NumPy ufuncs should be supported.

We don't currently have Numba as a dependency and it is somewhat heavy-weight as it also pulls in LLVM, so l would prefer to try implementing using CuPy kernels.
I would initially try converting all components to the GPU except for `_prune_blobs`, as the other parts look reasonably straightforward. That should still give a performance benefit over scikit-image even if we have to fall back to the CPU for pruning for now. 
Hi @grlee77,

yes, you are right: `_prune_blobs` is the function which needs a lot of time on our Jetson Xavier NX device with `scikit-image` version 1.1 (Python 3.6). However, this old code used `itertools.combinations` and not `cKDTree` in [`_prune_blobs`](https://github.com/scikit-image/scikit-image/blob/ceb448b586e3a15f921c417095ae307c5d0eef3f/skimage/feature/blob.py#L73). The Elementwise Kernel for `_prune_blobs` which was derived from `scikit-image` v1.1.x is working fine and now it is pretty fast on the Jetson. It is a kind of wasted time - but at least we learned more about `ElementwiseKernel`.

However, we are going to use the latest scikit-image version 1.9 now. 1.9 and 1.1 differs a lot and basically we are rewriting everything.
Thank you for your input! I will have a look on cupy's `fuse`. 



For those who are interested; that is the new cupy code `_prune_blobs` derived from `scikit-image` v1.1:
```python
import cupy as cp
import itertools
from math import sqrt, hypot, acos
import math

prune_blobs_kernel = cp.ElementwiseKernel(
    "raw T blob, int64 n_array, float32 overlap",
    "raw T blob_out",
    r"""
    
    // *************************************************************************
    // This function is derived from Scikit-Image _blob_overlap (v0.11.x):
    // https://github.com/scikit-image/scikit-image/blob/dc10d457d1c36f18651f72b8a21dbeb55fcd34aa/skimage/feature/blob.py#L73
    // *************************************************************************
    
    float d = 0.0;
    float r1 = 0.0;
    float r2 = 0.0;
    float a = 0.0;
    float b = 0.0;
    float c = 0.0;
    float ratio1 = 0.0;
    float ratio2 = 0.0;
    float value = 0.0;
    const float pi = 3.141592653589793;
    const float root2 = 1.4142135623730951;

    //Debug purpose
    //printf("array at index %7d: %f %f %f\n", i, blob[i * 3 + 0],  blob[i * 3 + 1],  blob[i * 3 + 2]);

    for(int k=0; k<(int)(n_array/3); k++)
    {
        // blob[i] --> blob1
        // blob[k] --> blob2
        
        value = 0.0;

        if (i >= k)
        {
            // skip calculation, it is done already
            continue;
        }

        // extent of the blob is given by sqrt(2)*scale
        r1 = blob[3*i + 2] * root2;
        r2 = blob[3*k + 2] * root2;
        d = hypot(blob[3*i + 0] - blob[3*k + 0], blob[3*i + 1] - blob[3*k + 1]);
        
        if (d > r1 + r2)
        {
            value = 0.0;
        }
        else if (d <= abs(r1 - r2))
        {
            // one blob is inside the other, the smaller blob must die
            value = 1.0;
        }
        else
        {
            ratio1 = (d*d + r1*r1 - r2*r2) / (2 * d * r1);
            if (ratio1 > 1)
            {
                ratio1 = 1;
            }
            else if (ratio1 < -1)
            {
                ratio1 = -1;
            }

            ratio2 = (d*d + r2*r2 - r1*r1) / (2 * d * r2);
            if (ratio2 > 1)
            {
                ratio2 = 1;
            }
            else if (ratio2 < -1)
            {
                ratio2 = -1;
            }

            a = -d + r2 + r1;
            b = d - r2 + r1;
            c = d + r2 - r1;
            d = d + r2 + r1;

            float r = min(r1, r2);
            value = (r1*r1 * acos(ratio1) * r2*r2 * acos(ratio2) - 0.5 * sqrt(abs(a * b * c * d))) / (pi * r*r);
            
        }
        
        if (value > overlap)
        {                
            if (blob[3*i + 2] <= blob[k*3 + 2])
            {
                blob_out[3*i + 2] = -1.0;
            }
            else
            {
                blob_out[3*k + 2] = -1.0;
            }
        }
    }
    """,
    "prune_blobs_kernel")


def prune_blobs_cu(blobs_array: cp.ndarray, overlap: float) -> cp.ndarray:
    """Eliminated blobs with area overlap.

    Parameters
    ----------
    blobs_array : ndarray
        A 2d array with each row representing 3 (or 4) values,
        ``(row, col, sigma)`` or ``(pln, row, col, sigma)`` in 3D,
        where ``(row, col)`` (``(pln, row, col)``) are coordinates of the blob
        and ``sigma`` is the standard deviation of the Gaussian kernel which
        detected the blob.
        This array must not have a dimension of size 0.
    overlap : float
        A value between 0 and 1. If the fraction of area overlapping for 2
        blobs is greater than `overlap` the smaller blob is eliminated.
    Returns
    -------
    A : ndarray
        `array` with overlapping blobs removed.
    """

    prune_blobs_kernel(blobs_array, blobs_array.size, overlap, blobs_array, size=blobs_array.shape[0])
    _return = blobs_array[blobs_array[:, 2] > 0, :]
    return _return


def prune_blobs_cu_old(blobs_array, overlap):
    """Eliminated blobs with area overlap"""

    # more iterations might eliminate more blobs, but one suffices
    # for most cases
    for blob1, blob2 in itertools.combinations(blobs_array, 2):
        value = _prune_blobs_blob_overlap_old(blob1, blob2)
        if value > overlap:
            if blob1[2] > blob2[2]:
                blob2[2] = -1
            else:
                blob1[2] = -1

    # return blobs_array[blobs_array[:,2] > 0]
    return cp.array([b for b in blobs_array if b[2] > 0])


def _prune_blobs_blob_overlap_old(blob1, blob2):
    """finds the overlapping area fraction between two blobs"""

    root2 = sqrt(2)

    # extent of the blob is given by sqrt(2)*scale
    r1 = blob1[2] * root2
    r2 = blob2[2] * root2

    d = hypot(blob1[0] - blob2[0], blob1[1] - blob2[1])
    # print("d=%1.1f" % d)
    if d > r1 + r2:
        return 0

    # one blob is inside the other, the smaller blob must die
    if d <= abs(r1 - r2):
        return 1

    ratio1 = (d ** 2 + r1 ** 2 - r2 ** 2) / (2 * d * r1)
    ratio1 = cp.clip(ratio1, -1, 1)
    acos1 = cp.arccos(ratio1)

    ratio2 = (d ** 2 + r2 ** 2 - r1 ** 2) / (2 * d * r2)
    ratio2 = cp.clip(ratio2, -1, 1)
    acos2 = cp.arccos(ratio2)

    a = -d + r2 + r1
    b = d - r2 + r1
    c = d + r2 - r1
    d = d + r2 + r1
    area = r1 ** 2 * acos1 * r2 ** 2 * acos2 - 0.5 * sqrt(abs(a * b * c * d))

    return area / (math.pi * (min(r1, r2) ** 2))


if __name__ == "__main__":
    # just a small example how to use

    array = cp.array([[401., 123., 10.666667],
                      [402., 126., 1.],
                      [416., 123., 7.4444447],
                      [499., 386., 4.2222223],
                      [166., 254., 4.2222223],
                      [218., 146., 4.2222223]], dtype=cp.float32)

    print("\n")
    old = prune_blobs_cu_old(array.copy(), 0.5)
    new = prune_blobs_cu(array.copy(), 0.5)

    print("input:\n", array)
    print("old:\n", old)
    print("new:\n", new)
    print("True? -> ", cp.array_equal(old, new))

```
Oh, great, I didn't realize you already had something working. Please make a PR when you are ready and we can help review.

 > scikit-image version 1.1 (Python 3.6).

1.1 sounds like it may be a scikit-learn version? The latest is 0.19.3 for scikit-image :)

Yes, I meant `scikit-image` [v0.11](https://github.com/scikit-image/scikit-image/tree/v0.11.x) x) 

Fix patch:
diff --git a/benchmarks/skimage/cucim_feature_bench.py b/benchmarks/skimage/cucim_feature_bench.py
index 05f61493f..ac55f9f7d 100644
--- a/benchmarks/skimage/cucim_feature_bench.py
+++ b/benchmarks/skimage/cucim_feature_bench.py
@@ -1,18 +1,52 @@
+import math
 import os
 import pickle
 import argparse
 
-import cucim.skimage
-import cucim.skimage.feature
 import cupy as cp
 import numpy as np
 import pandas as pd
 import skimage
 import skimage.feature
+from skimage import data, draw
+
+import cucim.skimage
+import cucim.skimage.feature
+from cucim.skimage import exposure
 
 from _image_bench import ImageBench
 
 
+class BlobDetectionBench(ImageBench):
+    def set_args(self, dtype):
+        ndim = len(self.shape)
+        if ndim == 2:
+            # create 2D image by tiling the coins image
+            img = cp.array(data.coins())
+            img = exposure.equalize_hist(img)  # improves detection
+            n_tile = (math.ceil(s/s0) for s, s0 in zip(self.shape, img.shape))
+            img = cp.tile(img, n_tile)
+            img = img[tuple(slice(s) for s in img.shape)]
+            img = img.astype(dtype, copy=False)
+        elif ndim == 3:
+            # create 3D volume with randomly positioned ellipses
+            e = cp.array(draw.ellipsoid(*(max(s // 20, 1) for s in self.shape)))
+            e = e.astype(dtype, copy=False)
+            img = cp.zeros(self.shape, dtype=dtype)
+            rng = np.random.default_rng(5)
+            num_ellipse = 64
+            offsets = rng.integers(0, np.prod(img.shape), num_ellipse)
+            locs = np.unravel_index(offsets, img.shape)
+            for loc in zip(*locs):
+                loc = tuple(min(l, s - es) for l, s, es in zip(loc, img.shape, e.shape))
+                sl = tuple(slice(l, l + es) for l, es in zip(loc, e.shape))
+                img[sl] = e
+        else:
+            raise NotImplementedError("unsupported ndim")
+        self.args_gpu = (img,)
+        self.args_cpu = (cp.asnumpy(img),)
+
+
 class MatchTemplateBench(ImageBench):
     def set_args(self, dtype):
         rstate = cp.random.RandomState(5)
@@ -27,6 +61,7 @@ def set_args(self, dtype):
         self.args_cpu = (image, template)
         self.args_gpu = (imaged, templated)
 
+
 def main(args):
 
     pfile = "cucim_feature_results.pickle"
@@ -58,11 +93,15 @@ def main(args):
         ("corner_shi_tomasi", dict(sigma=1), dict(), False, False),
         ("corner_foerstner", dict(sigma=1), dict(), False, False),
         ("corner_peaks", dict(), dict(min_distance=(2, 3, 5)), False, True),
-        ("match_template", dict(), dict(pad_input=[False], mode=["reflect"]), False, True)
+        ("match_template", dict(), dict(pad_input=[False], mode=["reflect"]), False, True),
+        # blob detectors. fixed kwargs here are taken from the docstring examples
+        ("blob_dog", dict(threshold=.05, min_sigma=10, max_sigma=40), dict(), False, True),
+        ("blob_log", dict(threshold=.3), dict(), False, True),
+        ("blob_doh", dict(), dict(), False, False),
     ]:
 
         if function_name == args.func_name:
-            shape = tuple(list(map(int,(args.img_size.split(',')))))
+            shape = tuple(list(map(int,(args.img_size.split(",")))))
         else:
             continue
 
@@ -70,18 +109,36 @@ def main(args):
             # skip any large sizes that take too long
         ndim = len(shape)
         run_cpu = not args.no_cpu
+        if not allow_nd:
+            if not allow_color:
+                if ndim > 2:
+                    continue
+            else:
+                if ndim > 3 or (ndim == 3 and shape[-1] not in [3, 4]):
+                    continue
+
+        if shape[-1] == 3 and not allow_color:
+            continue
 
-        if function_name != "match_template":
-            if not allow_nd:
-                if not allow_color:
-                    if ndim > 2:
-                        continue
+        if function_name.startswith("blob"):
+            if ndim == 3:
+                # set more reasonable threshold for the synthetic 3d data
+                if function_name == "blob_log":
+                    fixed_kwargs = {"threshold": 0.6}
                 else:
-                    if ndim > 3 or (ndim == 3 and shape[-1] not in [3, 4]):
-                        continue
+                    fixed_kwargs = {}
 
-            if shape[-1] == 3 and not allow_color:
-                continue
+            B = BlobDetectionBench(
+                function_name=function_name,
+                shape=shape,
+                dtypes=dtypes,
+                fixed_kwargs=fixed_kwargs,
+                var_kwargs=var_kwargs,
+                module_cpu=skimage.feature,
+                module_gpu=cucim.skimage.feature,
+                run_cpu=run_cpu,
+            )
+        elif function_name != "match_template":
 
             if function_name == "multiscale_basic_features":
                 fixed_kwargs["channel_axis"] = -1 if shape[-1] == 3 else None
@@ -100,15 +157,6 @@ def main(args):
                 run_cpu=run_cpu,
             )
         else:
-            if not allow_nd:
-                if allow_color:
-                    if ndim > 2:
-                        continue
-                else:
-                    if ndim > 3 or (ndim == 3 and shape[-1] not in [3, 4]):
-                        continue
-            if shape[-1] == 3 and not allow_color:
-                continue
 
             B = MatchTemplateBench(
                 function_name=function_name,
@@ -131,15 +179,40 @@ def main(args):
         f.write(all_results.to_markdown())
 
 
-if __name__ == '__main__':
-    parser = argparse.ArgumentParser(description='Benchmarking cuCIM Feature')
-    func_name_choices = ["multiscale_basic_features","canny","daisy","structure_tensor","hessian_matrix","hessian_matrix_det","shape_index","corner_kitchen_rosenfeld","corner_harris","corner_shi_tomasi","corner_foerstner","corner_peaks","match_template"]
-    dtype_choices = ['float16', 'float32', 'float64', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64']
-    parser.add_argument('-i','--img_size', type=str, help='Size of input image', required=True)
-    parser.add_argument('-d','--dtype', type=str, help='Dtype of input image', choices=dtype_choices, required=True)
-    parser.add_argument('-f','--func_name', type=str, help='function to benchmark', choices=func_name_choices, required=True)
-    parser.add_argument('-t','--duration', type=int, help='time to run benchmark', required=True)
-    parser.add_argument('--no_cpu', action='store_true', help='disable cpu measurements', default=False)
-
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Benchmarking cuCIM Feature")
+    # fmt: off
+    func_name_choices = [
+        "multiscale_basic_features", "canny", "daisy", "structure_tensor",
+        "hessian_matrix", "hessian_matrix_det", "shape_index",
+        "corner_kitchen_rosenfeld", "corner_harris", "corner_shi_tomasi",
+        "corner_foerstner", "corner_peaks", "match_template", "blob_dog",
+        "blob_log", "blob_doh",
+    ]
+    dtype_choices = [
+        "float16", "float32", "float64", "int8", "int16", "int32", "int64",
+        "uint8", "uint16", "uint32", "uint64"
+    ]
+    parser.add_argument(
+        "-i", "--img_size", type=str, help="Size of input image",
+        required=True
+    )
+    parser.add_argument(
+        "-d", "--dtype", type=str, help="Dtype of input image",
+        choices=dtype_choices, required=True
+    )
+    parser.add_argument(
+        "-f", "--func_name", type=str, help="function to benchmark",
+        choices=func_name_choices, required=True
+    )
+    parser.add_argument(
+        "-t", "--duration", type=int, help="time to run benchmark",
+        required=True
+    )
+    parser.add_argument(
+        "--no_cpu", action="store_true", help="disable cpu measurements",
+        default=False
+    )
+    # fmt: on
     args = parser.parse_args()
     main(args)
diff --git a/benchmarks/skimage/run-nv-bench-feature.sh b/benchmarks/skimage/run-nv-bench-feature.sh
index fbab913c6..bc572da2e 100755
--- a/benchmarks/skimage/run-nv-bench-feature.sh
+++ b/benchmarks/skimage/run-nv-bench-feature.sh
@@ -1,11 +1,11 @@
 #!/bin/bash
 param_shape=(512,512 3840,2160 3840,2160,3 192,192,192)
-param_filt=(multiscale_basic_features canny daisy structure_tensor hessian_matrix hessian_matrix_det shape_index corner_kitchen_rosenfeld corner_harris corner_shi_tomasi corner_foerstner corner_peaks match_template)
-param_dt=(float64 float32 float16)
+param_filt=(multiscale_basic_features canny daisy structure_tensor hessian_matrix hessian_matrix_det shape_index corner_kitchen_rosenfeld corner_harris corner_shi_tomasi corner_foerstner corner_peaks match_template blob_dog blob_log blob_doh)
+param_dt=(float64 float32)
 for shape in "${param_shape[@]}"; do
     for filt in "${param_filt[@]}"; do
         for dt in "${param_dt[@]}"; do
-            python cucim_feature_bench.py -f $filt -i $shape -d $dt -t 10
+            python cucim_feature_bench.py -f $filt -i $shape -d $dt -t 3
         done
     done
 done
diff --git a/python/cucim/src/cucim/skimage/feature/__init__.py b/python/cucim/src/cucim/skimage/feature/__init__.py
index 301816382..9f535a694 100644
--- a/python/cucim/src/cucim/skimage/feature/__init__.py
+++ b/python/cucim/src/cucim/skimage/feature/__init__.py
@@ -8,6 +8,7 @@
                      structure_tensor, structure_tensor_eigenvalues)
 from .match import match_descriptors
 from .peak import peak_local_max
+from .blob import blob_dog, blob_log, blob_doh
 from .template import match_template
 
 
@@ -32,4 +33,7 @@
            # 'corner_fast',
            # 'corner_orientations',
            'match_template',
-           'match_descriptors']
+           'match_descriptors',
+           'blob_dog',
+           'blob_log',
+           'blob_doh']
diff --git a/python/cucim/src/cucim/skimage/feature/_hessian_det_appx.py b/python/cucim/src/cucim/skimage/feature/_hessian_det_appx.py
new file mode 100644
index 000000000..28f893c9a
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/feature/_hessian_det_appx.py
@@ -0,0 +1,113 @@
+import math
+import os
+
+import cupy as cp
+
+
+def _dtype_to_cuda_float_type(dtype):
+    """Maps a float data type from cupy to cuda.
+
+    Returns a cuda (c++) data type.
+
+    Parameters
+    ----------
+    dtype : cupy dtype
+        A cupy dtype from type float.
+
+    Returns
+    -------
+    cpp_float_type : str
+        Supported cuda data type
+    """
+    cpp_float_types = {
+        cp.float32: 'float',
+        cp.float64: 'double',
+    }
+    dtype = cp.dtype(dtype)
+    if dtype.type not in cpp_float_types:
+        raise ValueError(f"unrecognized dtype: {dtype.type}")
+    return cpp_float_types[dtype.type]
+
+
+@cp.memoize()
+def _get_hessian_det_appx_kernel(dtype, large_int) -> cp.RawModule:
+    """Loads all kernels in cuda/_hessian_det_appx.cu.
+    Returns a cupy RawModule.
+
+    Parameters
+    ----------
+    dtype : cupy dtype
+        Only the cupy dtypes float32 and float64 are supported.
+
+    Returns
+    -------
+    out : RawModule
+        A cupy RawModule containing the __global__ functions
+        `_hessian_matrix_det`.
+    """
+    image_t = _dtype_to_cuda_float_type(dtype)
+
+    int_t = 'long long' if large_int else 'int'
+
+    _preamble = f"""
+#define IMAGE_T {image_t}
+#define INT_T {int_t}
+        """
+
+    kernel_directory = os.path.join(
+        os.path.normpath(os.path.dirname(__file__)), 'cuda')
+    cu_file = os.path.join(kernel_directory, "_hessian_det_appx.cu")
+    with open(cu_file, 'rt') as f:
+        _code = f.read()
+
+    return cp.RawModule(code=_preamble + _code,
+                        options=('--std=c++11',),
+                        name_expressions=["_hessian_matrix_det"])
+
+
+def _hessian_matrix_det(img: cp.ndarray, sigma) -> cp.ndarray:
+    """Compute the approximate Hessian Determinant over a 2D image.
+
+    This method uses box filters over integral images to compute the
+    approximate Hessian Determinant as described in [1]_.
+
+    Parameters
+    ----------
+    img : array
+        The integral image over which to compute Hessian Determinant.
+    sigma : float
+        Standard deviation used for the Gaussian kernel, used for the Hessian
+        matrix
+
+    Returns
+    -------
+    out : array
+        The array of the Determinant of Hessians.
+
+    References
+    ----------
+    .. [1] Herbert Bay, Andreas Ess, Tinne Tuytelaars, Luc Van Gool,
+           "SURF: Speeded Up Robust Features"
+           ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf
+
+    Notes
+    -----
+    The running time of this method only depends on size of the image. It is
+    independent of `sigma` as one would expect. The downside is that the
+    result for `sigma` less than `3` is not accurate, i.e., not similar to
+    the result obtained if someone computed the Hessian and took its
+    determinant.
+    """
+    rawmodule = _get_hessian_det_appx_kernel(img.dtype, max(img.shape) > 2**31)
+    _hessian_det_appx_kernel = rawmodule.get_function("_hessian_matrix_det")
+
+    out = cp.empty_like(img, dtype=img.dtype)
+
+    block_size = 64
+    grid_size = int(math.ceil(img.size / block_size))
+    _hessian_det_appx_kernel(
+        (grid_size,),
+        (block_size,),
+        (img.ravel(), img.shape[0], img.shape[1], float(sigma), out)
+    )
+    return out
diff --git a/python/cucim/src/cucim/skimage/feature/blob.py b/python/cucim/src/cucim/skimage/feature/blob.py
new file mode 100644
index 000000000..c2b121a03
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/feature/blob.py
@@ -0,0 +1,626 @@
+import math
+import os
+
+import cupy as cp
+import numpy as np
+from cupyx.scipy.ndimage import gaussian_laplace
+
+from .._shared.filters import gaussian
+from .._shared.utils import _supported_float_type, check_nD
+from ..transform import integral_image
+from ..util import img_as_float
+from ._hessian_det_appx import _hessian_matrix_det
+from .peak import peak_local_max
+
+# This basic blob detection algorithm is based on:
+# http://www.cs.utah.edu/~jfishbau/advimproc/project1/ (04.04.2013)
+# Theory behind: https://en.wikipedia.org/wiki/Blob_detection (04.04.2013)
+
+
+def _dtype_to_cuda_float_type(dtype):
+    """Maps a float data type from cupy to cuda.
+
+    Returns a cuda (c++) data type.
+
+    Parameters
+    ----------
+    dtype : cupy dtype
+        A cupy dtype from type float.
+
+    Returns
+    -------
+    cpp_float_type : cuda (c++) data type
+        Supported cuda data type
+    """
+    cpp_float_types = {
+        cp.float32: 'float',
+        cp.float64: 'double',
+    }
+    dtype = cp.dtype(dtype)
+    if dtype.type not in cpp_float_types:
+        raise ValueError(f"unrecognized dtype: {dtype.type}")
+    return cpp_float_types[dtype.type]
+
+
+@cp.memoize()
+def _get_prune_blob_rawmodule(dtype, large_int) -> cp.RawModule:
+    """Loads the kernel according to dtype /cuda/blob.cu
+    Returns a cupy RawModule.
+
+    Parameters
+    ----------
+    dtype : cupy dtype
+        Only the cupy dtypes float32 and float64 are supported.
+
+    Returns
+    -------
+    RawModule : cupy RawModule
+        A cupy RawModule containing the __global__ functions `_prune_blobs`.
+    """
+    blob_t = _dtype_to_cuda_float_type(dtype)
+    int_t = 'long long' if large_int else 'int'
+
+    _preamble = f"""
+#define BLOB_T {blob_t}
+#define INT_T {int_t}
+    """
+
+    kernel_directory = os.path.join(
+        os.path.normpath(os.path.dirname(__file__)), 'cuda'
+    )
+    with open(os.path.join(kernel_directory, "blob.cu"), 'rt') as f:
+        _code = f.read()
+
+    return cp.RawModule(
+        code=_preamble + _code,
+        options=('--std=c++11',),
+        name_expressions=["_prune_blobs"]
+    )
+
+
+def _prune_blobs(blobs_array, overlap, *, sigma_dim=1):
+    """Eliminated blobs with area overlap.
+
+    Parameters
+    ----------
+    blobs_array : ndarray
+        A 2d array with each row representing 3 (or 4) values,
+        ``(row, col, sigma)`` or ``(pln, row, col, sigma)`` in 3D,
+        where ``(row, col)`` (``(pln, row, col)``) are coordinates of the blob
+        and ``sigma`` is the standard deviation of the Gaussian kernel which
+        detected the blob.
+        This array must not have a dimension of size 0.
+    overlap : float
+        A value between 0 and 1. If the fraction of area overlapping for 2
+        blobs is greater than `overlap` the smaller blob is eliminated.
+    sigma_dim : int, optional
+        The number of columns in ``blobs_array`` corresponding to sigmas rather
+        than positions.
+
+    Returns
+    -------
+    A : ndarray
+        `blobs_array` with overlapping blobs removed.
+    """
+
+    # from here, the kernel does the calculation
+    blobs_module = _get_prune_blob_rawmodule(blobs_array.dtype,
+                                             max(blobs_array.shape) > 2**31)
+    _prune_blobs_kernel = blobs_module.get_function("_prune_blobs")
+
+    block_size = 64
+    grid_size = int(math.ceil(blobs_array.shape[0] / block_size))
+    _prune_blobs_kernel((grid_size,), (block_size,),
+                        (blobs_array.ravel(),
+                         int(blobs_array.shape[0]),
+                         int(blobs_array.shape[1]),
+                         float(overlap),
+                         int(sigma_dim))
+                        )
+    return blobs_array[blobs_array[:, -1] > 0, :]
+
+
+def _format_exclude_border(img_ndim, exclude_border):
+    """Format an ``exclude_border`` argument as a tuple of ints for calling
+    ``peak_local_max``.
+    """
+    if isinstance(exclude_border, tuple):
+        if len(exclude_border) != img_ndim:
+            raise ValueError(
+                "`exclude_border` should have the same length as the "
+                "dimensionality of the image.")
+        for exclude in exclude_border:
+            if not isinstance(exclude, int):
+                raise ValueError(
+                    "exclude border, when expressed as a tuple, must only "
+                    "contain ints.")
+        return exclude_border
+    elif isinstance(exclude_border, int):
+        return (exclude_border,) * img_ndim + (0,)
+    elif exclude_border is True:
+        raise ValueError("exclude_border cannot be True")
+    elif exclude_border is False:
+        return (0,) * (img_ndim + 1)
+    else:
+        raise ValueError(
+            f'Unsupported value ({exclude_border}) for exclude_border'
+        )
+
+
+def _prep_sigmas(ndim, min_sigma, max_sigma):
+    # if both min and max sigma are scalar, function returns only one sigma
+    scalar_max = np.isscalar(max_sigma)
+    scalar_min = np.isscalar(min_sigma)
+    scalar_sigma = scalar_max and scalar_min
+
+    # Gaussian filter requires that sequence-type sigmas have same
+    # dimensionality as image. This broadcasts scalar kernels
+    if scalar_max:
+        max_sigma = (max_sigma,) * ndim
+    if scalar_min:
+        min_sigma = (min_sigma,) * ndim
+    return scalar_sigma, min_sigma, max_sigma
+
+
+def blob_dog(image, min_sigma=1, max_sigma=50, sigma_ratio=1.6, threshold=0.5,
+             overlap=.5, *, threshold_rel=None, exclude_border=False):
+    r"""Finds blobs in the given grayscale image.
+    Blobs are found using the Difference of Gaussian (DoG) method [1]_, [2]_.
+    For each blob found, the method returns its coordinates and the standard
+    deviation of the Gaussian kernel that detected the blob.
+    Parameters
+    ----------
+    image : ndarray
+        Input grayscale image, blobs are assumed to be light on dark
+        background (white on black).
+    min_sigma : scalar or sequence of scalars, optional
+        The minimum standard deviation for Gaussian kernel. Keep this low to
+        detect smaller blobs. The standard deviations of the Gaussian filter
+        are given for each axis as a sequence, or as a single number, in
+        which case it is equal for all axes.
+    max_sigma : scalar or sequence of scalars, optional
+        The maximum standard deviation for Gaussian kernel. Keep this high to
+        detect larger blobs. The standard deviations of the Gaussian filter
+        are given for each axis as a sequence, or as a single number, in
+        which case it is equal for all axes.
+    sigma_ratio : float, optional
+        The ratio between the standard deviation of Gaussian Kernels used for
+        computing the Difference of Gaussians
+    threshold : float or None, optional
+        The absolute lower bound for scale space maxima. Local maxima smaller
+        than `threshold` are ignored. Reduce this to detect blobs with lower
+        intensities. If `threshold_rel` is also specified, whichever threshold
+        is larger will be used. If None, `threshold_rel` is used instead.
+    overlap : float, optional
+        A value between 0 and 1. If the area of two blobs overlaps by a
+        fraction greater than `threshold`, the smaller blob is eliminated.
+    threshold_rel : float or None, optional
+        Minimum intensity of peaks, calculated as
+        ``max(dog_space) * threshold_rel``, where ``dog_space`` refers to the
+        stack of Difference-of-Gaussian (DoG) images computed internally. This
+        should have a value between 0 and 1. If None, `threshold` is used
+        instead.
+    exclude_border : tuple of ints, int, or False, optional
+        If tuple of ints, the length of the tuple must match the input array's
+        dimensionality.  Each element of the tuple will exclude peaks from
+        within `exclude_border`-pixels of the border of the image along that
+        dimension.
+        If nonzero int, `exclude_border` excludes peaks from within
+        `exclude_border`-pixels of the border of the image.
+        If zero or False, peaks are identified regardless of their
+        distance from the border.
+
+    Returns
+    -------
+    A : (n, image.ndim + sigma) ndarray
+        A 2d array with each row representing 2 coordinate values for a 2D
+        image, or 3 coordinate values for a 3D image, plus the sigma(s) used.
+        When a single sigma is passed, outputs are:
+        ``(r, c, sigma)`` or ``(p, r, c, sigma)`` where ``(r, c)`` or
+        ``(p, r, c)`` are coordinates of the blob and ``sigma`` is the standard
+        deviation of the Gaussian kernel which detected the blob. When an
+        anisotropic gaussian is used (sigmas per dimension), the detected sigma
+        is returned for each dimension.
+
+    See also
+    --------
+    cucim.skimage.filters.difference_of_gaussians
+
+    References
+    ----------
+    .. [1] https://en.wikipedia.org/wiki/Blob_detection#The_difference_of_Gaussians_approach  # noqa
+    .. [2] Lowe, D. G. "Distinctive Image Features from Scale-Invariant
+        Keypoints." International Journal of Computer Vision 60, 91–110 (2004).
+        https://www.cs.ubc.ca/~lowe/papers/ijcv04.pdf
+        :DOI:`10.1023/B:VISI.0000029664.99615.94`
+
+    Examples
+    --------
+    >>> import cupy as cp
+    >>> from skimage import data
+    >>> from cucim.skimage import feature
+    >>> coins = cp.array(data.coins())
+    >>> feature.blob_dog(coins, threshold=.05, min_sigma=10, max_sigma=40)
+    array([[128., 155.,  10.],
+           [198., 155.,  10.],
+           [124., 338.,  10.],
+           [127., 102.,  10.],
+           [193., 281.,  10.],
+           [126., 208.,  10.],
+           [267., 115.,  10.],
+           [197., 102.,  10.],
+           [198., 215.,  10.],
+           [123., 279.,  10.],
+           [126.,  46.,  10.],
+           [259., 247.,  10.],
+           [196.,  43.,  10.],
+           [ 54., 276.,  10.],
+           [267., 358.,  10.],
+           [ 58., 100.,  10.],
+           [259., 305.,  10.],
+           [185., 347.,  16.],
+           [261., 174.,  16.],
+           [ 46., 336.,  16.],
+           [ 54., 217.,  10.],
+           [ 55., 157.,  10.],
+           [ 57.,  41.,  10.],
+           [260.,  47.,  16.]])
+
+    Notes
+    -----
+    The radius of each blob is approximately :math:`\sqrt{2}\sigma` for
+    a 2-D image and :math:`\sqrt{3}\sigma` for a 3-D image.
+    """
+    image = img_as_float(image)
+    float_dtype = _supported_float_type(image.dtype)
+    image = image.astype(float_dtype, copy=False)
+
+    # Determine if provixed sigmas are scalar and broadcast to image.ndim
+    scalar_sigma, min_sigma, max_sigma = _prep_sigmas(image.ndim, min_sigma,
+                                                      max_sigma)
+
+    if sigma_ratio <= 1.0:
+        raise ValueError('sigma_ratio must be > 1.0')
+
+    # k such that min_sigma*(sigma_ratio**k) > max_sigma
+    log_ratio = math.log(sigma_ratio)
+    k = sum(math.log(max_s / min_s) / log_ratio + 1
+            for max_s, min_s in zip(max_sigma, min_sigma))
+    k /= len(min_sigma)
+    k = int(k)
+
+    # a geometric progression of standard deviations for gaussian kernels
+    ratio_powers = tuple(sigma_ratio ** i for i in range(k + 1))
+    sigma_list = tuple(tuple(s * p for s in min_sigma)
+                       for p in ratio_powers)
+
+    # computing difference between two successive Gaussian blurred images
+    # to obtain an approximation of the scale invariant Laplacian of the
+    # Gaussian operator
+    dog_image_cube = cp.empty(image.shape + (k,), dtype=float_dtype)
+    gaussian_previous = gaussian(image, sigma_list[0], mode='reflect')
+    for i, s in enumerate(sigma_list[1:]):
+        gaussian_current = gaussian(image, s, mode='reflect')
+        dog_image_cube[..., i] = gaussian_previous - gaussian_current
+        gaussian_previous = gaussian_current
+
+    # normalization factor for consistency in DoG magnitude
+    sf = 1 / (sigma_ratio - 1)
+    dog_image_cube *= sf
+
+    exclude_border = _format_exclude_border(image.ndim, exclude_border)
+    local_maxima = peak_local_max(
+        dog_image_cube,
+        threshold_abs=threshold,
+        threshold_rel=threshold_rel,
+        exclude_border=exclude_border,
+        footprint=cp.ones((3,) * (image.ndim + 1)),
+    )
+
+    # Catch no peaks
+    if local_maxima.size == 0:
+        return cp.empty((0, image.ndim + (1 if scalar_sigma else image.ndim)))
+
+    # Convert local_maxima to float64
+    lm = local_maxima.astype(float_dtype)
+
+    # translate final column of lm, which contains the index of the
+    # sigma that produced the maximum intensity value, into the sigma
+    sigma_list = cp.asarray(sigma_list)
+    sigmas_of_peaks = sigma_list[local_maxima[:, -1]]
+
+    if scalar_sigma:
+        # select one sigma column, keeping dimension
+        sigmas_of_peaks = sigmas_of_peaks[:, 0:1]
+
+    # Remove sigma index and replace with sigmas
+    lm = cp.hstack([lm[:, :-1], sigmas_of_peaks])
+
+    sigma_dim = sigmas_of_peaks.shape[1]
+
+    return _prune_blobs(lm, overlap, sigma_dim=sigma_dim)
+
+
+def blob_log(image, min_sigma=1, max_sigma=50, num_sigma=10, threshold=.2,
+             overlap=.5, log_scale=False, *, threshold_rel=None,
+             exclude_border=False):
+    r"""Finds blobs in the given grayscale image.
+    Blobs are found using the Laplacian of Gaussian (LoG) method [1]_.
+    For each blob found, the method returns its coordinates and the standard
+    deviation of the Gaussian kernel that detected the blob.
+    Parameters
+    ----------
+    image : ndarray
+        Input grayscale image, blobs are assumed to be light on dark
+        background (white on black).
+    min_sigma : scalar or sequence of scalars, optional
+        the minimum standard deviation for Gaussian kernel. Keep this low to
+        detect smaller blobs. The standard deviations of the Gaussian filter
+        are given for each axis as a sequence, or as a single number, in
+        which case it is equal for all axes.
+    max_sigma : scalar or sequence of scalars, optional
+        The maximum standard deviation for Gaussian kernel. Keep this high to
+        detect larger blobs. The standard deviations of the Gaussian filter
+        are given for each axis as a sequence, or as a single number, in
+        which case it is equal for all axes.
+    num_sigma : int, optional
+        The number of intermediate values of standard deviations to consider
+        between `min_sigma` and `max_sigma`.
+    threshold : float or None, optional
+        The absolute lower bound for scale space maxima. Local maxima smaller
+        than `threshold` are ignored. Reduce this to detect blobs with lower
+        intensities. If `threshold_rel` is also specified, whichever threshold
+        is larger will be used. If None, `threshold_rel` is used instead.
+    overlap : float, optional
+        A value between 0 and 1. If the area of two blobs overlaps by a
+        fraction greater than `threshold`, the smaller blob is eliminated.
+    log_scale : bool, optional
+        If set intermediate values of standard deviations are interpolated
+        using a logarithmic scale to the base `10`. If not, linear
+        interpolation is used.
+    threshold_rel : float or None, optional
+        Minimum intensity of peaks, calculated as
+        ``max(log_space) * threshold_rel``, where ``log_space`` refers to the
+        stack of Laplacian-of-Gaussian (LoG) images computed internally. This
+        should have a value between 0 and 1. If None, `threshold` is used
+        instead.
+    exclude_border : tuple of ints, int, or False, optional
+        If tuple of ints, the length of the tuple must match the input array's
+        dimensionality.  Each element of the tuple will exclude peaks from
+        within `exclude_border`-pixels of the border of the image along that
+        dimension.
+        If nonzero int, `exclude_border` excludes peaks from within
+        `exclude_border`-pixels of the border of the image.
+        If zero or False, peaks are identified regardless of their
+        distance from the border.
+
+    Returns
+    -------
+    A : (n, image.ndim + sigma) ndarray
+        A 2d array with each row representing 2 coordinate values for a 2D
+        image, or 3 coordinate values for a 3D image, plus the sigma(s) used.
+        When a single sigma is passed, outputs are:
+        ``(r, c, sigma)`` or ``(p, r, c, sigma)`` where ``(r, c)`` or
+        ``(p, r, c)`` are coordinates of the blob and ``sigma`` is the standard
+        deviation of the Gaussian kernel which detected the blob. When an
+        anisotropic gaussian is used (sigmas per dimension), the detected sigma
+        is returned for each dimension.
+
+    References
+    ----------
+    .. [1] https://en.wikipedia.org/wiki/Blob_detection#The_Laplacian_of_Gaussian  # noqa
+
+    Examples
+    --------
+    >>> import cupy as cp
+    >>> from skimage import data
+    >>> from cucim.skimage import feature, exposure
+    >>> img = cp.array(data.coins())
+    >>> img = exposure.equalize_hist(img)  # improves detection
+    >>> feature.blob_log(img, threshold = .3)
+    array([[124.        , 336.        ,  11.88888889],
+           [198.        , 155.        ,  11.88888889],
+           [194.        , 213.        ,  17.33333333],
+           [121.        , 272.        ,  17.33333333],
+           [263.        , 244.        ,  17.33333333],
+           [194.        , 276.        ,  17.33333333],
+           [266.        , 115.        ,  11.88888889],
+           [128.        , 154.        ,  11.88888889],
+           [260.        , 174.        ,  17.33333333],
+           [198.        , 103.        ,  11.88888889],
+           [126.        , 208.        ,  11.88888889],
+           [127.        , 102.        ,  11.88888889],
+           [263.        , 302.        ,  17.33333333],
+           [197.        ,  44.        ,  11.88888889],
+           [185.        , 344.        ,  17.33333333],
+           [126.        ,  46.        ,  11.88888889],
+           [113.        , 323.        ,   1.        ]])
+    Notes
+    -----
+    The radius of each blob is approximately :math:`\sqrt{2}\sigma` for
+    a 2-D image and :math:`\sqrt{3}\sigma` for a 3-D image.
+    """
+    image = img_as_float(image)
+    float_dtype = _supported_float_type(image.dtype)
+    image = image.astype(float_dtype, copy=False)
+
+    # Determine if provixed sigmas are scalar and broadcast to image.ndim
+    scalar_sigma, min_sigma, max_sigma = _prep_sigmas(image.ndim, min_sigma,
+                                                      max_sigma)
+
+    if log_scale:
+        start = tuple(math.log10(s) for s in min_sigma)
+        stop = tuple(math.log10(s) for s in max_sigma)
+        sigma_list = np.logspace(start, stop, num_sigma)
+    else:
+        sigma_list = np.linspace(min_sigma, max_sigma, num_sigma)
+
+    def _mean_sq(s):
+        # mean-squared for small tuples (avoid numpy/cupy overhead)
+        m = sum(s) / len(s)
+        return m * m
+
+    # computing gaussian laplace
+    image_cube = cp.empty(image.shape + (len(sigma_list),), dtype=float_dtype)
+    for i, s in enumerate(sigma_list):
+        # average s**2 provides scale invariance
+        image_cube[..., i] = -gaussian_laplace(image, s) * _mean_sq(s)
+
+    exclude_border = _format_exclude_border(image.ndim, exclude_border)
+    local_maxima = peak_local_max(
+        image_cube,
+        threshold_abs=threshold,
+        threshold_rel=threshold_rel,
+        exclude_border=exclude_border,
+        footprint=cp.ones((3,) * (image.ndim + 1))
+    )
+
+    # Catch no peaks
+    if local_maxima.size == 0:
+        return cp.empty((0, image.ndim + (1 if scalar_sigma else image.ndim)))
+
+    # Convert local_maxima to float64
+    lm = local_maxima.astype(float_dtype)
+
+    # translate final column of lm, which contains the index of the
+    # sigma that produced the maximum intensity value, into the sigma
+    sigma_list = cp.asarray(sigma_list)
+    sigmas_of_peaks = sigma_list[local_maxima[:, -1]]
+
+    if scalar_sigma:
+        # select one sigma column, keeping dimension
+        sigmas_of_peaks = sigmas_of_peaks[:, 0:1]
+
+    # Remove sigma index and replace with sigmas
+    lm = cp.hstack([lm[:, :-1], sigmas_of_peaks])
+
+    sigma_dim = sigmas_of_peaks.shape[1]
+
+    return _prune_blobs(lm, overlap, sigma_dim=sigma_dim)
+
+
+def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=0.01,
+             overlap=.5, log_scale=False, *, threshold_rel=None):
+    """Finds blobs in the given grayscale image.
+
+    Blobs are found using the Determinant of Hessian method [1]_. For each blob
+    found, the method returns its coordinates and the standard deviation
+    of the Gaussian Kernel used for the Hessian matrix whose determinant
+    detected the blob. Determinant of Hessians is approximated using [2]_.
+
+    Parameters
+    ----------
+    image : 2D ndarray
+        Input grayscale image.Blobs can either be light on dark or vice versa.
+    min_sigma : float, optional
+        The minimum standard deviation for Gaussian Kernel used to compute
+        Hessian matrix. Keep this low to detect smaller blobs.
+    max_sigma : float, optional
+        The maximum standard deviation for Gaussian Kernel used to compute
+        Hessian matrix. Keep this high to detect larger blobs.
+    num_sigma : int, optional
+        The number of intermediate values of standard deviations to consider
+        between `min_sigma` and `max_sigma`.
+    threshold : float or None, optional
+        The absolute lower bound for scale space maxima. Local maxima smaller
+        than `threshold` are ignored. Reduce this to detect blobs with lower
+        intensities. If `threshold_rel` is also specified, whichever threshold
+        is larger will be used. If None, `threshold_rel` is used instead.
+    overlap : float, optional
+        A value between 0 and 1. If the area of two blobs overlaps by a
+        fraction greater than `threshold`, the smaller blob is eliminated.
+    log_scale : bool, optional
+        If set intermediate values of standard deviations are interpolated
+        using a logarithmic scale to the base `10`. If not, linear
+        interpolation is used.
+    threshold_rel : float or None, optional
+        Minimum intensity of peaks, calculated as
+        ``max(doh_space) * threshold_rel``, where ``doh_space`` refers to the
+        stack of Determinant-of-Hessian (DoH) images computed internally. This
+        should have a value between 0 and 1. If None, `threshold` is used
+        instead.
+
+    Returns
+    -------
+    A : (n, 3) ndarray
+        A 2d array with each row representing 3 values, ``(y,x,sigma)``
+        where ``(y,x)`` are coordinates of the blob and ``sigma`` is the
+        standard deviation of the Gaussian kernel of the Hessian Matrix whose
+        determinant detected the blob.
+
+    References
+    ----------
+    .. [1] https://en.wikipedia.org/wiki/Blob_detection#The_determinant_of_the_Hessian  # noqa
+    .. [2] Herbert Bay, Andreas Ess, Tinne Tuytelaars, Luc Van Gool,
+           "SURF: Speeded Up Robust Features"
+           ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf
+
+    Examples
+    --------
+    >>> import cupy as cp
+    >>> from skimage import data
+    >>> from cucim.skimage import feature
+    >>> img = cp.array(data.coins())
+    >>> feature.blob_doh(img)
+    array([[197.      , 153.      ,  20.333334],
+           [124.      , 336.      ,  20.333334],
+           [126.      , 153.      ,  20.333334],
+           [195.      , 100.      ,  23.555555],
+           [192.      , 212.      ,  23.555555],
+           [121.      , 271.      ,  30.      ],
+           [126.      , 101.      ,  20.333334],
+           [193.      , 275.      ,  23.555555],
+           [123.      , 205.      ,  20.333334],
+           [270.      , 363.      ,  30.      ],
+           [265.      , 113.      ,  23.555555],
+           [262.      , 243.      ,  23.555555],
+           [185.      , 348.      ,  30.      ],
+           [156.      , 302.      ,  30.      ],
+           [123.      ,  44.      ,  23.555555],
+           [260.      , 173.      ,  30.      ],
+           [197.      ,  44.      ,  20.333334]], dtype=float32)
+
+
+    Notes
+    -----
+    The radius of each blob is approximately `sigma`.
+    Computation of Determinant of Hessians is independent of the standard
+    deviation. Therefore detecting larger blobs won't take more time. In
+    methods line :py:meth:`blob_dog` and :py:meth:`blob_log` the computation
+    of Gaussians for larger `sigma` takes more time. The downside is that
+    this method can't be used for detecting blobs of radius less than `3px`
+    due to the box filters used in the approximation of Hessian Determinant.
+    """
+    check_nD(image, 2)
+
+    image = img_as_float(image)
+    float_dtype = _supported_float_type(image.dtype)
+    image = image.astype(float_dtype, copy=False)
+
+    image = integral_image(image)
+
+    if log_scale:
+        start, stop = math.log(min_sigma, 10), math.log(max_sigma, 10)
+        sigma_list = np.logspace(start, stop, num_sigma)
+    else:
+        sigma_list = np.linspace(min_sigma, max_sigma, num_sigma)
+
+    image_cube = cp.empty(image.shape + (len(sigma_list),), dtype=float_dtype)
+    for i, s in enumerate(sigma_list):
+        image_cube[..., i] = _hessian_matrix_det(image, s)
+
+    local_maxima = peak_local_max(image_cube,
+                                  threshold_abs=threshold,
+                                  threshold_rel=threshold_rel,
+                                  exclude_border=False,
+                                  footprint=cp.ones((3,) * image_cube.ndim))
+
+    # Catch no peaks
+    if local_maxima.size == 0:
+        return cp.empty((0, 3))
+    # Convert local_maxima to float type of input image
+    lm = local_maxima.astype(float_dtype)
+    # Convert the last index to its corresponding scale value
+    sigma_list = cp.asarray(sigma_list)
+    lm[:, -1] = sigma_list[local_maxima[:, -1]]
+    return _prune_blobs(lm, overlap)
diff --git a/python/cucim/src/cucim/skimage/feature/corner.py b/python/cucim/src/cucim/skimage/feature/corner.py
index a07fcfaa0..0863f94c6 100644
--- a/python/cucim/src/cucim/skimage/feature/corner.py
+++ b/python/cucim/src/cucim/skimage/feature/corner.py
@@ -12,6 +12,8 @@
 # from ..transform import integral_image
 from .._shared._gradient import gradient
 from .._shared.utils import _supported_float_type, warn
+from ..transform import integral_image
+from ._hessian_det_appx import _hessian_matrix_det
 from .peak import peak_local_max
 from .util import _prepare_grayscale_input_nD
 
@@ -383,9 +385,8 @@ def hessian_matrix_det(image, sigma=1, approximate=True):
     float_dtype = _supported_float_type(image.dtype)
     image = image.astype(float_dtype, copy=False)
     if image.ndim == 2 and approximate:
-        raise NotImplementedError("approximate=True case not implemented")
-        # integral = integral_image(image)
-        # return cp.asarray(_hessian_matrix_det(integral, sigma))
+        integral = integral_image(image)
+        return cp.asarray(_hessian_matrix_det(integral, sigma))
     else:  # slower brute-force implementation for nD images
         if image.ndim in [2, 3]:
             # Compute determinant as the product of the eigenvalues.
@@ -752,7 +753,7 @@ def hessian_matrix_eigvals(H_elems):
     --------
     >>> import cupy as cp
     >>> from cucim.skimage.feature import (hessian_matrix,
-    ...                                      hessian_matrix_eigvals)
+    ...                                    hessian_matrix_eigvals)
     >>> square = cp.zeros((5, 5))
     >>> square[2, 2] = 4
     >>> H_elems = hessian_matrix(square, sigma=0.1, order='rc',
diff --git a/python/cucim/src/cucim/skimage/feature/cuda/_hessian_det_appx.cu b/python/cucim/src/cucim/skimage/feature/cuda/_hessian_det_appx.cu
new file mode 100644
index 000000000..ace160602
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/feature/cuda/_hessian_det_appx.cu
@@ -0,0 +1,191 @@
+/***********************************************************************
+
+_hessian_det_appx.cu has the following scikit-image (v0.19.3) functions as cuda functions:
+    - _clip
+    - _integ
+    - _hessian_matrix_det
+
+Helpful example to get into ElementWiseKernel:
+    - https://github.com/yuyu2172/cupy/blob/a651541668c8dd024973526dbd5ed672cc482260/examples/gemm/README.md
+***********************************************************************/
+
+
+/* @func _clip
+    Clip coordinate between low and high values.
+    This method was created so that `hessian_det_appx` does not have to make
+    a Python call.
+    Parameters
+    ----------
+    x : int
+        Coordinate to be clipped.
+    low : int
+        The lower bound.
+    high : int
+        The higher bound.
+    Returns
+    -------
+    x : int
+        `x` clipped between `high` and `low`.
+    """
+    assert 0 <= low <= high
+
+    if x > high:
+        return high
+    elif x < low:
+        return low
+    else:
+        return x
+*/
+__device__
+inline int _clip(const int x, const int low, const int high)
+{
+//     assert((int)((0 <= low) && (low <= high));
+
+    if (x > high)
+    {
+        return high;
+    }
+    else if (x < low)
+    {
+        return low;
+    }
+    else
+    {
+        return x;
+    }
+}
+
+
+/*@func _integ
+    Integrate over the 2D integral image in the given window.
+    This method was created so that `hessian_det_appx` does not have to make
+    a Python call.
+    Parameters
+    ----------
+    img : array
+        The integral image over which to integrate.
+    r : int
+        The row number of the top left corner.
+    c : int
+        The column number of the top left corner.
+    rl : int
+        The number of rows over which to integrate.
+    cl : int
+        The number of columns over which to integrate.
+    Returns
+    -------
+    ans : double
+        The integral over the given window.
+*/
+__device__
+inline IMAGE_T _integ(const IMAGE_T * img,
+                     const INT_T img_rows,
+                     const INT_T img_cols,
+                     int r,
+                     int c,
+                     const int rl,
+                     const int cl)
+{
+
+    r = _clip(r, 0, img_rows - 1);
+    c = _clip(c, 0, img_cols - 1);
+
+    const int r2 = _clip(r + rl, 0, img_rows - 1);
+    const int c2 = _clip(c + cl, 0, img_cols - 1);
+
+    IMAGE_T ans = img[r * img_cols + c] + img[r2 * img_cols + c2] - img[r * img_cols + c2] - img[r2 * img_cols + c];
+
+    return max(0., ans);
+}
+
+
+/*@func _hessian_matrix_det
+    Compute the approximate Hessian Determinant over a 2D image.
+    This method uses box filters over integral images to compute the
+    approximate Hessian Determinant as described in [1]_.
+    Parameters
+    ----------
+    img : array
+        The integral image over which to compute Hessian Determinant.
+    sigma : double
+        Standard deviation used for the Gaussian kernel, used for the Hessian
+        matrix
+    Returns
+    -------
+    out : array
+        The array of the Determinant of Hessians.
+    References
+    ----------
+    .. [1] Herbert Bay, Andreas Ess, Tinne Tuytelaars, Luc Van Gool,
+           "SURF: Speeded Up Robust Features"
+           ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf
+    Notes
+    -----
+    The running time of this method only depends on size of the image. It is
+    independent of `sigma` as one would expect. The downside is that the
+    result for `sigma` less than `3` is not accurate, i.e., not similar to
+    the result obtained if someone computed the Hessian and took its
+    determinant.
+*/
+extern "C" __global__
+void _hessian_matrix_det(
+        const IMAGE_T* img,
+        const INT_T img_rows,
+        const INT_T img_cols,
+        const double sigma,
+        IMAGE_T* out)
+{
+    // Some note her: Input type from cupy must be a cp.float64, which maps to cuda's double
+    // If someone wants to fasten this function with cp.float32, this function must be re-written with cuda's float
+
+    // *************************************************************************
+    // This function is derived from Scikit-Image _hessian_matrix_det (v0.19.3):
+    // *************************************************************************
+
+    int tid = blockDim.x * blockIdx.x + threadIdx.x;
+
+    if (tid >= img_rows*img_cols)
+    {
+        return;
+    }
+
+    const int r = (int) floorf(tid / img_cols);
+    const int c = (int) tid % img_cols;
+
+    int size = (int)(3.0 * sigma);
+
+    const int s2 = floorf((size - 1) / 2);
+    const int s3 = floorf(size / 3);
+    const int l = floorf(size / 3);
+    const int w = size;
+    const int b = floorf((size - 1) / 2);
+
+    const IMAGE_T w_i = 1.0 / size / size;
+
+    if ((size % 2) == 0)
+    {
+        size += 1;
+    }
+
+    const IMAGE_T tl = _integ(img, img_rows, img_cols, r - s3, c - s3, s3, s3);  // top left
+    const IMAGE_T br = _integ(img, img_rows, img_cols, r + 1, c + 1, s3, s3);  // bottom right
+    const IMAGE_T bl = _integ(img, img_rows, img_cols, r - s3, c + 1, s3, s3);  // bottom left
+    const IMAGE_T tr = _integ(img, img_rows, img_cols, r + 1, c - s3, s3, s3);  // top right
+
+    IMAGE_T dxy = bl + tr - tl - br;
+    dxy = -dxy * w_i;
+
+    IMAGE_T mid = _integ(img, img_rows, img_cols, r - s3 + 1, c - s2, 2 * s3 - 1, w);  // middle box
+    IMAGE_T side = _integ(img, img_rows, img_cols, r - s3 + 1, c - (int)floorf(s3 / 2), 2 * s3 - 1, s3);  // sides
+
+    IMAGE_T dxx = mid - 3 * side;
+    dxx = -dxx * w_i;
+
+    mid = _integ(img, img_rows, img_cols, r - s2, c - s3 + 1, w, 2 * s3 - 1);
+    side = _integ(img, img_rows, img_cols, r - (int)floorf(s3 / 2), c - s3 + 1, s3, 2 * s3 - 1);
+
+    IMAGE_T dyy = mid - 3 * side;
+    dyy = -dyy * w_i;
+
+    out[tid] = (dxx * dyy - 0.81 * (dxy * dxy));
+}
diff --git a/python/cucim/src/cucim/skimage/feature/cuda/blob.cu b/python/cucim/src/cucim/skimage/feature/cuda/blob.cu
new file mode 100644
index 000000000..eddae9775
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/feature/cuda/blob.cu
@@ -0,0 +1,254 @@
+/***********************************************************************
+
+blob.cu has the following scikit-image functions as cuda functions:
+    - _compute_disk_overlap
+    - _compute_sphere_overlap
+    - _blob_overlap
+    - _prune_blobs
+
+Helpful example to get into ElementWiseKernel:
+    - https://github.com/yuyu2172/cupy/blob/a651541668c8dd024973526dbd5ed672cc482260/examples/gemm/README.md
+***********************************************************************/
+
+
+/*@func _compute_disk_overlap
+    Compute fraction of surface overlap between two disks of radii
+    ``r1`` and ``r2``, with centers separated by a distance ``d``.
+
+    Parameters
+    ----------
+    d : float
+        Distance between centers.
+    r1 : float
+        Radius of the first disk.
+    r2 : float
+        Radius of the second disk.
+
+    Returns
+    -------
+    fraction: float
+        Fraction of area of the overlap between the two disks.
+*/
+__device__
+inline BLOB_T _compute_disk_overlap(const BLOB_T d, const BLOB_T r1, const BLOB_T r2)
+{
+
+    const BLOB_T pi = 3.141592653589793;
+    BLOB_T ratio1 = 0.0;
+    BLOB_T ratio2 = 0.0;
+    BLOB_T r = 0.0;
+
+
+    ratio1 = (d*d + r1*r1 - r2*r2) / (2 * d * r1);
+    if (ratio1 > 1.0)
+    {
+        ratio1 = 1.0;
+    }
+    else if (ratio1 < -1.0)
+    {
+        ratio1 = -1.0;
+    }
+
+    ratio2 = (d*d + r2*r2 - r1*r1) / (2 * d * r2);
+    if (ratio2 > 1.0)
+    {
+        ratio2 = 1.0;
+    }
+    else if (ratio2 < -1.0)
+    {
+        ratio2 = -1.0;
+    }
+
+    BLOB_T a = -d + r2 + r1;
+    BLOB_T b = d - r2 + r1;
+    BLOB_T c = d + r2 - r1;
+    BLOB_T _d = d + r2 + r1;
+
+    r = min(r1, r2);
+    return (r1*r1 * acos(ratio1) + r2*r2 * acos(ratio2) - 0.5 * sqrt(abs(a * b * c * _d))) / (pi * r*r);
+}
+
+
+/*@func _compute_sphere_overlap
+    Compute volume overlap fraction between two spheres of radii
+    ``r1`` and ``r2``, with centers separated by a distance ``d``.
+
+    Parameters
+    ----------
+    d : float
+        Distance between centers.
+    r1 : float
+        Radius of the first sphere.
+    r2 : float
+        Radius of the second sphere.
+
+    Returns
+    -------
+    fraction: float
+        Fraction of volume of the overlap between the two spheres.
+
+    Notes
+    -----
+    See for example http://mathworld.wolfram.com/Sphere-SphereIntersection.html
+    for more details.
+*/
+__device__
+inline BLOB_T _compute_sphere_overlap(const BLOB_T d, const BLOB_T r1, const BLOB_T r2) {
+
+    const BLOB_T pi = 3.141592653589793;
+    const BLOB_T r = min(r1, r2);
+    const BLOB_T vol = (pi / (12 * d) * (r1 + r2 - d) * (r1 + r2 - d) * \
+                      (d*d + 2 * d * (r1 + r2) - 3 * (r1*r1 + r2*r2) + 6 * r1 * r2));
+    return vol / (4.0/3.0 * pi * r*r*r);
+}
+
+
+/*@func _blob_overlap
+    Finds the overlapping area fraction between two blobs.
+
+    Returns a float representing fraction of overlapped area. Note that 0.0
+    is *always* returned for dimension greater than 3.
+
+    Parameters
+    ----------
+    blob1 : sequence of arrays
+        A sequence of ``(row, col, sigma)`` or ``(pln, row, col, sigma)``,
+        where ``row, col`` (or ``(pln, row, col)``) are coordinates
+        of blob and ``sigma`` is the standard deviation of the Gaussian kernel
+        which detected the blob.
+    blob2 : sequence of arrays
+        A sequence of ``(row, col, sigma)`` or ``(pln, row, col, sigma)``,
+        where ``row, col`` (or ``(pln, row, col)``) are coordinates
+        of blob and ``sigma`` is the standard deviation of the Gaussian kernel
+        which detected the blob.
+    sigma_dim : int, optional
+        The dimensionality of the sigma value. Can be 1 or the same as the
+        dimensionality of the blob space (2 or 3).
+
+    Returns
+    -------
+    f : float
+        Fraction of overlapped area (or volume in 3D).
+
+*/
+__device__
+inline BLOB_T _blob_overlap(BLOB_T* blob1, BLOB_T* blob2, const INT_T sigma_dim, const INT_T n) {
+
+    const INT_T ndim = n - sigma_dim;
+    BLOB_T pos1 = 0.0;
+    BLOB_T pos2 = 0.0;
+    BLOB_T _sum = 0.0;
+    BLOB_T r1 = 0.0;
+    BLOB_T r2 = 0.0;
+    BLOB_T d = 0.0;
+
+    //printf("blob1: [%3.1f %3.1f %3.1f]\n", blob1[0], blob1[1], blob1[2]);
+    //printf("blob2: [%3.1f %3.1f %3.1f]\n", blob2[0], blob2[1], blob2[2]);
+
+    if (ndim > 3)
+    {
+        return 0.0;
+    }
+    const BLOB_T root_ndim = sqrt((BLOB_T)ndim);
+
+    if (blob1[n-1]== 0.0  && blob2[n-1]== 0.0)
+    {
+        return 0.0;
+    }
+    else if (blob1[n-1] > blob2[n-1])
+    {
+        r1 = 1.0;
+        r2 = blob2[n-1] / blob1[n-1];
+        for (INT_T _i=0; _i < ndim; _i++)
+        {
+            pos1 = blob1[_i] / (blob1[min(n-sigma_dim+_i, n-1)] * root_ndim);
+            pos2 = blob2[_i] / (blob1[min(n-sigma_dim+_i, n-1)] * root_ndim);
+            _sum += ((pos2 - pos1) * (pos2 - pos1));
+        }
+    }
+    else
+    {
+
+        r1 = blob1[n-1] / blob2[n-1];
+        r2 = 1.0;
+        for (INT_T _i=0; _i < ndim; _i++)
+        {
+            pos1 = blob1[_i] / (blob2[min(n-sigma_dim+_i, n-1)] * root_ndim);
+            pos2 = blob2[_i] / (blob2[min(n-sigma_dim+_i, n-1)] * root_ndim);
+            _sum += ((pos2 - pos1) * (pos2 - pos1));
+        }
+    }
+
+    d = sqrt(_sum);
+
+    if (d > (r1 + r2)) // centers farther than sum of radii, so no overlap
+    {
+        return 0.0;
+    }
+
+    // one blob is inside the other
+    if (d <= abs(r1 - r2))
+    {
+        return 1.0;
+    }
+    if (ndim == 2)
+    {
+        return _compute_disk_overlap(d, r1, r2);
+    }
+    else // ndim=3 http://mathworld.wolfram.com/Sphere-SphereIntersection.html
+    {
+        return _compute_sphere_overlap(d, r1, r2);
+    }
+}
+
+
+extern "C" __global__
+void _prune_blobs(
+        BLOB_T * blobs_array,
+        const INT_T n_rows,
+        const INT_T n_cols,
+        const double overlap,
+        INT_T sigma_dim)
+{
+
+    // *************************************************************************
+    // This function is derived from Scikit-Image _prune_blobs (v0.11.x):
+    // *************************************************************************
+
+    INT_T tid = blockDim.x * blockIdx.x + threadIdx.x;
+    BLOB_T *blob1;
+    BLOB_T *blob2;
+
+    if (tid >= n_rows)
+    {
+        return;  // all done
+    }
+
+    for(INT_T k=0; k<n_rows; k++)
+    {
+        // blob[tid] --> blob1
+        // blob[k  ] --> blob2
+
+        if (tid >= k)
+        {
+            // skip calculation, it is done already
+            continue;
+        }
+
+        blob1 = &blobs_array[(INT_T)tid*n_cols];
+        blob2 = &blobs_array[(INT_T)k*n_cols];
+
+        if (_blob_overlap(blob1, blob2, sigma_dim, n_cols) > overlap)
+        {
+            // note: this test works even in the anisotropic case because all sigmas increase together.
+            if (blob1[n_cols-1] > blob2[n_cols-1])
+            {
+                blob2[n_cols-1] = 0.0;
+            }
+            else
+            {
+                blob1[n_cols-1] = 0.0;
+            }
+        }
+    }
+}


Test patch:
diff --git a/python/cucim/src/cucim/skimage/feature/tests/test_blob.py b/python/cucim/src/cucim/skimage/feature/tests/test_blob.py
new file mode 100644
index 000000000..d1a4041cc
--- /dev/null
+++ b/python/cucim/src/cucim/skimage/feature/tests/test_blob.py
@@ -0,0 +1,529 @@
+import math
+
+import cupy as cp
+import pytest
+
+from skimage.draw import disk
+from skimage.draw.draw3d import ellipsoid
+
+from cucim.skimage import feature
+from cucim.skimage.feature import blob_dog, blob_doh, blob_log
+
+
+@pytest.mark.parametrize(
+    'dtype', [cp.uint8, cp.float16, cp.float32, cp.float64]
+)
+@pytest.mark.parametrize('threshold_type', ['absolute', 'relative'])
+def test_blob_dog(dtype, threshold_type):
+    img = cp.ones((512, 512), dtype=dtype)
+
+    xs, ys = cp.asarray(disk((400, 130), 5))
+    img[xs, ys] = 255
+
+    xs, ys = cp.asarray(disk((100, 300), 25))
+    img[xs, ys] = 255
+
+    xs, ys = cp.asarray(disk((200, 350), 45))
+    img[xs, ys] = 255
+
+    if threshold_type == 'absolute':
+        threshold = 2.0
+        if img.dtype.kind != 'f':
+            # account for internal scaling to [0, 1] by img_as_float
+            threshold /= img.ptp()
+        threshold_rel = None
+    elif threshold_type == 'relative':
+        threshold = None
+        threshold_rel = 0.5
+
+    blobs = blob_dog(
+        img,
+        min_sigma=4,
+        max_sigma=50,
+        threshold=threshold,
+        threshold_rel=threshold_rel,
+    )
+
+    def radius(x):
+        return math.sqrt(2) * x[2]
+
+    s = sorted(blobs, key=radius)
+    thresh = 5
+    ratio_thresh = 0.25
+
+    b = s[0]
+    assert abs(b[0] - 400) <= thresh
+    assert abs(b[1] - 130) <= thresh
+    assert abs(radius(b) - 5) <= ratio_thresh * 5
+
+    b = s[1]
+    assert abs(b[0] - 100) <= thresh
+    assert abs(b[1] - 300) <= thresh
+    assert abs(radius(b) - 25) <= ratio_thresh * 25
+
+    b = s[2]
+    assert abs(b[0] - 200) <= thresh
+    assert abs(b[1] - 350) <= thresh
+    assert abs(radius(b) - 45) <= ratio_thresh * 45
+
+    # Testing no peaks
+    img_empty = cp.zeros((100, 100), dtype=dtype)
+    assert blob_dog(img_empty).size == 0
+
+
+@pytest.mark.parametrize(
+    'dtype', [cp.uint8, cp.float16, cp.float32, cp.float64]
+)
+@pytest.mark.parametrize('threshold_type', ['absolute', 'relative'])
+def test_blob_dog_3d(dtype, threshold_type):
+    # Testing 3D
+    r = 10
+    pad = 10
+    im3 = cp.asarray(ellipsoid(r, r, r))
+    im3 = cp.pad(im3, pad, mode='constant')
+
+    if threshold_type == 'absolute':
+        threshold = 0.001
+        threshold_rel = 0
+    elif threshold_type == 'relative':
+        threshold = 0
+        threshold_rel = 0.5
+
+    blobs = blob_dog(
+        im3,
+        min_sigma=3,
+        max_sigma=10,
+        sigma_ratio=1.2,
+        threshold=threshold,
+        threshold_rel=threshold_rel,
+    )
+    b = blobs[0]
+
+    assert b.shape == (4,)
+    assert b[0] == r + pad + 1
+    assert b[1] == r + pad + 1
+    assert b[2] == r + pad + 1
+    assert abs(math.sqrt(3) * b[3] - r) < 1.1
+
+
+@pytest.mark.parametrize(
+    'dtype', [cp.uint8, cp.float16, cp.float32, cp.float64]
+)
+@pytest.mark.parametrize('threshold_type', ['absolute', 'relative'])
+def test_blob_dog_3d_anisotropic(dtype, threshold_type):
+    # Testing 3D anisotropic
+    r = 10
+    pad = 10
+    im3 = cp.asarray(ellipsoid(r / 2, r, r))
+    im3 = cp.pad(im3, pad, mode='constant')
+
+    if threshold_type == 'absolute':
+        threshold = 0.001
+        threshold_rel = None
+    elif threshold_type == 'relative':
+        threshold = None
+        threshold_rel = 0.5
+
+    blobs = blob_dog(
+        im3.astype(dtype, copy=False),
+        min_sigma=[1.5, 3, 3],
+        max_sigma=[5, 10, 10],
+        sigma_ratio=1.2,
+        threshold=threshold,
+        threshold_rel=threshold_rel,
+    )
+    b = blobs[0]
+
+    assert b.shape == (6,)
+    assert b[0] == r / 2 + pad + 1
+    assert b[1] == r + pad + 1
+    assert b[2] == r + pad + 1
+    assert abs(math.sqrt(3) * b[3] - r / 2) < 1.1
+    assert abs(math.sqrt(3) * b[4] - r) < 1.1
+    assert abs(math.sqrt(3) * b[5] - r) < 1.1
+
+
+def test_blob_dog_excl_border():
+    # Testing exclude border
+
+    # image where blob is 5 px from borders, radius 5
+    img = cp.ones((512, 512))
+    xs, ys = cp.asarray(disk((5, 5), 5))
+    img[xs, ys] = 255
+    blobs = blob_dog(
+        img,
+        min_sigma=1.5,
+        max_sigma=5,
+        sigma_ratio=1.2,
+    )
+    assert blobs.shape[0] == 1
+    b = blobs[0]
+    assert b[0] == b[1] == 5, "blob should be 5 px from x and y borders"
+
+    blobs = blob_dog(
+        img,
+        min_sigma=1.5,
+        max_sigma=5,
+        sigma_ratio=1.2,
+        exclude_border=6,
+    )
+    msg = "zero blobs should be detected, as only blob is 5 px from border"
+    assert blobs.shape[0] == 0, msg
+
+
+@pytest.mark.parametrize('anisotropic', [False, True])
+@pytest.mark.parametrize('ndim', [1, 2, 3, 4])
+@pytest.mark.parametrize('function_name', ['blob_dog', 'blob_log'])
+def test_nd_blob_no_peaks_shape(function_name, ndim, anisotropic):
+    # uniform image so no blobs will be found
+    z = cp.zeros((16,) * ndim, dtype=cp.float32)
+    if anisotropic:
+        max_sigma = 8 + cp.arange(ndim)
+    else:
+        max_sigma = 8
+    blob_func = getattr(feature, function_name)
+    blobs = blob_func(z, max_sigma=max_sigma)
+    # z.ndim +  (z.ndim sigmas if anisotropic, only one sigma otherwise)
+    expected_shape = 2 * z.ndim if anisotropic else z.ndim + 1
+    assert blobs.shape == (0, expected_shape)
+
+
+@pytest.mark.parametrize(
+    'dtype', [cp.uint8, cp.float16, cp.float32, cp.float64]
+)
+@pytest.mark.parametrize('threshold_type', ['absolute', 'relative'])
+def test_blob_log(dtype, threshold_type):
+    img = cp.ones((256, 256), dtype=dtype)
+
+    xs, ys = cp.asarray(disk((200, 65), 5))
+    img[xs, ys] = 255
+
+    xs, ys = cp.asarray(disk((80, 25), 15))
+    img[xs, ys] = 255
+
+    xs, ys = cp.asarray(disk((50, 150), 25))
+    img[xs, ys] = 255
+
+    xs, ys = cp.asarray(disk((100, 175), 30))
+    img[xs, ys] = 255
+
+    if threshold_type == 'absolute':
+        threshold = 1
+        if img.dtype.kind != 'f':
+            # account for internal scaling to [0, 1] by img_as_float
+            threshold /= img.ptp()
+        threshold_rel = None
+    elif threshold_type == 'relative':
+        threshold = None
+        threshold_rel = 0.5
+
+    blobs = blob_log(img, min_sigma=5, max_sigma=20, threshold=threshold,
+                     threshold_rel=threshold_rel)
+
+    def radius(x):
+        return math.sqrt(2) * x[2]
+
+    s = sorted(blobs, key=radius)
+    thresh = 3
+
+    b = s[0]
+    assert abs(b[0] - 200) <= thresh
+    assert abs(b[1] - 65) <= thresh
+    assert abs(radius(b) - 5) <= thresh
+
+    b = s[1]
+    assert abs(b[0] - 80) <= thresh
+    assert abs(b[1] - 25) <= thresh
+    assert abs(radius(b) - 15) <= thresh
+
+    b = s[2]
+    assert abs(b[0] - 50) <= thresh
+    assert abs(b[1] - 150) <= thresh
+    assert abs(radius(b) - 25) <= thresh
+
+    b = s[3]
+    assert abs(b[0] - 100) <= thresh
+    assert abs(b[1] - 175) <= thresh
+    assert abs(radius(b) - 30) <= thresh
+
+    # Testing log scale
+    blobs = blob_log(
+        img,
+        min_sigma=5,
+        max_sigma=20,
+        threshold=threshold,
+        threshold_rel=threshold_rel,
+        log_scale=True)
+
+    b = s[0]
+    assert abs(b[0] - 200) <= thresh
+    assert abs(b[1] - 65) <= thresh
+    assert abs(radius(b) - 5) <= thresh
+
+    b = s[1]
+    assert abs(b[0] - 80) <= thresh
+    assert abs(b[1] - 25) <= thresh
+    assert abs(radius(b) - 15) <= thresh
+
+    b = s[2]
+    assert abs(b[0] - 50) <= thresh
+    assert abs(b[1] - 150) <= thresh
+    assert abs(radius(b) - 25) <= thresh
+
+    b = s[3]
+    assert abs(b[0] - 100) <= thresh
+    assert abs(b[1] - 175) <= thresh
+    assert abs(radius(b) - 30) <= thresh
+
+    # Testing no peaks
+    img_empty = cp.zeros((100, 100))
+    assert blob_log(img_empty).size == 0
+
+
+def test_blob_log_no_warnings():
+    img = cp.ones((11, 11))
+
+    xs, ys = cp.asarray(disk((5, 5), 2))
+    img[xs, ys] = 255
+
+    xs, ys = cp.asarray(disk((7, 6), 2))
+    img[xs, ys] = 255
+
+    blob_log(img, max_sigma=20, num_sigma=10, threshold=.1)
+
+
+def test_blob_log_3d():
+    # Testing 3D
+    r = 6
+    pad = 10
+    im3 = cp.asarray(ellipsoid(r, r, r))
+    im3 = cp.pad(im3, pad, mode='constant')
+
+    blobs = blob_log(im3, min_sigma=3, max_sigma=10)
+    b = blobs[0]
+
+    assert b.shape == (4,)
+    assert b[0] == r + pad + 1
+    assert b[1] == r + pad + 1
+    assert b[2] == r + pad + 1
+    assert abs(math.sqrt(3) * b[3] - r) < 1
+
+
+def test_blob_log_3d_anisotropic():
+    # Testing 3D anisotropic
+    r = 6
+    pad = 10
+    im3 = cp.asarray(ellipsoid(r / 2, r, r))
+    im3 = cp.pad(im3, pad, mode='constant')
+
+    blobs = blob_log(
+        im3,
+        min_sigma=[1, 2, 2],
+        max_sigma=[5, 10, 10],
+    )
+
+    b = blobs[0]
+    assert b.shape == (6,)
+    assert b[0] == r / 2 + pad + 1
+    assert b[1] == r + pad + 1
+    assert b[2] == r + pad + 1
+    assert abs(math.sqrt(3) * b[3] - r / 2) < 1
+    assert abs(math.sqrt(3) * b[4] - r) < 1
+    assert abs(math.sqrt(3) * b[5] - r) < 1
+
+
+def test_blob_log_exclude_border():
+    # image where blob is 5 px from borders, radius 5
+    img = cp.ones((512, 512))
+    xs, ys = cp.asarray(disk((5, 5), 5))
+    img[xs, ys] = 255
+
+    blobs = blob_log(
+        img,
+        min_sigma=1.5,
+        max_sigma=5,
+    )
+    assert blobs.shape[0] == 1
+    b = blobs[0]
+    assert b[0] == b[1] == 5, "blob should be 5 px from x and y borders"
+
+    blobs = blob_log(
+        img,
+        min_sigma=1.5,
+        max_sigma=5,
+        exclude_border=6,
+    )
+    msg = "zero blobs should be detected, as only blob is 5 px from border"
+    assert blobs.shape[0] == 0, msg
+
+
+@pytest.mark.parametrize("dtype", [cp.uint8, cp.float16, cp.float32])
+@pytest.mark.parametrize('threshold_type', ['absolute', 'relative'])
+def test_blob_doh(dtype, threshold_type):
+    img = cp.ones((512, 512), dtype=dtype)
+
+    xs, ys = cp.asarray(disk((400, 130), 20))
+    img[xs, ys] = 255
+
+    xs, ys = cp.asarray(disk((460, 50), 30))
+    img[xs, ys] = 255
+
+    xs, ys = cp.asarray(disk((100, 300), 40))
+    img[xs, ys] = 255
+
+    xs, ys = cp.asarray(disk((200, 350), 50))
+    img[xs, ys] = 255
+
+    if threshold_type == 'absolute':
+        # Note: have to either scale up threshold or rescale the image to the
+        #       range [0, 1] internally.
+        threshold = 0.05
+        if img.dtype.kind == 'f':
+            # account for lack of internal scaling to [0, 1] by img_as_float
+            ptp = img.ptp()
+            threshold *= ptp ** 2
+        threshold_rel = None
+    elif threshold_type == 'relative':
+        threshold = None
+        threshold_rel = 0.5
+
+    blobs = blob_doh(
+        img,
+        min_sigma=1,
+        max_sigma=60,
+        num_sigma=10,
+        threshold=threshold,
+        threshold_rel=threshold_rel)
+
+    def radius(x):
+        return x[2]
+
+    s = sorted(blobs, key=radius)
+    thresh = 4
+
+    b = s[0]
+    assert abs(b[0] - 400) <= thresh
+    assert abs(b[1] - 130) <= thresh
+    assert abs(radius(b) - 20) <= thresh
+
+    b = s[1]
+    assert abs(b[0] - 460) <= thresh
+    assert abs(b[1] - 50) <= thresh
+    assert abs(radius(b) - 30) <= thresh
+
+    b = s[2]
+    assert abs(b[0] - 100) <= thresh
+    assert abs(b[1] - 300) <= thresh
+    assert abs(radius(b) - 40) <= thresh
+
+    b = s[3]
+    assert abs(b[0] - 200) <= thresh
+    assert abs(b[1] - 350) <= thresh
+    assert abs(radius(b) - 50) <= thresh
+
+
+def test_blob_doh_log_scale():
+    img = cp.ones((512, 512), dtype=cp.uint8)
+
+    xs, ys = cp.asarray(disk((400, 130), 20))
+    img[xs, ys] = 255
+
+    xs, ys = cp.asarray(disk((460, 50), 30))
+    img[xs, ys] = 255
+
+    xs, ys = cp.asarray(disk((100, 300), 40))
+    img[xs, ys] = 255
+
+    xs, ys = cp.asarray(disk((200, 350), 50))
+    img[xs, ys] = 255
+
+    blobs = blob_doh(
+        img,
+        min_sigma=1,
+        max_sigma=60,
+        num_sigma=10,
+        log_scale=True,
+        threshold=.05)
+
+    def radius(x):
+        return x[2]
+
+    s = sorted(blobs, key=radius)
+    thresh = 10
+
+    b = s[0]
+    assert abs(b[0] - 400) <= thresh
+    assert abs(b[1] - 130) <= thresh
+    assert abs(radius(b) - 20) <= thresh
+
+    b = s[2]
+    assert abs(b[0] - 460) <= thresh
+    assert abs(b[1] - 50) <= thresh
+    assert abs(radius(b) - 30) <= thresh
+
+    b = s[1]
+    assert abs(b[0] - 100) <= thresh
+    assert abs(b[1] - 300) <= thresh
+    assert abs(radius(b) - 40) <= thresh
+
+    b = s[3]
+    assert abs(b[0] - 200) <= thresh
+    assert abs(b[1] - 350) <= thresh
+    assert abs(radius(b) - 50) <= thresh
+
+
+def test_blob_doh_no_peaks():
+    # Testing no peaks
+    img_empty = cp.zeros((100, 100))
+    assert blob_doh(img_empty).size == 0
+
+
+def test_blob_doh_overlap():
+    img = cp.ones((256, 256), dtype=cp.uint8)
+
+    xs, ys = cp.asarray(disk((100, 100), 20))
+    img[xs, ys] = 255
+
+    xs, ys = cp.asarray(disk((120, 100), 30))
+    img[xs, ys] = 255
+
+    blobs = blob_doh(
+        img,
+        min_sigma=1,
+        max_sigma=60,
+        num_sigma=10,
+        threshold=.05
+    )
+
+    assert len(blobs) == 1
+
+
+def test_blob_log_overlap_3d():
+    r1, r2 = 7, 6
+    pad1, pad2 = 11, 12
+    blob1 = cp.asarray(ellipsoid(r1, r1, r1))
+    blob1 = cp.pad(blob1, pad1, mode='constant')
+    blob2 = cp.asarray(ellipsoid(r2, r2, r2))
+    blob2 = cp.pad(blob2, [(pad2, pad2), (pad2 - 9, pad2 + 9),
+                           (pad2, pad2)],
+                   mode='constant')
+    im3 = cp.logical_or(blob1, blob2)
+
+    blobs = blob_log(im3, min_sigma=2, max_sigma=10, overlap=0.1)
+    assert len(blobs) == 1
+
+
+def test_blob_log_anisotropic():
+    image = cp.zeros((50, 50))
+    image[20, 10:20] = 1
+    isotropic_blobs = blob_log(image, min_sigma=0.5, max_sigma=2, num_sigma=3)
+    assert len(isotropic_blobs) > 1  # many small blobs found in line
+    ani_blobs = blob_log(image, min_sigma=[0.5, 5], max_sigma=[2, 20],
+                         num_sigma=3)  # 10x anisotropy, line is 1x10
+    assert len(ani_blobs) == 1  # single anisotropic blob found
+
+
+def test_no_blob():
+    im = cp.zeros((10, 10))
+    blobs = blob_log(im, min_sigma=2, max_sigma=5, num_sigma=4)
+    assert len(blobs) == 0
diff --git a/python/cucim/src/cucim/skimage/feature/tests/test_corner.py b/python/cucim/src/cucim/skimage/feature/tests/test_corner.py
index 6d2a5ae05..0c24ae05f 100644
--- a/python/cucim/src/cucim/skimage/feature/tests/test_corner.py
+++ b/python/cucim/src/cucim/skimage/feature/tests/test_corner.py
@@ -275,11 +275,12 @@ def test_custom_eigvals_kernels_vs_linalg_eigvalsh(shape, dtype):
     cp.testing.assert_allclose(evs1, evs2, atol=atol)
 
 
-def test_hessian_matrix_det():
+@pytest.mark.parametrize('approximate', [False, True])
+def test_hessian_matrix_det(approximate):
     image = cp.zeros((5, 5))
     image[2, 2] = 1
     # TODO: approximate=True case not implemented
-    det = hessian_matrix_det(image, 5, approximate=False)
+    det = hessian_matrix_det(image, 5, approximate=approximate)
     assert_array_almost_equal(det, 0, decimal=3)