分散共分散行列で定義された楕円体の断面算出
ある生物の神経細胞画像データを機械学習解析で使える形として抽出する。
アノテーションのデータがMATLABで解析されたファイルで与えられていた。
github.com
figshare.com
.matファイルから辞書データとして抜き出したところ、細胞の形が中心の座標と三次元の分散共分散行列で表されていた。
# 例 # 中心座標 center = [50, 50, 0] # 分散共分散行列 # [x_var, xy_covar, zx_covar], # [xy_covar, y_var, yz_covar], # [zx_covar, yz_covar, z_var], cov_matrix = np.array( [ [37.67218941864915, 9.97966493551884, 2.1263842045250754], [9.97966493551884, 14.353317888367672, 5.0364599497266935], [2.1263842045250754, 5.0364599497266935, 15.805485755173178], ] )
描画するとこんな感じ。
def show_ellipsoid(center, cov_matrix): # 固有値と固有ベクトルを計算 eigenvalues, eigenvectors = np.linalg.eig(cov_matrix) # 楕円体を描画するためのパラメータ u = np.linspace(0, 2 * np.pi, 100) v = np.linspace(0, np.pi, 100) x = np.outer(np.cos(u), np.sin(v)) y = np.outer(np.sin(u), np.sin(v)) z = np.outer(np.ones(np.size(u)), np.cos(v)) # 固有値に応じて楕円体を変形 for i in range(len(eigenvalues)): x = np.sqrt(eigenvalues[i]) * (eigenvectors[0, i] * x) y = np.sqrt(eigenvalues[i]) * (eigenvectors[1, i] * y) z = np.sqrt(eigenvalues[i]) * (eigenvectors[2, i] * z) # 中心座標を追加 x += center[0] y += center[1] z += center[2] # 3Dプロットを作成 fig = plt.figure() ax = fig.add_subplot(111, projection="3d") # 楕円体をプロット ax.plot_surface(x, y, z, color="b", alpha=0.2) # 中心座標をプロット ax.scatter(center[0], center[1], center[2], color="r", s=100) # 軸ラベルの設定 ax.set_xlabel("X") ax.set_ylabel("Y") ax.set_zlabel("Z") # グラフの表示 plt.show()
顕微鏡で撮影された画像では上記のような楕円体の細胞が任意のzでスライスした楕円として観察される。
ChatGPTに助けてもらって関数を実装した。
def get_ellipse_intersection(center, cov_matrix, z): # 平面との交点で分散共分散行列を更新 shifted_cov_matrix = cov_matrix.copy() shifted_cov_matrix[:2, :2] -= (z - center[2]) ** 2 * np.linalg.inv( cov_matrix[:2, :2] ) # 2次元の分散共分散行列を抽出 cov_2d = shifted_cov_matrix[:2, :2] # 固有値と固有ベクトルを計算 eigenvalues, eigenvectors = np.linalg.eig(cov_2d) # 楕円のパラメータを計算 semi_major_axis = np.sqrt(eigenvalues[0]) semi_minor_axis = np.sqrt(eigenvalues[1]) angle = np.arctan2(eigenvectors[1, 0], eigenvectors[0, 0]) return semi_major_axis, semi_minor_axis, angle def add_ellipse(canvas, x, y, semi_major_axis, semi_minor_axis, angle): center = (int(x), int(y)) # 中心座標を整数に変換 axes = (int(semi_major_axis), int(semi_minor_axis)) # 軸の長さを整数に変換 angle_degrees = int(np.degrees(angle)) # 角度を度数法に変換 cv2.ellipse( canvas, center, axes, angle_degrees, 0, 360, 255, thickness=-1, )
結果
多分あってるはず...
コードの全体図は以下の通りです。
import numpy as np import matplotlib.pyplot as plt import cv2 import math import os def show_ellipsoid(center, cov_matrix): # 固有値と固有ベクトルを計算 eigenvalues, eigenvectors = np.linalg.eig(cov_matrix) # 楕円体を描画するためのパラメータ u = np.linspace(0, 2 * np.pi, 100) v = np.linspace(0, np.pi, 100) x = np.outer(np.cos(u), np.sin(v)) y = np.outer(np.sin(u), np.sin(v)) z = np.outer(np.ones(np.size(u)), np.cos(v)) # 固有値に応じて楕円体を変形 for i in range(len(eigenvalues)): x = np.sqrt(eigenvalues[i]) * (eigenvectors[0, i] * x) y = np.sqrt(eigenvalues[i]) * (eigenvectors[1, i] * y) z = np.sqrt(eigenvalues[i]) * (eigenvectors[2, i] * z) # 中心座標を追加 x += center[0] y += center[1] z += center[2] # 3Dプロットを作成 fig = plt.figure() ax = fig.add_subplot(111, projection="3d") # 楕円体をプロット ax.plot_surface(x, y, z, color="b", alpha=0.2) # 中心座標をプロット ax.scatter(center[0], center[1], center[2], color="r", s=100) # 軸ラベルの設定 ax.set_xlabel("X") ax.set_ylabel("Y") ax.set_zlabel("Z") # グラフの表示 plt.show() def get_ellipse_intersection(center, cov_matrix, z): # 平面との交点で分散共分散行列を更新 shifted_cov_matrix = cov_matrix.copy() shifted_cov_matrix[:2, :2] -= (z - center[2]) ** 2 * np.linalg.inv( cov_matrix[:2, :2] ) # 2次元の分散共分散行列を抽出 cov_2d = shifted_cov_matrix[:2, :2] # 固有値と固有ベクトルを計算 eigenvalues, eigenvectors = np.linalg.eig(cov_2d) # 楕円のパラメータを計算 semi_major_axis = np.sqrt(eigenvalues[0]) semi_minor_axis = np.sqrt(eigenvalues[1]) angle = np.arctan2(eigenvectors[1, 0], eigenvectors[0, 0]) return semi_major_axis, semi_minor_axis, angle def add_ellipse(canvas, x, y, semi_major_axis, semi_minor_axis, angle): center = (int(x), int(y)) # 中心座標を整数に変換 axes = (int(semi_major_axis), int(semi_minor_axis)) # 軸の長さを整数に変換 angle_degrees = int(np.degrees(angle)) # 角度を度数法に変換 cv2.ellipse( canvas, center, axes, angle_degrees, 0, 360, 255, thickness=-1, ) if __name__ == "__main__": os.makedirs("test_figs", exist_ok=True) os.makedirs("test_imgs", exist_ok=True) # 中心座標 center = [50, 50, 0] # 分散共分散行列 # [x_var, xy_covar, zx_covar], # [xy_covar, y_var, yz_covar], # [zx_covar, yz_covar, z_var], cov_matrix = np.array( [ [37.67218941864915, 9.97966493551884, 2.1263842045250754], [9.97966493551884, 14.353317888367672, 5.0364599497266935], [2.1263842045250754, 5.0364599497266935, 15.805485755173178], ] ) show_ellipsoid(center, cov_matrix) # 任意のz座標で交わる断面の楕円を求める for i, z in enumerate(range(-15, 15)): semi_major_axis, semi_minor_axis, angle = get_ellipse_intersection( center, cov_matrix, z ) # 楕円を描画 theta = np.linspace(0, 2 * np.pi, 100) x = semi_major_axis * np.cos(theta) y = semi_minor_axis * np.sin(theta) x_rotated = x * np.cos(angle) - y * np.sin(angle) y_rotated = x * np.sin(angle) + y * np.cos(angle) plt.figure(figsize=(5, 5)) plt.plot(x_rotated + center[0], y_rotated + center[1]) plt.xlim(40, 60) plt.ylim(40, 60) plt.xlabel("X") plt.ylabel("Y") plt.title(f"Ellipse Intersection with z={z} Plane") plt.grid(True) plt.savefig(f"test_figs/{i:03}.png") plt.close() xy_canvas = np.zeros((100, 100), dtype="uint8") if ( math.isnan(semi_major_axis) == False and math.isnan(semi_minor_axis) == False ): add_ellipse( xy_canvas, center[0], center[1], semi_major_axis, semi_minor_axis, angle ) cv2.imwrite(f"test_imgs/{i:03}.png", xy_canvas)
