kaggleのBiohub-cell tracking during development-（https://www.kaggle.com/competitions/biohub-cell-tracking-during-development ）に参加しています。
ゼブラフィッシュの受精卵の細胞分裂過程において、撮影された3Dの時系列画像から、細胞の位置の特定と細胞分裂イベントの特定を行う、というものです。
細胞の位置の特定方法に、LoG(Laplacian of Gaussian)やDoG(Difference of Gaussians)というものがあります。
今回はこれらの技術について簡単に理解するために記事を書いています。
Blob Detectionとは
Blobとは、画像中に存在する、周囲とは明るさや色が異なるまとまりのある領域を指します。
今回扱うのは蛍光顕微鏡画像で、細胞核が暗い背景の中にある明るい塊として観測されます。そのため、細胞核をBlobとして捉え、その中心座標を検出することができます。
Blob Detectionでは、単に明るい部分を探すのではなく、ノイズによる1画素だけの明るさや、画像全体に広がる緩やかな明るさの変化を除外し、特定の大きさを持つまとまりを検出する必要があります。
本記事では、以下の人工二次元画像を使用します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import (
gaussian_filter,
laplace,
maximum_filter,
)
height, width = 128, 128
y, x = np.mgrid[:height, :width]
image = np.zeros((height, width), dtype=np.float64)
cells = [
(35, 40, 5, 1.0),
(80, 55, 8, 0.8),
(70, 95, 12, 0.9),
]
for cy, cx, cell_sigma, intensity in cells:
image += intensity * np.exp(
-(
(x - cx) ** 2
+ (y - cy) ** 2
)
/ (2 * cell_sigma**2)
)
rng = np.random.default_rng(42)
image += rng.normal(
loc=0.0,
scale=0.05,
size=image.shape,
)
plt.figure(figsize=(6, 6))
plt.imshow(image, cmap="gray")
plt.title("Artificial fluorescence image")
plt.axis("off")
plt.show()
以降では、この画像にGaussianフィルタ、LoG、DoG、Local Maxima検出を適用し、それぞれの処理が何をしているのかを確認します。
LoGについて
LoGは大きく3つのstepからなります。
- Gaussianフィルタによる畳み込み
- Laplacian演算&LoG応答
- Local Maximaと細胞中心
順に追っていきます。
1.Gaussianフィルタによる畳み込み
二次元では
G(x,y;\sigma)=\frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}
三次元だと
G(x,y,z;\sigma)=\frac{1}{(2\pi)^{3/2}\sigma^3}e^{-\frac{x^2+y^2+z^2}{2\sigma^2}}
のガウス分布を元にしたkernelによる画像の畳み込みを行うことによって、画像をぼかして、ノイズを軽減させます。
sigmas = [0, 1, 2, 4, 8]
fig, axes = plt.subplots(
1,
len(sigmas),
figsize=(20, 4),
)
for ax, sigma in zip(axes, sigmas):
if sigma == 0:
filtered_image = image
title = "Original"
else:
filtered_image = gaussian_filter(
image,
sigma=sigma,
)
title = rf"$\sigma={sigma}$"
ax.imshow(
filtered_image,
cmap="gray",
vmin=image.min(),
vmax=image.max(),
)
ax.set_title(title)
ax.axis("off")
plt.tight_layout()
plt.show()
2.Laplacian演算&LoG応答
Gaussianフィルタによってノイズを軽減した画像に対して、次に二階微分を行います。
Blobの中心付近では、輝度分布が山形になるため、二階微分は負になります。
sigmas = [1, 2, 4, 8]
fig, axes = plt.subplots(
1,
len(sigmas) + 1,
figsize=(20, 4),
)
axes[0].imshow(image, cmap="gray")
axes[0].set_title("Original")
axes[0].axis("off")
for ax, sigma in zip(axes[1:], sigmas):
smoothed = gaussian_filter(
image,
sigma=sigma,
)
log_response = laplace(smoothed)
max_abs = np.max(np.abs(log_response))
ax.imshow(
log_response,
cmap="seismic",
vmin=-max_abs,
vmax=max_abs,
)
ax.set_title(
rf"LoG response ($\sigma={sigma}$)"
)
ax.axis("off")
plt.tight_layout()
plt.show()
LoG応答を見やすく色付けした画像が以下のものです。
Gaussianフィルタのσの値によって、LoGの応答や細胞の検出結果が変化していることが見て取れます。
3.Local Maximaと細胞中心
前節で確認したように、明るいBlobの中心ではLoGの応答が負になりますので、LoG応答の符号を反転します。
R(x,y)=-\nabla^2(G_\sigma * I)(x,y)
実際に、先ほどの人工画像から細胞中心を検出してみます。
sigmas = [1, 2, 4, 8]
threshold_ratio = 0.2
maximum_filter_size = 9
fig, axes = plt.subplots(
1,
len(sigmas),
figsize=(16, 4),
)
for ax, sigma in zip(axes, sigmas):
smoothed = gaussian_filter(
image,
sigma=sigma,
)
# 明るいBlobの中心が正になるように符号を反転
log_response = -laplace(smoothed)
neighborhood_max = maximum_filter(
log_response,
size=maximum_filter_size,
)
local_maxima = (
log_response == neighborhood_max
)
threshold = (
threshold_ratio
* log_response.max()
)
detected = (
local_maxima
& (log_response > threshold)
)
detected_y, detected_x = np.where(detected)
ax.imshow(
image,
cmap="gray",
vmin=image.min(),
vmax=image.max(),
)
ax.scatter(
detected_x,
detected_y,
s=80,
c="red",
)
ax.set_title(rf"$\sigma={sigma}$")
ax.axis("off")
plt.tight_layout()
plt.show()
さて、σ=4,8ではうまく細胞中心が検出されていますが、特にσ=1の時には細胞中心が検出されていません。
今回は、各σにおけるLoG応答の最大値の20%を閾値として設定しました。
σ=1ではGaussianフィルタによる平滑化が弱いため、背景ノイズによる小さなピークが多く残り、細胞以外の点も検出されています。一方、σを大きくすると細かなノイズが抑えられ、細胞中心に対応するピークが検出されやすくなります。
閾値を高くするとノイズによる誤検出は減りますが、応答の弱い細胞も検出できなくなる可能性があります。このように、Blob Detectionではσと閾値の両方を対象画像に合わせて調整する必要があります。
DoGについて
DoG(Difference of Gaussians)は、異なる2つのGaussianフィルタの差分を利用してLoGを近似する手法です。
LoGではGaussianフィルタを適用した後にLaplacian演算を行いますが、DoGではLaplacianを計算する代わりに、ぼかし具合の異なる2枚の画像の差分を計算します。
DoGは次の3つのステップからなります。
- 異なるσでのGaussianフィルタを適用
- 2枚の画像の差分を計算
- Local Maximaと細胞中心
1.異なる2つのGaussianフィルタを適用
例えば今回はσ=2,4を利用してみます。
sigma_small = 2
sigma_large = 4
gaussian_small = gaussian_filter(
image,
sigma=sigma_small,
)
gaussian_large = gaussian_filter(
image,
sigma=sigma_large,
)
fig, axes = plt.subplots(
1,
2,
figsize=(8, 4),
)
axes[0].imshow(
gaussian_small,
cmap="gray",
vmin=image.min(),
vmax=image.max(),
)
axes[0].set_title(
rf"Gaussian $\sigma={sigma_small}$"
)
axes[0].axis("off")
axes[1].imshow(
gaussian_large,
cmap="gray",
vmin=image.min(),
vmax=image.max(),
)
axes[1].set_title(
rf"Gaussian $\sigma={sigma_large}$"
)
axes[1].axis("off")
plt.tight_layout()
plt.show()
こちらのσ=2とσ=4を利用します。
2.二枚の画像の差分を計算
dog_response = (
gaussian_small
- gaussian_large
)
max_abs = np.max(np.abs(dog_response))
fig, axes = plt.subplots(
1,
3,
figsize=(12, 4),
)
axes[0].imshow(
gaussian_small,
cmap="gray",
vmin=image.min(),
vmax=image.max(),
)
axes[0].set_title(
rf"Gaussian $\sigma={sigma_small}$"
)
axes[0].axis("off")
axes[1].imshow(
gaussian_large,
cmap="gray",
vmin=image.min(),
vmax=image.max(),
)
axes[1].set_title(
rf"Gaussian $\sigma={sigma_large}$"
)
axes[1].axis("off")
dog_plot = axes[2].imshow(
dog_response,
cmap="seismic",
vmin=-max_abs,
vmax=max_abs,
)
axes[2].set_title(
rf"DoG: $\sigma={sigma_small} - \sigma={sigma_large}$"
)
axes[2].axis("off")
fig.colorbar(
dog_plot,
ax=axes[2],
fraction=0.046,
pad=0.04,
)
plt.tight_layout()
plt.show()
σ=2でぼかした画像から、σ=4でぼかした画像を引くことでDoG応答を計算しました。
小さいσでは細胞の局所的な明るさが比較的強く残りますが、大きいσでは明るさがより広範囲に拡散します。この二つの差分を取ることで、画像全体の緩やかな明るさの変化を抑えながら、特定の大きさを持つBlobを強調できます。
Gaussian関数には、σに関する微分とLaplacianの間に、次の関係があります。
\frac{\partial G}{\partial \sigma}
=
\sigma \nabla^2 G
次に、2つのスケールをσとkσとします。kが1に近い場合、Taylor展開によって次のように近似できます。
G(x,y;k\sigma)
\approx
G(x,y;\sigma)
+
(k\sigma-\sigma)
\frac{\partial G}{\partial \sigma}
ここで
k\sigma-\sigma=(k-1)\sigma
なので、
G(x,y;k\sigma)-G(x,y;\sigma)
\approx
(k-1)\sigma
\frac{\partial G}{\partial \sigma}
となります。
先ほどの関係
\frac{\partial G}{\partial \sigma}
=
\sigma\nabla^2G
を代入すると、
G(x,y;k\sigma)-G(x,y;\sigma)
\approx
(k-1)\sigma^2\nabla^2G
が得られます。
したがって、異なる2つのGaussian関数の差分は、定数倍を除けば、スケール正規化されたLoG
\sigma^2\nabla^2G
に近い応答になります。
3.Local Maximaと細胞中心
得られたDoG応答から、各DoG応答の最大値の20%を閾値として細胞中心の検出を行います。
sigma_pairs = [
(1, 2),
(2, 4),
(4, 8),
(1, 4),
(2, 8),
]
threshold_ratio = 0.2
maximum_filter_size = 9
fig, axes = plt.subplots(
1,
len(sigma_pairs),
figsize=(20, 4),
)
for ax, (sigma_small, sigma_large) in zip(
axes,
sigma_pairs,
):
gaussian_small = gaussian_filter(
image,
sigma=sigma_small,
)
gaussian_large = gaussian_filter(
image,
sigma=sigma_large,
)
dog_response = (
gaussian_small
- gaussian_large
)
neighborhood_max = maximum_filter(
dog_response,
size=maximum_filter_size,
)
local_maxima = (
dog_response == neighborhood_max
)
threshold = (
threshold_ratio
* dog_response.max()
)
detected = (
local_maxima
& (dog_response > threshold)
)
detected_y, detected_x = np.where(detected)
ax.imshow(
image,
cmap="gray",
vmin=image.min(),
vmax=image.max(),
)
ax.scatter(
detected_x,
detected_y,
s=45,
c="red",
)
ax.set_title(
rf"$\sigma={sigma_small}, {sigma_large}$"
f"\nDetected: {len(detected_x)}"
)
ax.axis("off")
plt.tight_layout()
plt.show()
今回の画像では、σ=2,4、σ=4,8、σ=2,8の組み合わせで正しく中心検出を行うことができました。
まとめ
本記事では、人工的に作成した蛍光顕微鏡画像を用いて、LoGとDoGによるBlob Detectionの基本的な流れを確認しました。
LoGでは、Gaussianフィルタによって画像のノイズを抑えた後、Laplacianを計算することでBlobを強調します。一方、DoGでは、異なるσでぼかした2枚の画像の差分を取ることで、LoGに近い応答を得ます。
どちらの手法でも、得られた応答のLocal Maximaを探し、一定の閾値を超えた点を細胞中心の候補として検出できます。