Skip to content
Snippets Groups Projects
Commit 25675ca4 authored by fima's avatar fima :beers:
Browse files

Merge branch 'blob_detection_refactoring' into 'main'

Refactoring of blob detection

See merge request !95
parents 527d197c ad5fb3f0
No related branches found
No related tags found
1 merge request!95Refactoring of blob detection
%% Cell type:code id: tags:
``` python
import qim3d
```
%% Output
WARNING:root:Could not load CuPy: No module named 'cupy'
%% Cell type:markdown id: tags:
### Blob detection notebook
%% Cell type:markdown id: tags:
This notebook shows how to do **blob detection** in a 3D volume using the `qim3d` library.
Blob detection is done by initializing a `qim3d.processing.Blob` object, and then calling the `qim3d.processing.Blob.detect` method. The `qim3d.processing.Blob.detect` method detects blobs by using the Difference of Gaussian (DoG) blob detection method, and returns an array `blobs` with the blobs found in the volume stored as `(p, r, c, radius)`. Subsequently, a binary mask of the volume can be retrieved with the `qim3d.processing.get_mask` method, in which the found blobs are marked as `True`.
Blob detection is done by using the `qim3d.processing.blob_detection` method, which detects blobs by using the Difference of Gaussian (DoG) blob detection method, and returns two arrays:
- `blobs`: The blobs found in the volume stored as `(p, r, c, radius)`
- `binary_volume`: A binary mask of the volume with the blobs marked as `True`
%% Cell type:markdown id: tags:
### **Example 1**: Blob detection in cement volume
%% Cell type:markdown id: tags:
**Applying Gaussian filter to volume**
%% Cell type:code id: tags:
``` python
# Import 3D volume of cement
cement = qim3d.examples.cement_128x128x128
# Visualize slices of the original cement volume
qim3d.viz.slices(cement, n_slices = 5, show = True)
# Apply Gaussian filter to the cement volume
cement_filtered = qim3d.processing.gaussian(cement, sigma = 2)
# Visualize slices of the filtered cement volume
qim3d.viz.slices(cement_filtered)
```
%% Output
<Figure size 1000x200 with 5 Axes>
%% Cell type:markdown id: tags:
**Detecting blobs in volume**
%% Cell type:code id: tags:
``` python
# Initialize blob detector
blob_detector = qim3d.processing.Blob(
background = "bright",
min_sigma = 1,
max_sigma = 8,
threshold = 0.001,
overlap = 0.1
# Detect blobs, and get binary mask
blobs, mask = qim3d.processing.blob_detection(
cement_filtered,
min_sigma=1,
max_sigma=8,
threshold=0.001,
overlap=0.1,
background="bright"
)
# Detect blobs in filtered volume
blobs = blob_detector.detect(vol = cement_filtered)
# Number of blobs found
print(f'Number of blobs found in the volume: {len(blobs)} blobs')
```
%% Output
Bright background selected, volume will be inverted.
Number of blobs found in the volume: 1813 blobs
%% Cell type:code id: tags:
``` python
# Visualize blobs on slices of cement volume
qim3d.viz.detection.circles(blobs, cement, show = True)
qim3d.viz.detection.circles(blobs, cement, alpha = 0.8, show = True, color = 'red')
```
%% Output
interactive(children=(IntSlider(value=64, description='Slice', max=127), Output()), layout=Layout(align_items=…
%% Cell type:markdown id: tags:
**Get binary mask of detected blobs**
**Binary mask of detected blobs**
%% Cell type:code id: tags:
``` python
# Get binary mask of detected blobs
mask = blob_detector.get_mask()
# Visualize mask
# Visualize binary mask
qim3d.viz.slicer(mask)
```
%% Output
interactive(children=(IntSlider(value=64, description='Slice', max=127), Output()), layout=Layout(align_items=…
......
......@@ -5,11 +5,16 @@ Here, we provide functionalities designed specifically for 3D image analysis and
::: qim3d.processing
options:
members:
- test_blob_detection
- blob_detection
- structure_tensor
- local_thickness
- get_3d_cc
- Pipeline
- Blob
- gaussian
- median
- maximum
- minimum
- tophat
::: qim3d.processing.operations
options:
......@@ -18,3 +23,7 @@ Here, we provide functionalities designed specifically for 3D image analysis and
- watershed
- edge_fade
- fade_mask
::: qim3d.processing
options:
members:
- Pipeline
\ No newline at end of file
"Testing docstring"
from .local_thickness_ import local_thickness
from .structure_tensor_ import structure_tensor
from .detection import blob_detection
from .filters import *
from .detection import *
from .operations import *
from .cc import get_3d_cc
""" Blob detection using Difference of Gaussian (DoG) method """
import numpy as np
from qim3d.io.logger import log
from skimage.feature import blob_dog
__all__ = ["Blob"]
class Blob:
def blob_detection(
vol: np.ndarray,
background: str = "dark",
min_sigma: float = 1,
max_sigma: float = 50,
sigma_ratio: float = 1.6,
threshold: float = 0.5,
overlap: float = 0.5,
threshold_rel: float = None,
exclude_border: bool = False,
) -> np.ndarray:
"""
Extract blobs from a volume using Difference of Gaussian (DoG) method
"""
def __init__(
self,
background="dark",
min_sigma=1,
max_sigma=50,
sigma_ratio=1.6,
threshold=0.5,
overlap=0.5,
threshold_rel=None,
exclude_border=False,
):
"""
Initialize the blob detection object
Args:
background: 'dark' if background is darker than the blobs, 'bright' if background is lighter than the blobs
min_sigma: The minimum standard deviation for Gaussian kernel
max_sigma: The maximum standard deviation for Gaussian kernel
sigma_ratio: The ratio between the standard deviation of Gaussian Kernels
threshold: The absolute lower bound for scale space maxima. Reduce this to detect blobs with lower intensities.
overlap: The fraction of area of two blobs that overlap
threshold_rel: The relative lower bound for scale space maxima
exclude_border: If True, exclude blobs that are too close to the border of the image
"""
self.background = background
self.min_sigma = min_sigma
self.max_sigma = max_sigma
self.sigma_ratio = sigma_ratio
self.threshold = threshold
self.overlap = overlap
self.threshold_rel = threshold_rel
self.exclude_border = exclude_border
self.vol_shape = None
self.blobs = None
def detect(self, vol):
"""
Detect blobs in the volume
Args:
vol: The volume to detect blobs in
Returns:
blobs: The blobs found in the volume as (p, r, c, radius)
Example:
Extract blobs from a volume using Difference of Gaussian (DoG) method, and retrieve a binary volume with the blobs marked as True
Args:
vol: The volume to detect blobs in
background: 'dark' if background is darker than the blobs, 'bright' if background is lighter than the blobs
min_sigma: The minimum standard deviation for Gaussian kernel
max_sigma: The maximum standard deviation for Gaussian kernel
sigma_ratio: The ratio between the standard deviation of Gaussian Kernels
threshold: The absolute lower bound for scale space maxima. Reduce this to detect blobs with lower intensities
overlap: The fraction of area of two blobs that overlap
threshold_rel: The relative lower bound for scale space maxima
exclude_border: If True, exclude blobs that are too close to the border of the image
Returns:
blobs: The blobs found in the volume as (p, r, c, radius)
binary_volume: A binary volume with the blobs marked as True
Example:
```python
import qim3d
......@@ -62,8 +42,9 @@ class Blob:
vol = qim3d.examples.cement_128x128x128
vol_blurred = qim3d.processing.gaussian(vol, sigma=2)
# Initialize Blob detector
blob_detector = qim3d.processing.Blob(
# Detect blobs, and get binary mask
blobs, mask = qim3d.processing.blob_detection(
vol_blurred,
min_sigma=1,
max_sigma=8,
threshold=0.001,
......@@ -71,88 +52,63 @@ class Blob:
background="bright"
)
# Detect blobs
blobs = blob_detector.detect(vol_blurred)
# Visualize results
qim3d.viz.circles(blobs,vol,alpha=0.8,color='blue')
# Visualize detected blobs
qim3d.viz.circles(blobs, vol, alpha=0.8, color='blue')
```
![blob detection](assets/screenshots/blob_detection.gif)
"""
self.vol_shape = vol.shape
if self.background == "bright":
log.info("Bright background selected, volume will be inverted.")
vol = np.invert(vol)
blobs = blob_dog(
vol,
min_sigma=self.min_sigma,
max_sigma=self.max_sigma,
sigma_ratio=self.sigma_ratio,
threshold=self.threshold,
overlap=self.overlap,
threshold_rel=self.threshold_rel,
exclude_border=self.exclude_border,
)
blobs[:, 3] = blobs[:, 3] * np.sqrt(3) # Change sigma to radius
self.blobs = blobs
return self.blobs
def get_mask(self):
'''
Retrieve a binary volume with the blobs marked as True
Returns:
binary_volume: A binary volume with the blobs marked as True
Example:
```python
import qim3d
# Get data
vol = qim3d.examples.cement_128x128x128
vol_blurred = qim3d.processing.gaussian(vol, sigma=2)
# Initialize Blob detector
blob_detector = qim3d.processing.Blob(
min_sigma=1,
max_sigma=8,
threshold=0.001,
overlap=0.1,
background="bright"
)
# Detect blobs
blobs = blob_detector.detect(vol_blurred)
# Get mask and visualize
mask = blob_detector.get_mask()
```python
# Visualize binary mask
qim3d.viz.slicer(mask)
```
![blob detection](assets/screenshots/blob_get_mask.gif)
'''
binary_volume = np.zeros(self.vol_shape, dtype=bool)
for z, y, x, radius in self.blobs:
# Calculate the bounding box around the blob
z_start = max(0, int(z - radius))
z_end = min(self.vol_shape[0], int(z + radius) + 1)
y_start = max(0, int(y - radius))
y_end = min(self.vol_shape[1], int(y + radius) + 1)
x_start = max(0, int(x - radius))
x_end = min(self.vol_shape[2], int(x + radius) + 1)
z_indices, y_indices, x_indices = np.indices((z_end - z_start, y_end - y_start, x_end - x_start))
z_indices += z_start
y_indices += y_start
x_indices += x_start
# Calculate distances from the center of the blob to voxels within the bounding box
dist = np.sqrt((x_indices - x)**2 + (y_indices - y)**2 + (z_indices - z)**2)
"""
binary_volume[z_start:z_end, y_start:y_end, x_start:x_end][dist <= radius] = True
if background == "bright":
log.info("Bright background selected, volume will be inverted.")
vol = np.invert(vol)
blobs = blob_dog(
vol,
min_sigma=min_sigma,
max_sigma=max_sigma,
sigma_ratio=sigma_ratio,
threshold=threshold,
overlap=overlap,
threshold_rel=threshold_rel,
exclude_border=exclude_border,
)
# Change sigma to radius
blobs[:, 3] = blobs[:, 3] * np.sqrt(3)
# Create binary mask of detected blobs
vol_shape = vol.shape
binary_volume = np.zeros(vol_shape, dtype=bool)
for z, y, x, radius in blobs:
# Calculate the bounding box around the blob
z_start = max(0, int(z - radius))
z_end = min(vol_shape[0], int(z + radius) + 1)
y_start = max(0, int(y - radius))
y_end = min(vol_shape[1], int(y + radius) + 1)
x_start = max(0, int(x - radius))
x_end = min(vol_shape[2], int(x + radius) + 1)
z_indices, y_indices, x_indices = np.indices(
(z_end - z_start, y_end - y_start, x_end - x_start)
)
z_indices += z_start
y_indices += y_start
x_indices += x_start
return binary_volume
# Calculate distances from the center of the blob to voxels within the bounding box
dist = np.sqrt(
(x_indices - x) ** 2 + (y_indices - y) ** 2 + (z_indices - z) ** 2
)
binary_volume[z_start:z_end, y_start:y_end, x_start:x_end][
dist <= radius
] = True
return blobs, binary_volume
......@@ -265,11 +265,13 @@ def minimum(vol, **kwargs):
def tophat(vol, **kwargs):
"""
Remove background from the volume
Remove background from the volume.
Args:
vol: The volume to remove background from
radius: The radius of the structuring element (default: 3)
background: color of the background, 'dark' or 'bright' (default: 'dark'). If 'bright', volume will be inverted.
Returns:
vol: The volume with background removed
"""
......
......@@ -15,7 +15,7 @@ def circles(blobs, vol, alpha=0.5, color="#ff9900", **kwargs):
Args:
blobs (array-like): An array-like object of blobs, where each blob is represented
as a 4-tuple (p, r, c, radius). Usally the result of `qim3d.processing.Blob().detect()`
as a 4-tuple (p, r, c, radius). Usually the result of `qim3d.processing.blob_detection(vol)`
vol (array-like): The 3D volume on which to plot the blobs.
alpha (float, optional): The transparency of the blobs. Defaults to 0.5.
color (str, optional): The color of the blobs. Defaults to "#ff9900".
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment