trackpy-particle-tracking

trackpy

概要

trackpy は、ビデオ顕微鏡における単一粒子追跡（SPT）用の Python ライブラリです。Crocker-Grier アルゴリズムを実装しており、各フレームで明るいスポットをサブピクセル精度で検出し、フレーム間で位置を連結して連続した軌跡を生成します。軌跡から、trackpy は平均二乗変位（MSD）、拡散係数、運動の分類（閉じ込められた運動、通常の拡散、方向性のある運動）を計算します。2D 蛍光動画、3D 共焦点 z スタック、大規模な画像シーケンスに対応しており、pims 画像読み込みライブラリによるメモリ効率の高いストリーミング処理が可能です。

使用する場面

• 標識した粒子（量子ドット、蛍光ビーズ、小胞、受容体）の蛍光顕微鏡動画があり、個々の軌跡と拡散係数を抽出する必要がある場合。
• 粒子の移動性を測定し、単一粒子の軌跡から MSD 曲線を計算して、ブラウン拡散、方向性のある運動、閉じ込められた運動を区別したい場合。
• 粒子ごとの軌跡データが必要なコロイドの動態、脂質膜の拡散、細胞内の物質輸送、ウイルスと細胞の相互作用を解析する場合。
• 共焦点 z スタックの時系列から 3D 追跡を行い、粒子や細胞小器官の面外運動を捉える必要がある場合。
• 固有の粒子運動統計量を計算する前に、ステージドリフトを除去するためのドリフト補正を適用したい場合。
• 数百の軌跡で平均したアンサンブル MSD が必要で、十分な統計的検出力をもって集団レベルの拡散挙動を抽出したい場合。
• グラフィカルインターフェース、軌跡の手動キュレーション、生物学的オブジェクトのセグメンター（Cellpose、StarDist）との統合が必要な場合は、代わりに TrackMate（Fiji/ImageJ プラグイン）を使用してください。
• 画像データと並べて軌跡をインタラクティブに可視化し、手動で編集したい場合は、代わりに napari と napari-trackpy を使用してください。

前提条件

• Python パッケージ：trackpy、pims、pandas、numpy、matplotlib、scipy
• データ要件：グレースケールまたは単一チャンネルの画像シーケンス（TIF スタック、AVI、または PNG/TIF フレームのディレクトリ）。粒子は暗い背景上の明るいガウス状スポットとして表示される必要があります（明るい背景上の暗いスポットには invert=True を使用してください）。
• 環境：Jupyter ノートブックとスクリプトで動作します。pims はほとんどの顕微鏡画像形式に対応しています。ND2 または CZI ファイルには、pims-nd2 または aicsimageio をインストールしてください。

pip install trackpy pims pandas numpy matplotlib scipy
# For reading multi-channel or proprietary formats:
pip install pims[bioformats]   # Bioformats via JPype
pip install aicsimageio        # ND2, CZI, LIF via AICSImageIO

クイックスタート

import trackpy as tp
import pims

# Load a TIF image stack (T frames × Y × X)
frames = pims.open("particles.tif")   # shape: (T, Y, X)

# Locate particles in all frames
f = tp.batch(frames, diameter=11, minmass=500)
print(f"Found {len(f)} particle detections across {f['frame'].nunique()} frames")

# Link into trajectories
t = tp.link(f, search_range=5, memory=3)

# Remove short-lived tracks (fewer than 10 frames)
t = tp.filter_stubs(t, threshold=10)
print(f"Retained {t['particle'].nunique()} trajectories")

# Compute ensemble MSD
imsd = tp.imsd(t, mpp=0.16, fps=10)   # mpp: microns per pixel, fps: frames per second
print(imsd.head())

コア API

モジュール 1：tp.locate() — 単一フレームの粒子検出

tp.locate() は、バンドパスフィルターを適用してから局所最大値を検出し、1 つの画像フレーム内の明るい円形の特徴を見つけます。検出した各粒子について、サブピクセル精度の x/y 位置、積分質量、信号、離心率を含む DataFrame を返します。

import trackpy as tp
import pims
import matplotlib.pyplot as plt

frames = pims.open("particles.tif")
frame0 = frames[0]   # single 2D array

