napariでPointを配置した輝点を２次元ガウシアンでフィッティングする

やりたいこと

輝点画像にPointを配置して、その輝点を２次元ガウシアンでフィッティングして、輝点の座標を出力したい。

コード

fit_signal_with_gaussian.py

import glob
import math
import napari
import numpy as np

from scipy.optimize import curve_fit
from skimage import io


dir_name = "./data"
file_names = glob.glob(dir_name + "/*.tif")
stack = np.array([io.imread(file_name) for file_name in file_names])

viewer = napari.view_image(stack)
points_layer = viewer.add_points(ndim=3)

ws = 9
def gaussian_2d(x, a, b, c, d):
    return a * np.exp(-(b * ((x % ws - c) ** 2 + ((x / ws).astype(int) - d) ** 2)))

@viewer.bind_key('f')
def fit_gaussian(viewer):
    for point in points_layer.data:
        frame = round(point[0])
        xc = round(point[2])
        yc = round(point[1])

        x0 = xc - int(ws / 2)
        y0 = yc - int(ws / 2)
        
        # Get signals around the specified point
        x = np.array([i for i in range(ws * ws)])
        s = np.array([stack[frame, y0 + int(i / ws), x0 + i % ws] for i in range(ws * ws)])

        # Estimate reciprocal variance for the fitting initial parameter
        v_sum = 0.0
        s_sum = 0.0
        for i in range(ws * ws):
            x_tmp = x0 + i % ws
            y_tmp = y0 + i / ws
            s_tmp = stack[frame, y0 + int(i / ws), x0 + i % ws]
            v_sum += s_tmp * ((x_tmp - xc) ** 2 + (y_tmp - yc) ** 2)
            s_sum += s_tmp 
        rv = s_sum / v_sum
        
        # Fitting
        popt, pcov = curve_fit(gaussian_2d, x, s, p0=[stack[frame, yc, xc],
                                                      rv,
                                                      math.floor(ws / 2),
                                                      math.floor(ws / 2)])
        
        # Output        
        x_fit = x0 + popt[2]
        y_fit = y0 + popt[3]
        print(f"frame, x, y: {frame}, {x_fit}, {y_fit}")
    
napari.run()

出力

Pointを配置して、「F」を押すと、、、

フィッティングした座標が出力された！

ポイント

フィッティング領域の切り出し

指定したPointの周りの領域の輝度を抽出。Fittingのために1次元化している。

s = np.array([stack[frame, y0 + int(i / ws), x0 + i % ws] for i in range(ws * ws)])

フィッティング

初期値を与えてscipyのcurve_fitでフィッティングする。

popt, pcov = curve_fit(gaussian_2d, x, s, p0=[stack[frame, yc, xc],
                                              rv,
                                              math.floor(ws / 2),
                                              math.floor(ws / 2)])

2次元ガウシアン関数

簡単のために、X方向とY方向のvarianceは等しいと仮定。

def gaussian_2d(x, a, b, c, d):
    return a * np.exp(-b * ((x % ws - c) ** 2 + ((x / ws).astype(int) - d) ** 2))