今回は線形最小二乗法を使用して、以下の細胞形状を楕円フィッテイングしてみます。
輪郭座標は非極大抑制、ヒステリシス処理を行い、つながった座標の集合を取得しています。
目次
楕円フィッティングの意義
細胞形状を楕円フィッテイングすることによって、以下のようなメリットがあります。
・細胞を数式で表すことにより、形態学的情報を一意に記述できる。(ただし、近似式である）
・パラメータから細胞の伸長方向、及び大きさ（主軸、副軸）を瞬時に計算することができる。
・細胞内生理状態解析の規格化ができる。
線形最小二乗法
フィッティングに関して、以下のように各輪郭座標と楕円の座標のL2ノルムの二乗誤差の和を最小化するようにn個のパラメータを決定します。
$$\begin{align} \:\:\:E &= \sum_{i=1}^N (f_i-\hat{f_i})^2 \\ &=\sum_{i=1}^N (f_i-\mathbf{c}^\mathrm{T}\phi(\mathbf{x_i}))^2\:\:(\mathbf{c}\in \mathbb{R}^n )\end{align}$$
上記の二乗誤差の和を最小化する方法として、勾配降下法、ニューラルネットワークを使用した二乗誤差の最小化、正規方程式等が挙げられますが、今回は解析的に解くことができる正規方程式を使用した方法を採用します。
よって、
$$\frac{\partial }{\partial \mathbf{c}}\sum_{i=1}^N (f_i-\mathbf{c}^\mathrm{T}\phi(\mathbf{x_i}))^2 = 0$$
について考えて、最小ノルム解を与えるパラメータを決定します。
線形最小二乗フィッティング
楕円の一般式は以下のように与えられます。
$Ax^2+Bxy+Cy^2+Dx+Ey+F=0$
ここで、
$$-C_1x^2-C_2xy+y^2-C_3x-C_4y-C_5=0$$
$$ \forall C_n\in \mathbb{R}$$
とすると、
$\mathbf{W} = \begin{pmatrix}x_1^2 &x_1y_1 & x_1&y_1&1\\\vdots & \vdots &\vdots&\vdots&\vdots\\x_n^2 &x_ny_n & x_n&y_n&1\\\end{pmatrix} \:\in \mathbb{R}^{n\times 5}$
$\mathbf{c} = \begin{pmatrix}C_1 &C_2&C_3&C_4&C_5\\\end{pmatrix}^\mathrm{T}\:\in \mathbb{R}^5$
$\mathbf{f} = \begin{pmatrix}y_1^2 &\dots&y_n^2 \\\end{pmatrix}^\mathrm{T}\:\in \mathbb{R}^{n}$
より、以下の正規方程式を得ます。
$\mathbf{W}\mathbf{c}= \mathbf{f}$
よって、
$\mathbf{c}=(\mathbf{W}^\mathrm{T}\mathbf{W})^{-1}\mathbf{W}^\mathrm{T}\mathbf{f} \:\: : \exists (\mathbf{W}^\mathrm{T}\mathbf{W})^{-1}$
より、
$\mathbf{W}^\mathrm{T}\mathbf{W} \in \mathbb{R}^{n\times n}$が正則の場合、n個の入力座標によりL2ノルムを正則化項として用いた楕円の最小二乗フィッティングを行うことができます。
偏微分を使用したパラメータの決定
上記では正規方程式を作成しましたが、$(\mathbf{W}^\mathrm{T}\mathbf{W})^{-1}$の計算に時間がかかるため、二乗誤差の和をパラメータベクトルで偏微分した勾配を0に固定する方法でアプローチします。
$$\frac{\partial }{\partial \mathbf{c}}\sum_{i=1}^N (f_i-\mathbf{c}^\mathrm{T}\phi(\mathbf{x_i}))^2 = 0$$
上記の式は以下の様に書き換えることができます。
$$\frac{\partial }{\partial \mathbf{c}}\sum_{i=1}^N ( x^2+C_1xy+C_2y^2+C_3x+C_4y+C_5)^2 = 0$$
この式を展開して行列表示すると、
$$\Sigma =\begin{pmatrix}\sum{x_i^2y_j^2} & \sum{x_iy_j^2} & \sum{x_i^2y_j} & \sum{x_iy_j^2} &\sum{x_iy_j}\\ \sum{x_iy_j^3} & \sum{y_j^4} & \sum{x_iy_j^2} & \sum{y_j^3} &\sum{y_j^2}\\ \sum{x_i^2y_j} &\sum{x_iy_j^2} &\sum{x_i^2} & \sum{x_iy_j} & \sum{x_i}\\ \sum{x_iy_j^2} & \sum{y_j^3} & \sum{x_iy_j} &\sum{y_j^2} & \sum{y_j}\\ \sum{x_iy_j}&\sum{y_j^2}&\sum{x_i}&\sum{y_j}&\sum{1}
\end{pmatrix}$$
$$\mathbf{f} = \begin{pmatrix}-\sum x_i^3y_i\\-\sum x_i^2y_i^2\\-\sum x_i^3\\-\sum x_i^1y_i\\-\sum x_i^2\end{pmatrix}$$
のもとで、
$$\mathbf{c}=\Sigma^{-1}\mathbf{f}$$
を解けば良いことになります。
スペクトル分解による幾何学情報の抽出
ここでは、固有値分解を用いて細胞輪郭座標から計算した楕円の近似式の基底を取り直すことにより、規格化を行い解析の効率化を目指します。
例えば、基底変換前の楕円の軸はxy軸に直交していないため、幾何学的な情報を取り出す際に数値的な手法が必要となってきますが、楕円の各軸がxy軸に直交するように基底を取り直した場合は主要な特性を解析的に求めることができます。
線形最小二乗法から、楕円の一般式の全てのパラメータ$\mathbf{c}\in \mathbb{R}^5$を得ることができたので、これらを使用して一般式を整理すると
$$-C_1x^2-C_2xy+y^2-C_3x-C_4y-C_5=0$$
よって、回転項に着目して実対称行列を以下のように定義します。
$\mathbf{A} = \begin{bmatrix}-C_1&-\frac{C_2}{2}\\-\frac{C_2}{2}&1 \\ \end{bmatrix}$
ここで、固有方程式
$\:\:\det(\mathbf{A}-\lambda \mathbf{I}) = 0$
を解き、以下のように得られた固有値に対応するそれぞれの固有空間
$$\mathbf{W}_{\lambda _i} = \{c\mathbf{v}_i|c\in \mathbb{R},\lambda_i\notin \mathbb{C}\}$$
のカーネルの基底を利用して、実対称行列Aのスペクトル分解を求めます。
$\mathbf{A} = \mathbf{Q}\mathbf{\Lambda} \mathbf{Q}^\mathrm{T}$
$\mathbf{Q} = (\mathbf{v_1}\:,\:\: \mathbf{v_2})$
$\mathbf{\Lambda} = \begin{bmatrix} \lambda_1 &0\\ 0&\lambda_2\\ \end{bmatrix}$
よって、二次形式は上記のスペクトル分解を用いて
$\mathbf{X}^\mathrm{T} \mathbf{A}\mathbf{X}- \begin{pmatrix}C_3&C_4\end{pmatrix}\mathbf{X} -C_5 =\mathbf{X}^\mathrm{T}\mathbf{Q}\mathbf{\Lambda} \mathbf{Q}^\mathrm{T}\mathbf{X}-\begin{pmatrix}C_3&C_4\end{pmatrix}\mathbf{X}-C_5= 0$
上記から、基底変換を行うと
$\mathbf{U}^\mathrm{T}\mathbf{\Lambda}\mathbf{U} – \begin{pmatrix} C_3&C_4\end{pmatrix}(\mathbf{Q}^\mathrm{T})^{-1}\mathbf{U}^\mathrm{T} -C_5= 0$
この時、
$ker\mathbf{W}_{\lambda _1}\perp ker\mathbf{W}_{\lambda _2}$
より、新規座標基底は直交します。
ここで、基底変換について、回転項のみを扱ったため平行移動について考慮すると、
基底変換前の楕円の中心座標は
$\begin{pmatrix} c_x&c_y\end{pmatrix}^\mathrm{T} = \mathbf{A}^{-1}\begin{pmatrix} C_3&C_4\end{pmatrix}^\mathrm{T}$
基底変換後の楕円の中心座標は線型結合を用いて、
$\begin{pmatrix} c_{u_1}&c_{u_2}\end{pmatrix}^\mathrm{T}=\begin{pmatrix} c_xv_{1_x}+c_yv_{1_y}&c_xv_{2_x}+c_yv_{2_y}\end{pmatrix}^\mathrm{T}$
と表されます。
よって、この中心座標に基づいて基底変換後のすべての座標を平行移動させることにより、中心を原点にすることが可能です。
また、その他の座標についても
$\mathbf{U} = \mathbf{X}^\mathrm{T}\mathbf{Q}$
より基底を変換することができます。
基底変換結果
フィッテイング後の楕円を細胞輪郭のプロット上に重ねたグラフを示します。
楕円が途切れているのは描画の問題であり、楕円の式の各種パラメータ は取得できているため問題ありません。また、⻑軸と短軸が直交していないように見えますが、これは 座標の軸の幅がxとyで異なることに起因します。
次に、スペクトル分解による基底変換後の楕円及び細胞輪郭座標を示します。
メッシュによるフィット楕円の細分化
細胞を一定条件で細分化することで、細胞内の蛍光輝度情報解析を規格化できる可能性があります。
まずは解析の規格化の可能性を探るため、基底変換後のフィット楕円をメッシュに分割してみます。
基底変換後のフィット楕円は
$\mathbf{U}^\mathrm{T}\mathbf{\Lambda}\mathbf{U} – \begin{pmatrix} C_3&C_4\end{pmatrix}(\mathbf{Q}^\mathrm{T})^{-1}\mathbf{U}^\mathrm{T} -C_5= 0$
で表されるため、各種パラメータを使用してメッシュ情報を計算します。
ここで、
$\mathbf{U} = \mathbf{X}^\mathrm{T}\mathbf{Q}$
を利用して規定の逆変換を行うと、
メッシュによる輪郭座標の細分化
上記でフィット楕円からメッシュ情報を計算して、細胞形状を細分化するアルゴリズムについて考えましたが、ここではそのアルゴリズムを応用して輪郭座標をメッシュで細分化します。
この時、 輪郭座標は解析的には記述できないため、基底変換後の座標のu1の範囲を利用してメッシュ情報を計算する方法について考えます。
例えば以下のようなアルゴリズムによって、u1の値に対して近似的にメッシュの両端座標を求めることができます。
@staticmethod --- Time Complexity: Θ(n)
def find_closest(t:float,l:list]) -> list]:
d:float, dst:list[float] = ∞, []
for i in l:
if abs(i[0]-t) < d: dst = i
d = abs(i[0]-t)
return dst上記のアルゴリズムを使用した結果、以下のようなメッシュを得ることができました。
このメッシュ情報を元に、基底を逆変換すると
細胞分裂検出への応用
ここまでで、輪郭座標を入力として５つのパラメータを決定して楕円フィッティングを行うことができたので、このアルゴリズムを細胞増殖検出へ応用できないか検討します。
はじめに、上記のアルゴリズムが細胞形状の経時変化を上手く追跡することができるかを検証し、その時の各フィッティングパラメータの変化を確認します。
ここでは、0から20までのフレームを追跡してみます。
import cv2
import matplotlib.pyplot as plt
import numpy as np
from main2 import LeastSquares
def get_contours(img_num:str):
try:
#! Read img
img = cv2.imread(f"images/{img_num}.png")
img = img[250:400,280:430]
#! Gray scale
img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
print(len(img_gray))
#! Bin
#cell6
ret,thresh = cv2.threshold(img_gray,80,255,cv2.THRESH_BINARY)
# cv2.imwrite("out1.png",thresh)
#! Canny edge detection
img_canny = cv2.Canny(thresh,0,55)
# cv2.imwrite("out1.png",img_canny)
#! Find contours
contours, hierarchy = cv2.findContours(img_canny,cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
#! Print out areas for each contours
for i in range(len(contours)):
print('Shape_{}'.format(i),'area', cv2.contourArea(contours[i]))
#! Exclude contours that are larer than 2000
contours = list(filter(lambda x: cv2.contourArea(x) >= 100 , contours))
#! Draw contours
img_contours = cv2.drawContours(img, contours, -1, (0, 255, 0),2)
cv2.imwrite(f"out_{img_num}.png",img_contours)
#! Retrieve coordinates ([[x_1,y_1],,,[x_n,y_n]])
contours_coordinates = [[i[0][0],i[0][1]] for i in contours[0]]
print(contours_coordinates)
return contours_coordinates
except Exception as e:
print(e)
import time
with open("result.txt","w") as fpout:
fpout.write("C1,C2,C3,C4,C5,R2\n")
for i in range(1,21):
time.sleep(2)
if __name__ == "__main__":
contours_coordinates = get_contours(f"{i}")
L = LeastSquares(contours_coordinates,i)
L._solve2()
L.base_conversion()
L.draw()
fpout.write(",".join([str(i) for i in L.get_params()])+ "," + str(L.get_R2()) +"\n")
print(L)
import matplotlib.pyplot as plt
with open("text.txt","r") as fp:
lines = [i.split(",") for i in fp.readlines()]
data = lines[1:]
print(data)
x = [i for i in range(len(data))]
C_1s = [float(i[0].replace("\n","")) for i in data]
C_2s = [float(i[1].replace("\n","")) for i in data]
C_3s = [float(i[2].replace("\n","")) for i in data]
C_4s = [float(i[3].replace("\n","")) for i in data]
C_5s = [float(i[4].replace("\n","")) for i in data]
R2 = [float(i[5].replace("\n","")) for i in data]
fig = plt.figure()
label = "R2"
plt.scatter(x,R2)
plt.plot(x,R2,label = label)
plt.grid()
plt.xlabel("Frame")
plt.ylabel(label)
plt.legend()
fig.savefig(f"{label}.png",dpi = 500)実行結果
画像が多すぎるので、特徴的な形状のフレームを３種類選びました。
上記の処理で、各パラメータの値を取得できているので、その値をプロットしてみます。
ついでに、決定係数R2の変遷も見てみます。
ちょうど細胞が分裂したフレーム11のときに特徴的なピークが確認できます。
細胞形状の多項式回帰
細胞形状を解析的に記述できない問題に対して、最小二乗法による多項式回帰を用いて解決を試みます。
例えば、細分化メッシュ群は各輪郭座標に基づいており、細胞の軸に直交していないので、定量的な解析を行いたい場合に影響が出る可能性があります。
よって、細胞輪郭を、多項式を用いて関数で記述できるようにすることを目指します。
この時、回帰後の多項式はベイズ情報量規準によって評価します。
近似多項式に関して、
$$\mathbf{w} = \begin{pmatrix} w_0 …w_n &\epsilon \end{pmatrix}^\mathrm{T}\forall{w_n}\in\mathbb{R}$$
$$\mathbf{x} = \begin{pmatrix}x_0…x^{n-1}\end{pmatrix}^\mathrm{T}$$
のもとで、
$$\hat{\theta}(\mathbf{w},\mathbf{x}) = \mathbf{w}^\mathrm{T}\mathbf{x} + \epsilon$$
と定義します。
この時、
$\mathbf{A} = \begin{pmatrix}x_{0_0} &\cdots & x_{n_0}^{n-1} & 1 \\\vdots &\cdots & \vdots & \vdots \\ x_{0_i} &\cdots & x_{n_i}^{n-1} & 1\end{pmatrix}$
$\mathbf{f}=\begin{pmatrix}\theta(\mathbf{w,\mathbf{x}_0})\\ \vdots \\ \theta(\mathbf{w,\mathbf{x}_i})\end{pmatrix}$
とすると、多項式の係数パラメータはそれぞれ
$\mathbf{w} =(\mathbf{A}^\mathrm{T}\mathbf{A})^{-1}\mathbf{A}^\mathrm{T}\mathbf{f}$
となります。
今回は次数19までの多項式を利用して細胞輪郭座標を解析的に記述できないか試します。
上記の式を計算すると、
$BIC = n\ln (\hat{\sigma}^2)+k\ln(n)$
ここで、上式で定義されるベイズ情報量規準を使用して多項式を評価すると、
それぞれ次数k=18,k=8の時が最適であるということがわかりました。
これらの情報を利用して細胞輪郭を再構築すると、以下の様になります。
よって、基底に直交するメッシュを解析的に生成することができるようになりました。k=8だと少しずれている気もしますが、BICによると最善のようです。
アルゴリズムの実装
計算速度より実装を優先したいため、Pythonを使用します。
取得した輪郭座標をコンストラクタに渡すことで、あとは全自動で解析を行うことができます。
import cv2
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import solve, eig, norm
import matplotlib.pyplot as plt
import math
def get_contours(img_num:str):
try:
#! Read img
img = cv2.imread(f"images/{img_num}.png")
img = img[250:400,280:430]
#! Gray scale
img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
print(len(img_gray))
#! Bin
#cell6
ret,thresh = cv2.threshold(img_gray,80,255,cv2.THRESH_BINARY)
# cv2.imwrite("out1.png",thresh)
#! Canny edge detection
img_canny = cv2.Canny(thresh,0,55)
# cv2.imwrite("out1.png",img_canny)
#! Find contours
contours, hierarchy = cv2.findContours(img_canny,cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
#! Print out areas for each contours
for i in range(len(contours)):
print('Shape_{}'.format(i),'area', cv2.contourArea(contours[i]))
#! Exclude contours that are larer than 2000
contours = list(filter(lambda x: cv2.contourArea(x) >= 100 , contours))
#! Draw contours
img_contours = cv2.drawContours(img, contours, -1, (0, 255, 0),2)
cv2.imwrite(f"out_{img_num}.png",img_contours)
#! Retrieve coordinates ([[x_1,y_1],,,[x_n,y_n]])
contours_coordinates = [[i[0][0],i[0][1]] for i in contours[0]]
print(contours_coordinates)
return contours_coordinates
except Exception as e:
print(e)
class LeastSquares:
@staticmethod
def find_closest(t:float,l) -> list[float]:
d:float = float("inf")
dst:list[float] = []
for i in l:
if abs(i[0]-t) < d:
dst = i
d = abs(i[0]-t)
return dst
def __init__(self, contours_cds,n:int) -> None:
self.n = n
#!輪郭座標
self.x, self.y = [i[0] for i in contours_cds], [i[1] for i in contours_cds]
#!正規方程式のためのコンポーネント
self.W: np.ndarray = np.array(
[np.array([i**2, i*j, i, j, 1.0]) for i, j in zip(self.x, self.y)])
self.W_t: np.ndarray = np.transpose(self.W)
self.f: np.ndarray = np.transpose(np.array([i**2 for i in self.y]))
#!0で初期化した各パラメータ
self.C_1, self.C_2, self.C_3, self.C_4, self.C_5 = 0.0, 0.0, 0.0, 0.0, 0.0
#!基底変換前楕円のxに対するy座標
self.y_p: np.ndarray
self.y_n: np.ndarray
#!基底変換前の楕円の中心座標
self.cx: float
self.cy: float
#!基底変換前の楕円と楕円の長軸と短軸の傾き
self.m_long: float
self.m_short: float
#!基底変換前の楕円と楕円の長軸との交点２つ
self.intersec_a1: list[float]
self.intersec_a2: list[float]
#!基底変換前の楕円と楕円の短軸との交点２つ
self.intersec_a3: list[float]
self.intersec_a4: list[float]
#!スペクトル分解のコンポーネント
self.Q: np.ndarray
#!基底変換後の楕円と楕円の長軸との交点２つ
self.intersec_a1_u: list[float]
self.intersec_a2_u: list[float]
#!基底変換後の楕円と楕円の短軸との交点２つ
self.intersec_a3_u: list[float]
self.intersec_a4_u: list[float]
#!基底変換後の楕円の座標
self.u1_s: list[float]
self.u2_ps: list[float]
self.u2_ns: list[float]
#!基底変換後の輪郭座標
self.c1s: list[float]
self.c2s: list[float]
#!細分化の本数
self.num_subdiv: int = 30
#!基底変換後の楕円の細分化後のメッシュ座標
self.u1_s_subdiv: list[float]
self.u2_ns_subdiv: list[float]
self.u2_ps_subdiv: list[float]
#!基底逆変換後の楕円の細分化後のメッシュ座標
self.x_ps_subdiv: list[float]
self.x_ns_subdiv: list[float]
self.y_ps_subdiv: list[float]
self.y_ns_subdiv: list[float]
#!基底変換後の輪郭の細分化後のメッシュ座標
self.c1u_s_subdiv: list[float]
self.c2u_ns_subdiv: list]
self.c2u_ps_subdiv: list]
#!基底逆変換後の輪郭の細分化後のメッシュ座標
self.cx_ps_subdiv: list[float]
self.cx_ns_subdiv: list[float]
self.cy_ps_subdiv: list[float]
self.cy_ns_subdiv: list[float]
#!R2係数計算のための数値 (全部リスト、リストfloatにする。)
self.yp_R2: list]
self.yn_R2: list]
self.yp_R2_hat:list]
self.yn_R2_hat:list]
self.R2: float
def _solve(self) -> np.ndarray:
left, right = self.W_t@self.W, self.W_t@self.f
d = solve(left, right)
self.C_1, self.C_2, self.C_3, self.C_4, self.C_5 = d[0], d[1], d[2], d[3], d[4]
#!中心座標の計算
A = np.array([[-self.C_1, -self.C_2/2], [-self.C_2/2, 1]])
A_inv = np.linalg.inv(A)
center = -A_inv@np.array([-self.C_3, -self.C_4])/2
self.cx, self.cy = center[0], center[1]
return d
def _solve2(self) -> np.ndarray:
S_11, S12, S13, S14, S15 = 0, 0, 0, 0, 0
S_21, S22, S23, S24, S25 = 0, 0, 0, 0, 0
S_31, S32, S33, S34, S35 = 0, 0, 0, 0, 0
S_41, S42, S43, S44, S45 = 0, 0, 0, 0, 0
S_51, S52, S53, S54, S55 = 0, 0, 0, 0, 0
f1, f2, f3, f4, f5 = 0, 0, 0, 0, 0
for i, j in zip(self.x, self.y):
S_11 += i**2*j**2
S_21 += i*j**3
S_31 += i**2*j
S_41 += i*j**2
S_51 += i*j
S12 += i*j**3
S22 += j**4
S32 += i*j**2
S42 += j**3
S52 += j**2
S13 += i**2*j
S23 += i*j**2
S33 += i**2
S43 += i*j
S53 += i
S14 += i*j**2
S24 += j**3
S34 += i*j
S44 += j**2
S54 += j
S15 += i*j
S25 += j**2
S35 += i
S45 += j
S55 += 1
f1 += i**3*j
f2 += i**2*j**2
f3 += i**3
f4 += i**2*j
f5 += i**2
Sigma = np.array([[S_11, S12, S13, S14, S15], [S_21, S22, S23, S24, S25], [S_31, S32, S33, S34, S35], [S_41, S42, S43, S44, S45], [S_51, S52, S53, S54, S55]])
f = np.transpose(np.array([-f1, -f2, -f3, -f4, -f5]))
d = np.linalg.inv(Sigma)@f
A, B, C, D, E = d[0], d[1], d[2], d[3], d[4]
self.C_1 = -1/B
self.C_2 = A*self.C_1
self.C_3 = C*self.C_1
self.C_4 = D*self.C_1
self.C_5 = E*self.C_1
print((A*D-2*B*C)/(4*B-A**2))
self.cx, self.cy = (A*D-2*B*C)/(4*B-A**2),(A*C-2*D)/(4*B-A**2)
return np.array([A, B, C, D, E])
def base_conversion(self) -> None:
fig = plt.figure(figsize=[5, 5])
a, b, c = -self.C_1, -self.C_2, 1
#!実対称行列と固有値計算
A = np.array([[a, b/2], [b/2, c]])
eigs = eig(A)
#!固有空間の基底
if eigs[0][0] < eigs[0][1]:
v1, v2 = np.array([eigs[1][0][0], eigs[1][1][0]]), np.array([eigs[1][0][1], eigs[1][1][1]])
LAMBDA_1, LAMBDA_2 = eigs[0][0], eigs[0][1]
self.m_long, self.m_short = v1[1]/v1[0], v2[1]/v2[0]
else:
v1, v2 = np.array([eigs[1][0][1], eigs[1][1][1]]), np.array([eigs[1][0][0], eigs[1][1][0]])
LAMBDA_1, LAMBDA_2 = eigs[0][1], eigs[0][0]
self.m_long, self.m_short = v1[1]/v1[0], v2[1]/v2[0]
v1_x, v1_y = v1
v2_x, v2_y = v2
#!基底変換前の楕円と角軸の交点
m = self.m_long
k = self.cy-self.cx*self.m_long
a = -self.C_1-self.C_2*self.m_long+self.m_long**2
b = -self.C_2*k + 2*k*m-self.C_3-self.C_4*m
c = k**2-self.C_4*k-self.C_5
a_x, b_x = (-b+np.sqrt(b**2-4*a*c))/2/a, (-b-np.sqrt(b**2-4*a*c))/2/a
a_y, b_y = m*a_x+k, m*b_x+k
self.intersec_a1, self.intersec_a2 = [a_x, a_y], [b_x, b_y]
m = self.m_short
k = self.cy-self.cx*m
a = -self.C_1-self.C_2*m+m**2
b = -self.C_2*k + 2*k*m-self.C_3-self.C_4*m
c = k**2-self.C_4*k-self.C_5
a_x, b_x = (-b+np.sqrt(b**2-4*a*c))/2/a, (-b-np.sqrt(b**2-4*a*c))/2/a
a_y, b_y = m*a_x+k, m*b_x+k
self.intersec_a3, self.intersec_a4 = [a_x, a_y], [b_x, b_y]
#!基底変換のための各種パラメータ
v1_x, v1_y = v1
v2_x, v2_y = v2
self.Q = np.array([[v1_x,v2_x],[v1_y,v2_y]])
a = LAMBDA_1
b = LAMBDA_2
c_u1, c_u2 = self.cx*v1_x + self.cy*v1_y, self.cx*v2_x + self.cy*v2_y
#!輪郭の基底変換（線型結合）
c1s = [(i*v1_x + j*v1_y - c_u1) for i, j in zip(self.x, self.y)]
c2s = [(i*v2_x + j*v2_y - c_u2) for i, j in zip(self.x, self.y)]
self.c1s,self.c2s = c1s, c2s
#!交点の基底変換 (線型結合）
self.intersec_a1_u = [self.intersec_a1[0]*v1_x + self.intersec_a1[1]* v1_y-c_u1, self.intersec_a1[0]*v2_x + self.intersec_a1[1]*v2_y-c_u2]
self.intersec_a2_u = [self.intersec_a2[0]*v1_x + self.intersec_a2[1]* v1_y-c_u1, self.intersec_a2[0]*v2_x + self.intersec_a2[1]*v2_y-c_u2]
self.intersec_a3_u = [self.intersec_a3[0]*v1_x + self.intersec_a3[1]* v1_y-c_u1, self.intersec_a3[0]*v2_x + self.intersec_a3[1]*v2_y-c_u2]
self.intersec_a4_u = [self.intersec_a4[0]*v1_x + self.intersec_a4[1]* v1_y-c_u1, self.intersec_a4[0]*v2_x + self.intersec_a4[1]*v2_y-c_u2]
#!基底変換後の楕円のプロット
a = max([self.intersec_a1_u[0],self.intersec_a2_u[0],self.intersec_a3_u[0],self.intersec_a4_u[0]])
b = max([self.intersec_a1_u[1],self.intersec_a2_u[1],self.intersec_a3_u[1],self.intersec_a4_u[1]])
self.u1_s = np.linspace(-a, a, 1000)
self.u2_ps = [b/a*np.sqrt(a**2-i**2) for i in self.u1_s]
self.u2_ns = [-b/a*np.sqrt(a**2-i**2) for i in self.u1_s]
#??基底変換後の楕円の細分化（メッシュ情報の計算）
self.u1_s_subdiv = np.linspace(-a, a, self.num_subdiv)
self.u2_ps_subdiv = [b/a*np.sqrt(a**2-i**2) for i in self.u1_s_subdiv]
self.u2_ns_subdiv = [-b/a*np.sqrt(a**2-i**2) for i in self.u1_s_subdiv]
#!楕円のメッシュによる細分化プロット （基底変換後)
# plt.plot([self.u1_s_subdiv[0],self.u1_s_subdiv[0]],[self.u2_ns_subdiv[0],self.u2_ps_subdiv[0]],color =[0.0, 1.0, 0.0],label = f"Subdivision Mesh (n={self.num_subdiv})")
# for i in range(1,len(self.u1_s_subdiv)-1):
# plt.scatter([self.u1_s_subdiv[i],self.u1_s_subdiv[i]],[self.u2_ns_subdiv[i],self.u2_ps_subdiv[i]],color =[0.0, 1.0, 0.0],s = 15)
# plt.plot([self.u1_s_subdiv[i],self.u1_s_subdiv[i]],[self.u2_ns_subdiv[i],self.u2_ps_subdiv[i]],color =[0.0, 1.0, 0.0])
#!基底変換前の楕円の細分化（メッシュ情報の計算)
a,b,c,d = v1_x,v1_y,v2_x,v2_y
self.x_ps_subdiv = [(d*i-b*j)/(a*d-b*c)+self.cx for i,j in zip(self.u1_s_subdiv,self.u2_ps_subdiv)]
self.x_ns_subdiv = [(d*i-b*j)/(a*d-b*c)+self.cx for i,j in zip(self.u1_s_subdiv,self.u2_ns_subdiv)]
self.y_ps_subdiv = [(c*i-a*j)/(c*b-a*d)+self.cy for i,j in zip(self.u1_s_subdiv,self.u2_ps_subdiv)]
self.y_ns_subdiv = [(c*i-a*j)/(c*b-a*d)+self.cy for i,j in zip(self.u1_s_subdiv,self.u2_ns_subdiv)]
#!基底変換後の輪郭の細分化（メッシュ情報の計算）計算量n**2→よくない
r = max(self.c1s)
i = 0
self.c1u_s_subdiv = [-r + 2*r/self.num_subdiv*i for i in range(1,self.num_subdiv)]
tmp_p = [[i,j] for i,j in zip(self.c1s,self.c2s) if j >= 0 ]
tmp_n = [[i,j] for i,j in zip(self.c1s,self.c2s) if j < 0 ]
self.c2u_ps_subdiv = [self.find_closest(i,tmp_p) for i in self.c1u_s_subdiv]
self.c2u_ns_subdiv = [self.find_closest(i,tmp_n) for i in self.c1u_s_subdiv]
#!輪郭のメッシュによる細分化プロット（基底変換後）
plt.scatter([i[0] for i in self.c2u_ps_subdiv],[i[1] for i in self.c2u_ps_subdiv ],color =[0.0, 1.0, 0.0],s = 30,label = f"Subdivision Mesh (n={self.num_subdiv})")
plt.scatter([i[0] for i in self.c2u_ns_subdiv],[i[1] for i in self.c2u_ns_subdiv],color =[0.0, 1.0, 0.0],s = 30)
for i,j in zip(self.c2u_ps_subdiv, self.c2u_ns_subdiv):
plt.plot([i[0],j[0]],[i[1],j[1]],color =[0.0, 1.0, 0.0])
#??基底変換前の輪郭のメッシュ情報の計算
a,b,c,d = v1_x,v1_y,v2_x,v2_y
self.cx_ps_subdiv = [(d*i[0]-b*i[1])/(a*d-b*c)+self.cx for i in self.c2u_ps_subdiv]
self.cx_ns_subdiv = [(d*i[0]-b*i[1])/(a*d-b*c)+self.cx for i in self.c2u_ns_subdiv]
self.cy_ps_subdiv = [(c*i[0]-a*i[1])/(c*b-a*d)+self.cy for i in self.c2u_ps_subdiv]
self.cy_ns_subdiv = [(c*i[0]-a*i[1])/(c*b-a*d)+self.cy for i in self.c2u_ns_subdiv]
#!R2取得のための数値計算
self.yp_R2 = [[i,j] for i,j in zip(self.c1s,self.c2s) if j>=0]
self.yn_R2= [[i,j] for i,j in zip(self.c1s,self.c2s) if j< 0]
a = max([self.intersec_a1_u[0],self.intersec_a2_u[0],self.intersec_a3_u[0],self.intersec_a4_u[0]])
b = max([self.intersec_a1_u[1],self.intersec_a2_u[1],self.intersec_a3_u[1],self.intersec_a4_u[1]])
self.yp_R2_hat = [b/a*np.sqrt(a**2-i**2) for i in [j[0] for j in self.yp_R2]]
self.yn_R2_hat = [-b/a*np.sqrt(a**2-i**2) for i in [j[0] for j in self.yn_R2]]
print(self.yp_R2)
print(self.yp_R2_hat)
mean_yp_R2 = sum([i[1] for i in self.yp_R2])/len(self.yp_R2)
mean_yn_R2 = sum([i[1] for i in self.yn_R2])/len(self.yp_R2)
RSS = 0
TSS = 0
for i,j,k,l in zip(self.yp_R2_hat,self.yn_R2_hat,self.yp_R2,self.yn_R2):
print("#################################################")
RSS += (i-k[1])**2
RSS += (j-l[1])**2
TSS += (k[1] - mean_yp_R2)**2
TSS += (l[1] - mean_yp_R2)**2
print(1-RSS/TSS)
self.R2 = 1-RSS/TSS
plt.scatter(c1s, c2s, color="black", label="Cell contour",s=10)
plt.scatter(0, 0, label="Center (0,0)", color="green")
plt.scatter([self.intersec_a1_u[0], self.intersec_a2_u[0]], [self.intersec_a1_u[1], self.intersec_a2_u[1]], color="blue",label = "intersection")
plt.scatter([self.intersec_a3_u[0], self.intersec_a4_u[0]], [self.intersec_a3_u[1], self.intersec_a4_u[1]], color="blue")
plt.plot([self.intersec_a1_u[0], self.intersec_a2_u[0]], [self.intersec_a1_u[1], self.intersec_a2_u[1]], color="blue")
plt.plot([self.intersec_a3_u[0], self.intersec_a4_u[0]], [self.intersec_a3_u[1], self.intersec_a4_u[1]], color="blue")
plt.plot(self.u1_s, self.u2_ps, color="red",label = "Fit ellipse")
plt.plot(self.u1_s, self.u2_ns, color="red")
plt.grid()
plt.xlabel(r"$u_1$", fontsize=18)
plt.ylabel(r"$u_2$", fontsize=18)
plt.legend(loc='upper right', borderaxespad=0)
fig.savefig(f"base_conv_{self.n}.png", dpi=500)
def draw(self) -> None:
fig = plt.figure(figsize=[7, 7])
#!基底変換前の輪郭座標のプロット
plt.scatter(self.x, self.y, s=4, color="black", label="Cell contour")
#!基底変換前の楕円の長軸と楕円との交点のプロット
plt.plot([self.intersec_a1[0], self.intersec_a2[0]], [self.intersec_a1[1], self.intersec_a2[1]], color="blue")
plt.scatter(self.intersec_a1[0], self.intersec_a1[1],
s=15, color="blue", label="Intersection")
plt.scatter(self.intersec_a2[0],
self.intersec_a2[1], s=15, color="blue")
#!基底変換前の楕円の短軸と楕円との交点のプロット
plt.plot([self.intersec_a3[0], self.intersec_a4[0]], [self.intersec_a3[1], self.intersec_a4[1]], color="blue")
plt.scatter(self.intersec_a3[0],
self.intersec_a3[1], s=15, color="blue")
plt.scatter(self.intersec_a4[0],
self.intersec_a4[1], s=15, color="blue")
#!中心座標のプロット
plt.scatter(self.cx, self.cy, color="green",label=f"Center ({math.floor(self.cx)},{math.floor(self.cy)})", s=10)
#!基底変換前の楕円のプロット
x = np.linspace(self.intersec_a1[0], self.intersec_a2[0], 1000)
self.y_p = np.array([(self.C_2*i+self.C_4)/2 + np.sqrt((self.C_2*i+self.C_4)** 2+4*(self.C_1*i**2+self.C_3*i+self.C_5))/2 for i in x])
self.y_n = np.array([(self.C_2*i+self.C_4)/2 - np.sqrt((self.C_2*i+self.C_4)** 2+4*(self.C_1*i**2+self.C_3*i+self.C_5))/2 for i in x])
#!基底変換前の楕円へのメッシュのプロット
# plt.scatter([self.x_ps_subdiv[0],self.x_ns_subdiv[0]],[self.y_ps_subdiv[0],self.y_ns_subdiv[0]],color =[0.0, 1.0, 0.0],s = 15,label = f"Subdivision Mesh (n={self.num_subdiv})")
# for i in range(1,len(self.x_ns_subdiv)-1):
# plt.scatter([self.x_ps_subdiv[i],self.x_ns_subdiv[i]],[self.y_ps_subdiv[i],self.y_ns_subdiv[i]],color =[0.0, 1.0, 0.0],s = 15)
# plt.plot([self.x_ns_subdiv[i],self.x_ps_subdiv[i]],[self.y_ns_subdiv[i],self.y_ps_subdiv[i]],color =[0.0, 1.0, 0.0])
#!基底変換前の輪郭へのメッシュのプロット
plt.scatter([self.cx_ps_subdiv[0],self.cx_ns_subdiv[0]],[self.cy_ps_subdiv[0],self.cy_ns_subdiv[0]],color =[0.0, 1.0, 0.0],s = 15,label = f"Subdivision Mesh (n={self.num_subdiv})")
for i in range(1,len(self.cx_ns_subdiv)-1):
plt.scatter([self.cx_ps_subdiv[i],self.cx_ns_subdiv[i]],[self.cy_ps_subdiv[i],self.cy_ns_subdiv[i]],color =[0.0, 1.0, 0.0],s = 15)
plt.plot([self.cx_ns_subdiv[i],self.cx_ps_subdiv[i]],[self.cy_ns_subdiv[i],self.cy_ps_subdiv[i]],color =[0.0, 1.0, 0.0])
plt.plot(x, self.y_p, color="red", label="Fit ellipse")
plt.plot(x, self.y_n, color="red")
plt.grid()
plt.xlabel(r"$x$", fontsize=18)
plt.ylabel(r"$y$", fontsize=18)
plt.legend()
plt.xlim(min(self.x)-20, max(self.x)+20)
plt.ylim(min(self.y_n)-20, max(self.y_p)+20)
fig.savefig(f"img_fit{self.n}.png", dpi=600)
def __repr__(self) -> str:
print("#################################################")
C_1: str = f"C1 = {self.C_1}\n"
C_2: str = f"C2 = {self.C_2}\n"
C_3: str = f"C3 = {self.C_3}\n"
C_4: str = f"C4 = {self.C_4}\n"
C_5: str = f"C5 = {self.C_5}\n"
equation: str = f"{-self.C_1}x^2+{-self.C_2}xy++y^2+{-self.C_3}x+{-self.C_4}y+{-self.C_5} = 0\n"
intersection_raw: str = f"{self.intersec_a1},{self.intersec_a2,self.intersec_a3,self.intersec_a4}"
intersection_u: str = f"{self.intersec_a1_u},{self.intersec_a2_u,self.intersec_a3_u,self.intersec_a4_u}"
return f"Parameters:\n{C_1+C_2+C_3+C_4+C_5}" + f"Equation:\n{equation}" + f"Intersec raw:\n{intersection_raw}" + f"\nIntersec u:\n{intersection_u}"
def get_params(self) -> list[float]:
return [self.C_1, self.C_2, self.C_3, self.C_4, self.C_5]
def get_R2(self) -> float:
return self.R2
def get_contours_u(self):
return self.c1s, self.c2s
contours_coordinates = get_contours(f"{1}")
L = LeastSquares(contours_coordinates,1)
L._solve2()
L.base_conversion()
contours_u = L.get_contours_u()
L.draw()画像データ