# Locate particles: diameter must be odd integer, roughly matching spot size in pixels
f0 = tp.locate(frame0, diameter=11, minmass=300, maxsize=None, separation=None)
print(f"Detected {len(f0)} particles in frame 0")
print(f0[['x', 'y', 'mass', 'size', 'ecc']].head())
# x, y: subpixel centroid; mass: integrated brightness; size: Gaussian width; ecc: eccentricity (0=circular)

# Diagnostic plot: annotate detected particles on the raw frame
fig, ax = plt.subplots(figsize=(8, 8))
tp.annotate(f0, frame0, ax=ax, imshow_style={"cmap": "gray"})
ax.set_title(f"Frame 0: {len(f0)} particles detected")
plt.tight_layout()
plt.savefig("locate_diagnostic.png", dpi=150)
print("Saved locate_diagnostic.png")

モジュール 2：tp.batch() — 複数フレームの検出

tp.batch() は、画像シーケンス内のすべてのフレームに tp.locate() を適用し、frame 列を含む単一の DataFrame に結果を連結します。pims 対応の画像読み込み機能、または 2D 配列のリストを受け付けます。

import trackpy as tp
import pims

frames = pims.open("particles.tif")

# Locate particles across all frames (same parameters as tp.locate)
f = tp.batch(frames, diameter=11, minmass=300, processes=1)
# processes=1 uses serial processing; set processes="auto" for multicore (requires joblib)
print(f"Total detections: {len(f)}")
print(f"Frames with data: {f['frame'].nunique()} / {len(frames)}")
print(f"Mean particles per frame: {len(f)/f['frame'].nunique():.1f}")
print(f.groupby('frame').size().describe())

# Mass histogram: use to choose minmass cutoff
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(6, 4))
f['mass'].hist(bins=40, ax=ax)
ax.axvline(300, color='red', linestyle='--', label='minmass=300')
ax.set_xlabel("Integrated mass")
ax.set_ylabel("Count")
ax.set_title("Mass distribution of detections")
ax.legend()
plt.tight_layout()
plt.savefig("mass_histogram.png", dpi=150)
print("Saved mass_histogram.png — use to refine minmass cutoff")

モジュール 3：tp.link() — 軌跡の連結

tp.link() は、二部割当問題（ハンガリー法）を解くことで、フレーム間の粒子検出結果を連結して軌跡を生成します。位置 DataFrame に particle 列（整数の軌跡 ID）を追加します。search_range（ピクセル）はフレーム間の最大変位です。memory を指定すると、粒子が除外されるまでに最大 N フレームのあいだ見失われても許容されます。

import trackpy as tp
import pims

frames = pims.open("particles.tif")
f = tp.batch(frames, diameter=11, minmass=300)

# Link: search_range in pixels; memory handles brief disappearances (blinking, out-of-focus)
t = tp.link(f, search_range=5, memory=3)
print(f"Number of unique trajectories: {t['particle'].nunique()}")
print(f"Trajectory length distribution:")
print(t.groupby('particle').size().describe())

# Visualize all trajectories overlaid on the first frame
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8, 8))
tp.plot_traj(t, superimpose=frames[0], ax=ax)
ax.set_title(f"{t['particle'].nunique()} trajectories")
plt.tight_layout()
plt.savefig("trajectories.png", dpi=150)
print("Saved trajectories.png")

モジュール 4：tp.filter_stubs() — 短い軌跡の除去

tp.filter_stubs() は、指定したフレーム数より短い軌跡を除去します。短い軌跡は、ノイズの検出、視野への粒子の進入や視野からの退出、または連結エラーによって生じます。短い軌跡を除去すると、長いラグ時間では短い軌跡が分散の大きい MSD 推定値をもたらすため、MSD の信頼性が向上します。

import trackpy as tp
import pims

frames = pims.open("particles.tif")
f = tp.batch(frames, diameter=11, minmass=300)
t = tp.link(f, search_range=5, memory=3)

before = t['particle'].nunique()
t_filt = tp.filter_stubs(t, threshold=10)   # keep only tracks with ≥10 frames
after = t_filt['particle'].nunique()
print(f"Tracks before filtering: {before}")
print(f"Tracks after filtering (≥10 frames): {after}")
print(f"Removed {before - after} short tracks ({100*(before-after)/before:.1f}%)")

