311 lines
11 KiB
Python
311 lines
11 KiB
Python
|
import numpy as np
|
|||
|
|
import SimpleITK as sitk
|
|||
|
|
from scipy.ndimage import center_of_mass
|
|||
|
|
import matplotlib.pyplot as plt
|
|||
|
|
|
|||
|
|
def azimuth_rotation(image, show_plt=False, save_plt=False, output_path=None):
|
|||
|
|
|
|||
|
|
img = sitk.ReadImage(image, sitk.sitkUInt8)
|
|||
|
|
arr_zyx = sitk.GetArrayFromImage(img) # (z, y, x)
|
|||
|
|
|
|||
|
|
max_proj = np.max(arr_zyx, axis=0) # -> (y, x)
|
|||
|
|
|
|||
|
|
binary_proj = (max_proj > 0).astype(np.uint8)
|
|||
|
|
|
|||
|
|
ys, xs = np.where(binary_proj > 0)
|
|||
|
|
if len(xs) < 10:
|
|||
|
|
raise ValueError("Not enough foreground points")
|
|||
|
|
|
|||
|
|
# 2) centroid ←←← 這裡一定會定義 cx, cy
|
|||
|
|
cy = ys.mean()
|
|||
|
|
cx = xs.mean()
|
|||
|
|
centroid = np.array([cy, cx])
|
|||
|
|
|
|||
|
|
cy = ys.mean()
|
|||
|
|
cx = xs.mean()
|
|||
|
|
centroid = np.array([cy, cx])
|
|||
|
|
|
|||
|
|
y_min = ys.min()
|
|||
|
|
top_row_mask = (ys == y_min)
|
|||
|
|
xs_top_row = xs[top_row_mask]
|
|||
|
|
|
|||
|
|
# 取這一排的中位數或平均值
|
|||
|
|
x_center = int(np.median(xs_top_row)) # 或用 np.mean()
|
|||
|
|
top_point = (y_min, x_center)
|
|||
|
|
# print(f"最上排中心點: {top_point}")
|
|||
|
|
|
|||
|
|
# 計算從 top_point 到 centroid 的向量
|
|||
|
|
dy = cy - y_min # y 方向的變化
|
|||
|
|
dx = cx - x_center # x 方向的變化
|
|||
|
|
|
|||
|
|
# 計算與 y 軸的夾角
|
|||
|
|
# 注意：影像座標系中 y 軸向下，所以要特別處理
|
|||
|
|
angle_rad = np.arctan2(dx, dy) # 弧度
|
|||
|
|
angle_deg = np.degrees(angle_rad) # 轉成角度
|
|||
|
|
|
|||
|
|
if show_plt:
|
|||
|
|
|
|||
|
|
fig = plt.figure(figsize=(12, 12))
|
|||
|
|
|
|||
|
|
# 視覺化時加上角度資訊
|
|||
|
|
plt.imshow(binary_proj, cmap='gray')
|
|||
|
|
|
|||
|
|
plt.scatter(x_center, y_min, c='red', s=60, label='Top')
|
|||
|
|
plt.scatter(cx, cy, c='yellow', s=60, label='Centroid')
|
|||
|
|
plt.plot([x_center, cx], [y_min, cy], 'c-', lw=2,
|
|||
|
|
label=f'Angle with y-axis: {angle_deg:.1f}°')
|
|||
|
|
|
|||
|
|
plt.legend()
|
|||
|
|
plt.title("Top point + centroid-directed line")
|
|||
|
|
plt.axis("off")
|
|||
|
|
|
|||
|
|
if save_plt:
|
|||
|
|
if output_path is None:
|
|||
|
|
output_path = "azimuth_rotation.png"
|
|||
|
|
fig.savefig(output_path, dpi=200, bbox_inches="tight")
|
|||
|
|
|
|||
|
|
plt.show()
|
|||
|
|
plt.close(fig)
|
|||
|
|
|
|||
|
|
return angle_deg
|
|||
|
|
|
|||
|
|
def split_spine_anterior_posterior(image_path, center_mode='com'):
|
|||
|
|
"""
|
|||
|
|
從 sagittal view 看，沿著 y 軸（前後方向）將脊椎切成前半部和後半部
|
|||
|
|
|
|||
|
|
Parameters:
|
|||
|
|
image_path (str): 影像路徑
|
|||
|
|
center_mode (str): 'com' 使用質心的 y 座標，'image' 使用圖片中心的 y 座標
|
|||
|
|
|
|||
|
|
Returns:
|
|||
|
|
anterior, posterior: 前半部和後半部的 binary mask
|
|||
|
|
"""
|
|||
|
|
|
|||
|
|
# Sagittal projection: 沿著 x 軸投影 -> (z, y)
|
|||
|
|
img = sitk.ReadImage(image_path, sitk.sitkUInt8)
|
|||
|
|
arr_zyx = sitk.GetArrayFromImage(img)
|
|||
|
|
|
|||
|
|
# Sagittal projection
|
|||
|
|
max_proj_sagittal = np.max(arr_zyx, axis=2)
|
|||
|
|
binary_proj = (max_proj_sagittal > 0).astype(np.uint8)
|
|||
|
|
|
|||
|
|
# 決定切割的 y 座標
|
|||
|
|
if center_mode == 'com':
|
|||
|
|
cz, cy = center_of_mass(binary_proj)
|
|||
|
|
split_y = int(round(cy))
|
|||
|
|
label = f'Center of Mass (y={split_y})'
|
|||
|
|
|
|||
|
|
elif center_mode == 'image':
|
|||
|
|
split_y = binary_proj.shape[1] // 2
|
|||
|
|
label = f'Image Center (y={split_y})'
|
|||
|
|
|
|||
|
|
# 檢查是否為數字，且範圍在 0 到 1 之間 (不含邊界)
|
|||
|
|
elif isinstance(center_mode, (int, float)) and 0 < center_mode < 1:
|
|||
|
|
ys = np.where(binary_proj > 0)[1]
|
|||
|
|
|
|||
|
|
if ys.size == 0: # 額外保險：如果投影是空的
|
|||
|
|
split_y = binary_proj.shape[1] // 2
|
|||
|
|
else:
|
|||
|
|
y_min, y_max = ys.min(), ys.max()
|
|||
|
|
split_y = int(round(y_min + center_mode * (y_max - y_min)))
|
|||
|
|
|
|||
|
|
label = f'Custom Ratio {center_mode} (y={split_y})'
|
|||
|
|
|
|||
|
|
else:
|
|||
|
|
raise ValueError("center_mode 必須是 'com'、'image' 或介於 0 到 1 之間的浮點數 (例如 0.8)")
|
|||
|
|
|
|||
|
|
# 切割：anterior (y < split_y) 和 posterior (y >= split_y)
|
|||
|
|
anterior = binary_proj.copy()
|
|||
|
|
posterior = binary_proj.copy()
|
|||
|
|
|
|||
|
|
anterior[:, :split_y] = 0 # 保留spine前半部(image後半部)（y >= split_y）
|
|||
|
|
posterior[:, split_y:] = 0 # 保留spine後半部(image前半部)（y < split_y）
|
|||
|
|
|
|||
|
|
return anterior, posterior, binary_proj
|
|||
|
|
|
|||
|
|
def analyze_vertebral_tilt_contour(image_path, edge_type='superior', show_plot=False, debug=False, save_plt=False, output_path=None):
|
|||
|
|
"""
|
|||
|
|
通過椎體前緣輪廓分析傾斜（可選上或下終板）
|
|||
|
|
|
|||
|
|
Parameters:
|
|||
|
|
edge_type: 'superior' 上終板, 'inferior' 下終板, 'both' 兩者都分析
|
|||
|
|
"""
|
|||
|
|
|
|||
|
|
# 切割出 anterior 部分
|
|||
|
|
anterior, posterior, binary_proj = split_spine_anterior_posterior(image_path, center_mode='com')
|
|||
|
|
|
|||
|
|
zs, ys = np.where(anterior > 0)
|
|||
|
|
|
|||
|
|
if len(zs) == 0:
|
|||
|
|
return None
|
|||
|
|
|
|||
|
|
from sklearn.linear_model import RANSACRegressor
|
|||
|
|
|
|||
|
|
results = {}
|
|||
|
|
|
|||
|
|
# === 根據 edge_type 決定要分析哪些邊 ===
|
|||
|
|
edges_to_analyze = []
|
|||
|
|
if edge_type == 'superior' or edge_type == 'both':
|
|||
|
|
edges_to_analyze.append('superior')
|
|||
|
|
if edge_type == 'inferior' or edge_type == 'both':
|
|||
|
|
edges_to_analyze.append('inferior')
|
|||
|
|
|
|||
|
|
all_edge_points = {}
|
|||
|
|
all_inliers = {}
|
|||
|
|
all_outliers = {}
|
|||
|
|
all_slopes = {}
|
|||
|
|
all_intercepts = {}
|
|||
|
|
all_angles = {}
|
|||
|
|
|
|||
|
|
for edge in edges_to_analyze:
|
|||
|
|
# 對每個 y，找對應的邊緣點
|
|||
|
|
edge_points = []
|
|||
|
|
unique_ys = np.unique(ys)
|
|||
|
|
|
|||
|
|
for y in unique_ys:
|
|||
|
|
z_at_y = zs[ys == y]
|
|||
|
|
|
|||
|
|
if edge == 'superior':
|
|||
|
|
z_edge = z_at_y.max() # 最上面的點（z 最小）
|
|||
|
|
else: # inferior
|
|||
|
|
z_edge = z_at_y.min() # 最下面的點（z 最大）
|
|||
|
|
|
|||
|
|
edge_points.append([y, z_edge])
|
|||
|
|
|
|||
|
|
edge_points = np.array(edge_points)
|
|||
|
|
all_edge_points[edge] = edge_points
|
|||
|
|
|
|||
|
|
if len(edge_points) < 10:
|
|||
|
|
continue
|
|||
|
|
|
|||
|
|
# RANSAC 擬合
|
|||
|
|
X = edge_points[:, 0].reshape(-1, 1)
|
|||
|
|
y_data = edge_points[:, 1]
|
|||
|
|
|
|||
|
|
ransac = RANSACRegressor(
|
|||
|
|
residual_threshold=5.0,
|
|||
|
|
random_state=42
|
|||
|
|
)
|
|||
|
|
ransac.fit(X, y_data)
|
|||
|
|
|
|||
|
|
inlier_mask = ransac.inlier_mask_
|
|||
|
|
outlier_mask = ~inlier_mask
|
|||
|
|
|
|||
|
|
edge_points_inliers = edge_points[inlier_mask]
|
|||
|
|
edge_points_outliers = edge_points[outlier_mask]
|
|||
|
|
|
|||
|
|
all_inliers[edge] = edge_points_inliers
|
|||
|
|
all_outliers[edge] = edge_points_outliers
|
|||
|
|
|
|||
|
|
# 獲取擬合結果
|
|||
|
|
slope = ransac.estimator_.coef_[0]
|
|||
|
|
intercept = ransac.estimator_.intercept_
|
|||
|
|
|
|||
|
|
all_slopes[edge] = slope
|
|||
|
|
all_intercepts[edge] = intercept
|
|||
|
|
|
|||
|
|
# 計算 R²
|
|||
|
|
y_pred = ransac.predict(edge_points_inliers[:, 0].reshape(-1, 1))
|
|||
|
|
ss_res = np.sum((edge_points_inliers[:, 1] - y_pred) ** 2)
|
|||
|
|
ss_tot = np.sum((edge_points_inliers[:, 1] - np.mean(edge_points_inliers[:, 1])) ** 2)
|
|||
|
|
r_squared = 1 - (ss_res / ss_tot) if ss_tot > 0 else 0
|
|||
|
|
|
|||
|
|
# 計算傾斜角度
|
|||
|
|
tilt_angle = np.degrees(np.arctan(slope))
|
|||
|
|
all_angles[edge] = tilt_angle
|
|||
|
|
|
|||
|
|
if debug:
|
|||
|
|
print(f"\n=== {edge.upper()} ENDPLATE ===")
|
|||
|
|
print(f"Total points: {len(edge_points)}")
|
|||
|
|
print(f"Inliers: {len(edge_points_inliers)}")
|
|||
|
|
print(f"Outliers: {len(edge_points_outliers)}")
|
|||
|
|
|
|||
|
|
print(f"{edge.capitalize()} 終板傾斜角度: {tilt_angle:.2f}°")
|
|||
|
|
print(f"斜率: {slope:.4f}, R²: {r_squared:.4f}")
|
|||
|
|
|
|||
|
|
results[edge] = {
|
|||
|
|
'tilt_angle_deg': tilt_angle,
|
|||
|
|
'slope': slope,
|
|||
|
|
'intercept': intercept,
|
|||
|
|
'r_squared': r_squared,
|
|||
|
|
'n_inliers': len(edge_points_inliers),
|
|||
|
|
'n_outliers': len(edge_points_outliers)
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
# Visualization
|
|||
|
|
if show_plot:
|
|||
|
|
n_edges = len(edges_to_analyze)
|
|||
|
|
fig, axes = plt.subplots(n_edges, 2, figsize=(16, 6*n_edges))
|
|||
|
|
|
|||
|
|
if n_edges == 1:
|
|||
|
|
axes = axes.reshape(1, -1)
|
|||
|
|
|
|||
|
|
colors = {'superior': 'red', 'inferior': 'cyan'}
|
|||
|
|
|
|||
|
|
for idx, edge in enumerate(edges_to_analyze):
|
|||
|
|
edge_points = all_edge_points[edge]
|
|||
|
|
edge_points_inliers = all_inliers[edge]
|
|||
|
|
edge_points_outliers = all_outliers[edge]
|
|||
|
|
slope = all_slopes[edge]
|
|||
|
|
intercept = all_intercepts[edge]
|
|||
|
|
tilt_angle = all_angles[edge]
|
|||
|
|
color = colors[edge]
|
|||
|
|
|
|||
|
|
# 左圖：scatter plot
|
|||
|
|
ax_left = axes[idx, 0]
|
|||
|
|
|
|||
|
|
if len(edge_points_outliers) > 0:
|
|||
|
|
ax_left.scatter(edge_points_outliers[:, 0], edge_points_outliers[:, 1],
|
|||
|
|
c='lightcoral', s=30, alpha=0.6, marker='x',
|
|||
|
|
label=f'Outliers ({len(edge_points_outliers)})', zorder=3)
|
|||
|
|
|
|||
|
|
ax_left.scatter(edge_points_inliers[:, 0], edge_points_inliers[:, 1],
|
|||
|
|
c=color, s=20, alpha=0.7,
|
|||
|
|
label=f'Inliers ({len(edge_points_inliers)})', zorder=4)
|
|||
|
|
|
|||
|
|
# 擬合線
|
|||
|
|
y_line = np.array([edge_points[:, 0].min(), edge_points[:, 0].max()])
|
|||
|
|
z_line = slope * y_line + intercept
|
|||
|
|
ax_left.plot(y_line, z_line, 'lime', linewidth=3,
|
|||
|
|
label=f'Angle: {tilt_angle:.1f}°\nR²: {results[edge]["r_squared"]:.3f}',
|
|||
|
|
zorder=5)
|
|||
|
|
|
|||
|
|
ax_left.set_xlabel('y')
|
|||
|
|
ax_left.set_ylabel('z')
|
|||
|
|
ax_left.set_title(f'{edge.capitalize()} Endplate Analysis\nTilt: {tilt_angle:.1f}°')
|
|||
|
|
ax_left.invert_yaxis()
|
|||
|
|
ax_left.legend()
|
|||
|
|
ax_left.grid(True, alpha=0.3)
|
|||
|
|
|
|||
|
|
# 右圖：原始影像
|
|||
|
|
ax_right = axes[idx, 1]
|
|||
|
|
ax_right.imshow(anterior, cmap='gray', aspect='equal')
|
|||
|
|
|
|||
|
|
if len(edge_points_outliers) > 0:
|
|||
|
|
ax_right.scatter(edge_points_outliers[:, 0], edge_points_outliers[:, 1],
|
|||
|
|
c='red', s=40, alpha=0.7, marker='x', label='Outliers', zorder=4)
|
|||
|
|
|
|||
|
|
ax_right.scatter(edge_points_inliers[:, 0], edge_points_inliers[:, 1],
|
|||
|
|
c=color, s=25, alpha=0.8, label='Inliers', zorder=3)
|
|||
|
|
|
|||
|
|
ax_right.plot(y_line, z_line, 'lime', linewidth=3, linestyle='--',
|
|||
|
|
label=f'{edge.capitalize()}: {tilt_angle:.1f}°', zorder=5)
|
|||
|
|
|
|||
|
|
ax_right.set_title(f'Anterior Half - {edge.capitalize()} Edge')
|
|||
|
|
ax_right.set_xlabel('y')
|
|||
|
|
ax_right.set_ylabel('z')
|
|||
|
|
ax_right.invert_yaxis()
|
|||
|
|
ax_right.legend()
|
|||
|
|
|
|||
|
|
plt.tight_layout()
|
|||
|
|
|
|||
|
|
if save_plt:
|
|||
|
|
if output_path is None:
|
|||
|
|
output_path = "analyze_vertebral_tilt_contour.png"
|
|||
|
|
fig.savefig(output_path, dpi=200, bbox_inches="tight")
|
|||
|
|
|
|||
|
|
plt.show()
|
|||
|
|
plt.close(fig)
|
|||
|
|
|
|||
|
|
return results
|