Finding The Average Pixel Values Of A List Of Blobs Identified By Scikit-image's Blob_log (laplacian Of Gaussian) Method

Input is a uint16 grayscale .tif-image, 512 x 512 pixels. As the title to this question implies, I would like to calculate the average pixel intensity of blobs identified by the bl

Solution 1:

The log blob detector returns an array in which each row represents a blob. It gives center coordinates (row and column) as well as the Gaussian blog width, sigma. The radius of a blob in 2D is approximately sqrt(2) * sigma.

import numpy as np
from skimage import draw, feature
import matplotlib.pyplot as plt

# Create a test image
image = np.zeros((200, 200))

# Parameters for test image blobs
positions_r, positions_c = (50, 100, 150), (50, 100, 150)
radii = (20, 15, 30)
values = (.5, .75, 1)

# We'll make the blobs vary a bit with noise, to make it# more realistic.  Here we'll store their known average values.
blob_avg = []

for (r, c, radius, val) inzip(positions_r, positions_c,
                               radii, values):

    # Get the coordinates for a circle at the specified coordinates
    r_ix, c_ix = draw.circle(r, c, radius)

    # Generate values for all pixels inside the circle, varying# between val/2 and val.
    noisy_vals = val * 0.5 + np.random.random(size=len(r_ix)) / 2# Put those values into our test image
    image[r_ix, c_ix] = noisy_vals

    # Save the average value of this blob
    blob_avg.append(noisy_vals.mean())

# Compute the blobs in the image, setting the desired sigma range,# and lowering the threshold so that we also grab our faintest blob
blobs_log = feature.blob_log(image, min_sigma=5, max_sigma=50,
                             threshold=0.3, num_sigma=50)

# `blob_log` returns the blobs in reverse order (in this case),# so to compare with our test data we reverse the list of blob# averages
blob_avg = blob_avg[::-1]

# Compute each blob's radius, by multiplying its sigma by sqrt(2)
blobs_log[:, 2] = blobs_log[:, 2] * np.sqrt(2)

# Create a plot, and display our test data
f, ax = plt.subplots(figsize=(15, 10))
ax.imshow(image, cmap='gray');

# Generate all row and column coordinates for our test image# For an `(N, M)` test image, `ixs` will have shape `(N, M, 2)`,# since it stores row and column coordinates.
ixs = np.indices(image.shape)

# Now, we plot each detected blob and estimate its average intensityfor i, blob inenumerate(blobs_log):
    y, x, r = blob
    c = plt.Circle((x, y), r, color='red', linewidth=2, fill=False)
    ax.add_patch(c)

    # Define an array of shape `[2, 1, 1]`, containing# the center of the blob
    blob_center = np.array([y, x])[:, np.newaxis, np.newaxis]

    # Using the formula for a circle, `x**2 + y**2 < r**2`,# generate a mask for this blob.
    mask = ((ixs - blob_center)**2).sum(axis=0) < r**2# Calculate the average intensity of pixels under the mask
    blob_avg_est = image[mask].mean()

    print(f'Blob {i} average value: true={blob_avg[i]:.2f}, estimated={blob_avg_est:.2f}')

Output:

Blob0 average value: true=0.75, estimated=0.75Blob1 average value: true=0.63, estimated=0.63Blob2 average value: true=0.50, estimated=0.49