モジュール 5：MSD 解析 — tp.imsd() と tp.emsd()

tp.imsd() は、ラグ時間に応じた粒子ごとの平均二乗変位を計算し、DataFrame（インデックスはラグ時間、列は粒子 ID）を返します。tp.emsd() は、すべての粒子のアンサンブル平均 MSD を計算します。どちらも、物理スケール（mpp、ピクセルあたりのマイクロメートル数）とフレームレート（fps）が必要です。

import trackpy as tp
import pims
import matplotlib.pyplot as plt

frames = pims.open("particles.tif")
f = tp.batch(frames, diameter=11, minmass=300)
t = tp.link(f, search_range=5, memory=3)
t = tp.filter_stubs(t, threshold=10)

mpp = 0.16    # microns per pixel (from microscope calibration)
fps = 10.0    # frames per second

# Individual MSD curves (one column per particle)
imsd = tp.imsd(t, mpp=mpp, fps=fps, max_lagtime=100)
print(f"IMSD shape: {imsd.shape}")  # (lag times) × (particles)

# Ensemble MSD
emsd = tp.emsd(t, mpp=mpp, fps=fps, max_lagtime=100)
print(f"EMSD at lag 1 s: {emsd.iloc[0]:.4f} µm²")

# Plot ensemble MSD and fit diffusion coefficient
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import linregress

mpp = 0.16
fps = 10.0

# Fit MSD = 4*D*t (2D Brownian) over first 10 lag times
lag_s = emsd.index.values[:10]   # lag times in seconds
msd_vals = emsd.values[:10]
slope, intercept, r, p, se = linregress(lag_s, msd_vals)
D = slope / 4   # diffusion coefficient in µm²/s
print(f"Diffusion coefficient D = {D:.4f} µm²/s  (R²={r**2:.3f})")

fig, ax = plt.subplots(figsize=(6, 5))
ax.plot(emsd.index, emsd.values, 'o-', label='Ensemble MSD')
ax.plot(lag_s, slope * lag_s + intercept, 'r--', label=f'Fit: D={D:.4f} µm²/s')
ax.set_xlabel("Lag time (s)")
ax.set_ylabel("MSD (µm²)")
ax.set_title("Ensemble Mean Squared Displacement")
ax.legend()
plt.tight_layout()
plt.savefig("emsd.png", dpi=150)
print("Saved emsd.png")

モジュール 6：運動解析 — 特性評価とドリフト補正

tp.motion.characterize() は、軌跡ごとの統計量（平均速度、正味変位、直線性）を計算します。tp.subtract_drift() は、MSD 解析の前に軌跡からステージ全体のドリフトを除去します。

import trackpy as tp
import pims

frames = pims.open("particles.tif")
f = tp.batch(frames, diameter=11, minmass=300)
t = tp.link(f, search_range=5, memory=3)
t = tp.filter_stubs(t, threshold=10)

# Estimate and subtract drift (bulk movement of the sample/stage)
drift = tp.compute_drift(t)
print("Drift (first 5 frames):")
print(drift.head())

t_corrected = tp.subtract_drift(t.copy(), drift)
print(f"Drift subtracted from {t_corrected['particle'].nunique()} trajectories")

import trackpy as tp

# Characterize individual trajectories (requires tp.motion module)
from trackpy import motion

# Per-particle summary statistics
char = motion.characterize(t, mpp=0.16, fps=10.0)
print(char.columns.tolist())
# Columns: 'alpha' (anomalous exponent), 'D_app' (apparent diffusion), 'r^2' (fit quality)
print(char[['alpha', 'D_app']].describe())
# alpha ~ 1.0: Brownian; alpha < 1: confined/subdiffusion; alpha > 1: directed/superdiffusion

よくあるワークフロー

ワークフロー 1：MSD と拡散係数を求める 2D トラッキングの全工程

目標：蛍光動画を読み込み、すべてのフレームで粒子の位置を特定して連結し、短い軌跡を除外して、MSD を計算し、拡散係数を求めます。

import trackpy as tp
import pims
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import linregress

# ── 1. Load image sequence ──────────────────────────────────────────────────
frames = pims.open("fluorescence_video.tif")   # (T, Y, X) grayscale TIF stack
print(f"Loaded {len(frames)} frames, frame shape: {frames.frame_shape}")

