测量圆形等距区域的径向强度分布

Question

嗨，我目前遇到一个问题，希望有人能够帮助我解决。我尝试四处寻找解决方案，但运气不佳。我是一名生物学家，开始将自动图像分析集成到我的工作流程中，所以如果这不是问这个问题的正确地方，请原谅我。我已经使用 skimage 进行了分割，以识别通过显微镜获得的一组图像（下面的示例）中的对象，如下所示：

然而，我现在想做的是获取这些物体中的每一个，并通过在物体内创建 n 个等间隔的环来测量它们各自的径向轮廓。由此，我想从每个环中提取平均强度和总强度，以量化径向强度的变化（下面的示例）。我遇到的最大问题实际上是找到一种方法，用 skimage 将这些对象重新分割成等距的环。我尝试从对象中心创建多个环，但是由于对象边界不是完美的圆，这并不能达到我们想要的效果。

仅供参考 - 这也是使用 CellProfiler 软件可用的功能，我过去曾使用该软件来完成此操作（可以在此处找到功能示例），但是我正在尝试独立于 CellProfiler 重新创建此功能，以便我们可以运行我们的分析更有效。任何正确方向的指导甚至推动都将不胜感激。

Answer 1

假设我们有对象中心的坐标 (x0, y0) 以及对象边缘的 M 坐标集 H={(hx1, hy1), (hx2, hy2), ..., (hxM ，hyM）}。
在这种情况下，我们的目标是找到对物体进行均等分割的 N 个环或闭合曲线的坐标。所有环都与位于 (x0, y0) 的对象具有相同的中心。
我建议我们可以通过计算 N 个均匀间隔的部分的坐标来实现，将位于对象中心和对象边缘之间的所有线分割为 L={ {(x0, y0), (hx1, hy1)}, {( x0, y0), (hx2, hy2)}, ..., {(x0, y0), (hxM, hyM)} }.
例如，利用这个公式https://www.geeksforgeeks.org/section-formula-point-divides-line-given-ratio/

Answer 2

我认为您分享的工具在幕后与您一样接近完美的圆形或椭圆形，请参见下面的屏幕截图： 测量物体强度分布

但是，如果您想遵循对象边界，最简单的方法可能是通过各种因素重新缩放图像。

这里是如何实现它的。 我将使用 skimage 中的 binary_blobs 来复制您拥有的图像：

from skimage import data
from scipy import ndimage

img = data.binary_blobs(length=100, blob_size_fraction=0.2, volume_fraction=0.2, rng=0).astype(np.uint8) # your binary original image
img_orig = ndimage.distance_transform_edt(img) # your original image

它将生成以下图像： 模拟二值图像 模拟原始图像 所以， img 是您的二值图像，而 img_orig 是您的原始图像（已使用距离变换）。

现在，让我们定义一个函数，它将采用二进制分段对象，并通过以不同的比率重新缩放二进制分段对象来生成具有各种半径的二进制掩码。比率从 0.4 开始，即图像将首先按该因子重新缩放，然后增加直至 1。如果需要，您可以更改其范围。

from skimage import transform, filters, measure
import numpy as np

# you can also use skimage.measure.centroid or skimage.measure.regionprops to get centroid of the binary object
def get_centroid(pattern):
    """
    Return centroid of an image.
    """
    coords = np.array(np.where(pattern)).T # coords = measure.centroid(pattern)
    return np.round(np.mean(coords, axis=0)).astype(int)

def get_masks(obj_mask, bins=4):
    """
    Return object radial masks with its "radial ratios" (fractions of object "radius")
    """
    obj_centroid = get_centroid(obj_mask)
    ratios = np.linspace(0.4, 1, bins)
    masks = []
    for r in ratios:
        rescaled = transform.rescale(obj_mask, r) # resizes image by the given factor

        # after rescaling, we need to make it back binary
        t = filters.threshold_otsu(rescaled)
        rescaled_binary = (rescaled > t).astype(np.uint8)

        # now, since we have smaller image, we need to put the mask back to original image frame
        # this is similar to padding with 0s, but here we need to align small image centroid with original image centroid
        # first we create blank image with 0s
        rescaled_mask = np.zeros_like(obj_mask)

        # then we need to find starting coordinates
        rescaled_centroid = get_centroid(rescaled_binary)
        r_start = obj_centroid[0] - rescaled_centroid[0]
        c_start = obj_centroid[1] - rescaled_centroid[1]
        
        # placing small (rescaled) mask into original image frame
        rescaled_mask[r_start:r_start+rescaled_binary.shape[0], c_start:c_start+rescaled_binary.shape[1]] = rescaled_binary

        masks.append(rescaled_mask)

    return ratios, masks

为了计算径向轮廓，我们将通过几个（箱数）因子重新缩放对象，并计算生成的掩模内的平均强度。

def get_radial_profile(obj, obj_mask, bins=4):
    """
    obj - image of the original object (bboxed)
    obj_mask - binary version of original image (bboxed)
    bins - number of bins

    Return radial profile of the image
    """
    ratios, masks = get_masks(obj_mask, bins=4)
    radial_profiles = [(obj*m).sum()/m.sum() for m in masks] # (obj*m).mean() can be used instead, depending of your needs
    return ratios, radial_profiles

现在，我们应该迭代所有分段对象，找到边界框，使用该边界框内的图像，然后计算其径向轮廓。

import matplotlib.pyplot as plt

img_segmented = measure.label(img)

for seg_num in np.unique(img_segmented)[1:]:
    obj = (img_segmented==seg_num).astype(np.uint8)
    props = measure.regionprops(obj)
    r_start, c_start, r_end, c_end = props[0].bbox
    obj_bbox = img_orig[r_start:r_end, c_start:c_end]
    obj_bbox_mask = obj[r_start:r_end, c_start:c_end]

    ratios, radial_profile = get_radial_profile(obj_bbox, obj_bbox_mask)
    _, axes = plt.subplots(1, 2, figsize=(10,5))
    axes[0].imshow(obj_bbox)
    axes[1].plot(ratios, radial_profile)
    axes[1].grid()

综合起来：

from skimage import data
import numpy as np
from skimage import transform, filters, measure
from scipy import ndimage
import matplotlib.pyplot as plt


img = data.binary_blobs(length=100, blob_size_fraction=0.2, volume_fraction=0.2, rng=0).astype(np.uint8) # your binary original image
img_orig = ndimage.distance_transform_edt(img) # your original image

# you can also use skimage.measure.centroid or skimage.measure.regionprops to get centroid of the binary object
def get_centroid(pattern):
    """
    Return centroid of an image.
    """
    coords = np.array(np.where(pattern)).T # coords = measure.centroid(pattern)
    return np.round(np.mean(coords, axis=0)).astype(int)

def get_masks(obj_mask, bins=4):
    """
    Return object radial masks with its "radial ratios" (fractions of object "radius")
    """
    obj_centroid = get_centroid(obj_mask)
    ratios = np.linspace(0.4, 1, bins)
    masks = []
    for r in ratios:
        rescaled = transform.rescale(obj_mask, r) # resizes image by the given factor

        # after rescaling, we need to make it back binary
        t = filters.threshold_otsu(rescaled)
        rescaled_binary = (rescaled > t).astype(np.uint8)

        # now, since we have smaller image, we need to put the mask back to original image frame
        # this is similar to padding with 0s, but here we need to align small image centroid with original image centroid
        # first we create blank image with 0s
        rescaled_mask = np.zeros_like(obj_mask)

        # then we need to find starting coordinates
        rescaled_centroid = get_centroid(rescaled_binary)
        r_start = obj_centroid[0] - rescaled_centroid[0]
        c_start = obj_centroid[1] - rescaled_centroid[1]

        # placing small (rescaled) mask into original image frame
        rescaled_mask[r_start:r_start+rescaled_binary.shape[0], c_start:c_start+rescaled_binary.shape[1]] = rescaled_binary

        masks.append(rescaled_mask)

    return ratios, masks

def get_radial_profile(obj, obj_mask, bins=4):
    """
    obj - image of the original object (bboxed)
    obj_mask - binary version of original image (bboxed)
    bins - number of bins

    Return radial profile of the image
    """
    ratios, masks = get_masks(obj_mask, bins=4)
    radial_profiles = [(obj*m).sum()/m.sum() for m in masks] # (obj*m).mean() can be used instead, depending of your needs
    return ratios, radial_profiles

img_segmented = measure.label(img)

for seg_num in np.unique(img_segmented)[1:]:
    obj = (img_segmented==seg_num).astype(np.uint8)
    props = measure.regionprops(obj)
    r_start, c_start, r_end, c_end = props[0].bbox
    obj_bbox = img_orig[r_start:r_end, c_start:c_end]
    obj_bbox_mask = obj[r_start:r_end, c_start:c_end]

    ratios, radial_profile = get_radial_profile(obj_bbox, obj_bbox_mask)
    _, axes = plt.subplots(1, 2, figsize=(10,5))
    axes[0].imshow(obj_bbox)
    axes[1].plot(ratios, radial_profile)
    axes[1].grid()

其中一个对象的结果： 对象径向轮廓

Alternatevely：如果您仍然想不近似对象边界的形状，而是近似椭圆，您可以使用以下 get_masks 函数代替：

from skimage import draw

def get_masks(obj_mask, bins=4):
    """
    Return object radial masks with its "radial ratios" (fractions of object "radius")
    """
    props = measure.regionprops(obj_mask)
    centroid = props[0].centroid
    orientation = props[0].orientation
    ratios = np.linspace(0.4, 1, bins)
    masks = []
    for r in ratios:
        rr, cc = draw.ellipse(centroid[0], centroid[1], r*obj_mask.shape[0]/2, r*obj_mask.shape[1]/2, shape=obj_mask.shape, rotation=orientation)
        blank = np.zeros_like(obj_bin_bbox)
        blank[rr, cc] = 1
        masks.append(blank)
    return ratios, masks

借助

skimage.draw

功能，我们正在创建具有各种比例的椭圆蒙版。

请也看看这个帖子： 计算径向轮廓的最有效方法 您可以在分段对象上使用here定义的radial_profile函数。它将为您提供“圆形”径向轮廓。但稍加修改，您可以获得“椭圆形”径向轮廓。像这样的东西：

import numpy as np

def elliptical_radial_profile(data, center):
    y, x = np.indices((data.shape))
    r = np.sqrt(data.shape[0]*data.shape[1])*np.sqrt((x - center[0])**2 + (y - center[1])**2)
    r = r.astype(np.int)

    tbin = np.bincount(r.ravel(), data.ravel())
    nr = np.bincount(r.ravel())
    radialprofile = tbin / nr
    return radialprofile 

注意： 这个

elliptical_radial_profile

函数与我在

get_masks

中使用的函数之间的区别在于，

elliptical_radial_profile

不会考虑对象方向，但在

get_masks

中会考虑对象方向。但另一方面

elliptical_radial_profile

会给你更多的粒度。您可能需要细化椭圆方程以添加对象的方向。