# ── 2. Tune detection on a single frame ─────────────────────────────────────
f0 = tp.locate(frames[0], diameter=11, minmass=400)
print(f"Frame 0: {len(f0)} particles detected")
# Adjust diameter (odd integer ≥ spot size) and minmass until count looks right

# ── 3. Batch detect across all frames ───────────────────────────────────────
f = tp.batch(frames, diameter=11, minmass=400, processes=1)
print(f"Total detections: {len(f)} across {f['frame'].nunique()} frames")

# ── 4. Link into trajectories ────────────────────────────────────────────────
t = tp.link(f, search_range=6, memory=3)
print(f"Unique trajectories before filtering: {t['particle'].nunique()}")

# ── 5. Remove short trajectories ─────────────────────────────────────────────
t = tp.filter_stubs(t, threshold=15)
print(f"Trajectories after filtering (≥15 frames): {t['particle'].nunique()}")

# ── 6. Subtract stage drift ───────────────────────────────────────────────────
drift = tp.compute_drift(t)
t = tp.subtract_drift(t.copy(), drift)

# ── 7. Compute MSD ────────────────────────────────────────────────────────────
mpp = 0.16   # µm/pixel — from microscope calibration
fps = 10.0   # frames per second

emsd = tp.emsd(t, mpp=mpp, fps=fps, max_lagtime=50)
imsd = tp.imsd(t, mpp=mpp, fps=fps, max_lagtime=50)

# ── 8. Fit diffusion coefficient from linear regime (first 10 points) ─────────
n_fit = 10
lag_s = emsd.index.values[:n_fit]
msd_v = emsd.values[:n_fit]
slope, intercept, r, _, _ = linregress(lag_s, msd_v)
D = slope / 4   # MSD = 4Dt for 2D Brownian
print(f"Diffusion coefficient D = {D:.4f} µm²/s  (R²={r**2:.3f})")

# ── 9. Plot ────────────────────────────────────────────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Individual + ensemble MSD
axes[0].plot(imsd.index, imsd.values, alpha=0.2, color='steelblue', linewidth=0.8)
axes[0].plot(emsd.index, emsd.values, 'k-', linewidth=2, label='Ensemble MSD')
axes[0].plot(lag_s, slope * lag_s + intercept, 'r--', label=f'D={D:.4f} µm²/s')
axes[0].set_xlabel("Lag time (s)")
axes[0].set_ylabel("MSD (µm²)")
axes[0].set_title("MSD: individual (blue) + ensemble (black)")
axes[0].legend()

# Trajectory overlay
tp.plot_traj(t, superimpose=frames[0], ax=axes[1])
axes[1].set_title(f"{t['particle'].nunique()} trajectories")

plt.tight_layout()
plt.savefig("tracking_results.png", dpi=150, bbox_inches="tight")
print("Saved tracking_results.png")

# ── 10. Export trajectories ───────────────────────────────────────────────────
t.to_csv("trajectories.csv", index=False)
emsd.to_csv("ensemble_msd.csv")
print("Exported trajectories.csv and ensemble_msd.csv")

ワークフロー 2：共焦点 Z スタックからの 3D 粒子トラッキング

目標：時系列の共焦点 Z スタック（T × Z × Y × X）から粒子を 3D でトラッキングし、3D で連結して、3D MSD を計算します。

import trackpy as tp
import pims
import numpy as np
import matplotlib.pyplot as plt

# ── 1. Load 4D image stack (T × Z × Y × X) ─────────────────────────────────
# pims opens multi-page TIF; reshape into (T, Z, Y, X) as needed
raw = pims.open("confocal_3d_timeseries.tif")
# Assume each "frame" in pims is one Z-slice; reshape to (T, Z, Y, X)
T, n_z = 50, 20   # adjust to match acquisition
frames_4d = np.array(raw).reshape(T, n_z, raw.frame_shape[0], raw.frame_shape[1])
print(f"4D stack shape: {frames_4d.shape}")  # (T, Z, Y, X)

# ── 2. Detect in 3D (tp.locate works on 3D arrays) ─────────────────────────
# For 3D, pass a single 3D volume; diameter can be (z_diam, y_diam, x_diam)
f0_3d = tp.locate(frames_4d[0], diameter=(7, 11, 11), minmass=2000)
print(f"3D detections in t=0: {len(f0_3d)}")
print(f0_3d[['x', 'y', 'z', 'mass']].head())

# ── 3. Batch detect across all time points ──────────────────────────────────
detections = []
for t_idx in range(T):
    frame_3d = frames_4d[t_idx]    # Z × Y × X
    detected = tp.locate(frame_3d, diameter=(7, 11, 11), minmass=2000)
    detected['frame'] = t_idx
    detections.append(detected)

import pandas as pd
f3d = pd.concat(detections, ignore_index=True)
print(f"Total 3D detections: {len(f3d)}")

# ── 4. Link in 3D ────────────────────────────────────────────────────────────
# search_range in pixels; use a 3-tuple (z, y, x) for anisotropic voxels
t3d = tp.link(f3d, search_range=(3, 6, 6), memory=2)
t3d = tp.filter_stubs(t3d, threshold=10)
print(f"3D trajectories: {t3d['particle'].nunique()}")

# ── 5. Compute 3D MSD ────────────────────────────────────────────────────────
mpp_xy = 0.16   # µm/pixel in x, y
mpp_z = 0.30    # µm/pixel in z (z-step size)
fps = 1.0       # z-stack volume rate

# Scale z coordinates to µm
t3d_um = t3d.copy()
t3d_um['x'] *= mpp_xy
t3d_um['y'] *= mpp_xy
t3d_um['z'] *= mpp_z

emsd_3d = tp.emsd(t3d_um, mpp=1.0, fps=fps, max_lagtime=20)  # mpp=1 since already in µm
print(f"3D ensemble MSD (lag=1 s): {emsd_3d.iloc[0]:.4f} µm²")

# ── 6. Fit 3D diffusion coefficient (MSD = 6Dt for 3D) ─────────────────────
from scipy.stats import linregress
lag_s = emsd_3d.index.values[:8]
slope, _, r, _, _ = linregress(lag_s, emsd_3d.values[:8])
D_3d = slope / 6
print(f"3D Diffusion coefficient D = {D_3d:.4f} µm²/s  (R²={r**2:.3f})")

t3d.to_csv("trajectories_3d.csv", index=False)
print("Saved trajectories_3d.csv")

主なパラメーター

パラメーター	モジュール	デフォルト	範囲／選択肢	効果
diameter	locate, batch	必須	3 以上の奇数（または 3D 用のタプル）	ピクセル単位での粒子直径の概算値です。奇数である必要があります。小さすぎると検出が分割され、大きすぎると検出が結合されます
minmass	locate, batch	100	0～∞	積分輝度の最小値です。ノイズを除外するための主なフィルターです。まず 0 に設定して質量のヒストグラムをプロットし、ノイズのピークを分離する値を設定します
search_range	link	必須	1～50 ピクセル	フレーム間の最大変位です。フレームごとの最大予想移動量の約 1.5 倍に設定します
memory	link	0	0～10 フレーム	粒子が検出されなくても軌跡が途切れないフレーム数です。点滅する蛍光色素に有用です
threshold	filter_stubs	1	1 以上の整数	軌跡の最小長をフレーム数で指定します。短い軌跡では MSD の信頼性が低くなります
max_lagtime	imsd, emsd	100	整数	MSD 計算の最大ラグ時間をフレーム数で指定します。信頼性を高めるには、全フレーム数の約 10～20％を使用します
mpp	imsd, emsd	1	0 より大きい浮動小数点数	ピクセルあたりのマイクロメートル数です。ピクセル単位を物理単位（µm）に変換します
fps	imsd, emsd	1	0 より大きい浮動小数点数	1 秒あたりのフレーム数です。フレーム単位のラグを秒に変換します
separation	locate, batch	diameter+1	整数	検出対象間の中心間距離の最小値です。密集した粒子の二重カウントを防ぎます
invert	locate, batch	False	True, False	明るい背景上の暗い粒子（透過光イメージング）では True に設定します

ベストプラクティス

1. まず1フレームでdiameterとminmassを必ず調整します：代表的なフレーム1枚に対してtp.locate()を実行し、tp.batch()を実行する前にtp.annotate()で検出結果を目視確認します。検出しすぎると時間が無駄になり、検出が足りないと粒子を見逃します。

   f0 = tp.locate(frames[0], diameter=11, minmass=200)
   tp.annotate(f0, frames[0])   # visual check in Jupyter
   

2. search_rangeは控えめに設定します：探索範囲が大きすぎると、密なサンプルで無関係な粒子同士が誤ってリンクされます。リンクする前に、tp.locate()の出力のばらつきから、フレームごとの典型的な変位を推定します。

3. MSDを計算する前にドリフトを差し引きます：ステージのドリフトはMSDを増加させ、Dを過大評価させます。tp.emsd()の前に必ずtp.compute_drift()とtp.subtract_drift()を呼び出します。

4. 拡散係数のフィッティングには線形領域のみを使用します：長いラグタイムでは、MSD曲線がノイズを含むようになります（寄与する軌跡ペアが少なくなります）。利用可能なラグタイムの最初の10～20％のみをフィッティングします。フィッティングの前に対数-対数プロットの傾きを使って、非ブラウン運動を検出します。

5. locateとMSDの処理でmppの単位を混在させないでください：tp.locate()は位置をピクセル単位で返します。mppはtp.imsd()／tp.emsd()でのみ適用されます。リンク処理の前に位置を手動で再スケーリングすると、ピクセル単位のsearch_rangeが機能しなくなるため、避けてください。

6. 異方性ボクセルで3D追跡を行う場合：diameterとsearch_rangeは、(z, y, x)の軸順に対応するタプルとして渡します。通常、z方向のステップはxy方向のピクセルサイズより2～5倍粗いため、それに応じてdiameterとsearch_rangeのz成分を設定します。

よく使うレシピ

レシピ：ドリフト補正と補正後MSDの比較

使用する場面：ステージのドリフトの寄与を評価するため、生のMSDとドリフト補正後のMSDを比較します。

import trackpy as tp
import pims
import matplotlib.pyplot as plt

frames = pims.open("particles.tif")
f = tp.batch(frames, diameter=11, minmass=400, processes=1)
t = tp.link(f, search_range=6, memory=3)
t = tp.filter_stubs(t, threshold=15)

mpp, fps = 0.16, 10.0

# MSD without drift correction
emsd_raw = tp.emsd(t, mpp=mpp, fps=fps, max_lagtime=50)

# Subtract drift
drift = tp.compute_drift(t)
t_corr = tp.subtract_drift(t.copy(), drift)
emsd_corr = tp.emsd(t_corr, mpp=mpp, fps=fps, max_lagtime=50)

fig, ax = plt.subplots(figsize=(6, 5))
ax.loglog(emsd_raw.index, emsd_raw.values, 'r--', label='Raw MSD')
ax.loglog(emsd_corr.index, emsd_corr.values, 'b-', label='Drift-corrected MSD')
ax.set_xlabel("Lag time (s)")
ax.set_ylabel("MSD (µm²)")
ax.set_title("Effect of drift correction on MSD")
ax.legend()
plt.tight_layout()
plt.savefig("drift_correction_comparison.png", dpi=150)
print("Saved drift_correction_comparison.png")

レシピ：拡散レジームによる粒子の分類

使用する場面：対数-対数MSDの傾き（異常指数alpha）に基づいて、粒子集団を閉じ込め、正常（ブラウン）運動、方向性運動に分類します。

import trackpy as tp
import pims
import numpy as np
import pandas as pd
from scipy.stats import linregress

frames = pims.open("particles.tif")
f = tp.batch(frames, diameter=11, minmass=400, processes=1)
t = tp.link(f, search_range=6, memory=3)
t = tp.filter_stubs(t, threshold=20)

drift = tp.compute_drift(t)
t = tp.subtract_drift(t.copy(), drift)

mpp, fps = 0.16, 10.0
imsd = tp.imsd(t, mpp=mpp, fps=fps, max_lagtime=30)

# Fit log-log slope (anomalous exponent alpha) for each particle
results = []
for pid in imsd.columns:
    curve = imsd[pid].dropna()
    if len(curve) < 5:
        continue
    log_lag = np.log(curve.index.values)
    log_msd = np.log(curve.values)
    slope, intercept, r, _, _ = linregress(log_lag[:10], log_msd[:10])
    D_app = np.exp(intercept) / 4   # apparent D from intercept
    results.append({'particle': pid, 'alpha': slope, 'D_app': D_app, 'r2': r**2})

df_char = pd.DataFrame(results)

# Classify by alpha
df_char['regime'] = pd.cut(
    df_char['alpha'],
    bins=[-np.inf, 0.7, 1.3, np.inf],
    labels=['confined', 'brownian', 'directed']
)
print(df_char['regime'].value_counts())
print(f"\nMean D by regime:\n{df_char.groupby('regime')['D_app'].mean()}")

df_char.to_csv("particle_classification.csv", index=False)
print("Saved particle_classification.csv")

レシピ：離心率によるフィルタリングで凝集体を除去

使用する場面：質量のしきい値は満たしているものの、細長い非円形の検出対象（二重粒子、凝集体、デブリ）を除外します。

import trackpy as tp
import pims

frames = pims.open("particles.tif")
f = tp.batch(frames, diameter=11, minmass=400, processes=1)

# Eccentricity: 0 = perfect circle, 1 = line
# Remove elongated features (likely aggregates or debris)
f_round = f[f['ecc'] < 0.3]
print(f"Before ecc filter: {len(f)} detections")
print(f"After ecc filter (ecc<0.3): {len(f_round)} detections")
print(f"Removed: {len(f)-len(f_round)} elongated features")

t = tp.link(f_round, search_range=6, memory=3)
t = tp.filter_stubs(t, threshold=10)
print(f"Trajectories after eccentricity filtering: {t['particle'].nunique()}")

トラブルシューティング

問題	原因	解決策
誤検出が多すぎる	minmassが低すぎる、またはdiameterがスポットサイズに合っていない	質量ヒストグラムをプロットし、ノイズのピークとシグナルのピークの間の谷までminmassを上げます。diameterが実際のスポット幅（ピクセル単位）に合っていることを確認します
検出が少ない、または検出されない	minmassが高すぎる、または粒子が暗い／焦点が合っていない	minmassを下げ、画像のコントラストを確認し、locateの前に背景差し引きを適用します
軌跡が非常に短い（すべての短い軌跡がフィルタリングされる）	粒子の速度に対してsearch_rangeが小さすぎる、または明滅があるのにmemory=0になっている	search_rangeをフレームごとの最大変位の1.5～2倍に増やします。明滅する色素ではmemory=2または3に設定します
長いラグタイムでMSD曲線がノイズを含む、または単調でない	軌跡が少なすぎる、またはフィッティングするラグポイントが多すぎる	フィッティングにはラグタイムの最初の10～20％のみを使用します。アンサンブルMSDには少なくとも50本以上の軌跡があることを確認します
ドリフト補正でMSDが悪化する	動かない基準粒子が少なすぎる、または動いている粒子からドリフトが推定されている	基準マーカーのビーズまたは動かない粒子を含めます。動かないことが分かっている粒子に対してのみtp.compute_drift()を使用します
tp.batch()中にMemoryErrorが発生する	すべてのフレームが一度にRAMに読み込まれている	pimsの遅延読み込みリーダー（デフォルト）を使用します。processes=1に設定し、tp.locate()をループで実行してフレームを分割処理します
3Dのlocateが2D位置のみを返す	3Dボリュームではなく2Dフレームを渡している	入力配列が3次元(Z, Y, X)であることを確認します。frames_4d[t_idx].ndim == 3であることを確認します
リンクされた軌跡が多数の短いセグメントに分断される	粒子がフレーム間でsearch_rangeを超える速度で動いている	search_rangeを増やし、memoryを増やします。フレームレートが非常に高い場合は、フレームの間引きを検討します

参考資料

• trackpy のドキュメント — 公式 API リファレンス、チュートリアル、ノートブック
• trackpy の GitHub リポジトリ — ソースコード、課題トラッカー、サンプルノートブック
• Crocker, J.C. & Grier, D.G. (1996). コロイド研究のためのデジタルビデオ顕微鏡法. J. Colloid Interface Sci. 179, 298–310 — オリジナルのアルゴリズム論文
• pims のドキュメント — trackpy と連携する画像シーケンス読み込みライブラリ
• trackpy のウォークスルーノートブック — 2D トラッキングのステップごとのチュートリアル