import numpy as np
import SimpleITK as sitk
from scipy.ndimage import center_of_mass, rotate
import matplotlib.pyplot as plt

def _best_vertical_split(mask2d):
    """
    binary 2D 陣列（最後一軸 = x）：找讓 mask 與其鏡射重疊最大的垂直線 x = t。
    c[s] = Σ_u A[u]·A[s-u] 是每列自卷積；批次 FFT 一次算出所有 s。
    回傳 (score, s_best)；s = 2·t（t 可能為半整數）。
    """
    n = mask2d.shape[-1]
    flat = np.asarray(mask2d).reshape(-1, n)
    rows = flat[flat.max(axis=1) > 0]
    if rows.size == 0:
        return 0.0, 0
    L = 1 << (2 * n - 1).bit_length()
    padded = np.zeros((rows.shape[0], L), dtype=np.float32)
    padded[:, :n] = rows.astype(np.float32)
    F = np.fft.rfft(padded, axis=1)
    conv = np.fft.irfft(F * F, axis=1)[:, :2 * n - 1]
    score = np.clip(conv.sum(axis=0), 0, None)
    s_best = int(np.argmax(score))
    return float(score[s_best]), s_best


def best_symmetry_axis_angle(proj, coarse_step=10.0, fine_step=1.0):
    """
    2D binary 投影（axis0=row y, axis1=col x）的左右鏡稱軸搜尋：
    將影像旋轉 θ 後，鏡稱軸變成垂直線（col = const），用 _best_vertical_split 評分；
    粗搜 0..180°（coarse_step）再在 winner 附近細搜（fine_step）。
    回傳 (theta_deg, score, s_best)。
    """
    base = (np.asarray(proj) > 0).astype(np.float32)
    best_th, best, best_s = 0.0, -1.0, 0
    for th in np.arange(0.0, 180.0, coarse_step):
        sc, s = _best_vertical_split(rotate(base, th, reshape=False, order=0) > 0.5)
        if sc > best:
            best_th, best, best_s = float(th), sc, s
    for th in np.arange(best_th - coarse_step, best_th + coarse_step, fine_step):
        th2 = float(th) % 180.0
        sc, s = _best_vertical_split(rotate(base, th2, reshape=False, order=0) > 0.5)
        if sc > best:
            best_th, best, best_s = th2, sc, s
    return best_th, best, best_s


def _mirror_hit_count(X, Y, Z, n, d, m, sz, sy, sx):
    """每個骨 voxel 對平面 n·p=d 鏡射後四捨五入到最近 voxel，回傳命中骨 voxel 的數量"""
    nx, ny, nz = n
    dist = X * nx + Y * ny + Z * nz - d
    rx = (X - 2.0 * dist * nx).round().astype(np.int32)
    ry = (Y - 2.0 * dist * ny).round().astype(np.int32)
    rz = (Z - 2.0 * dist * nz).round().astype(np.int32)
    ok = (rx >= 0) & (rx < sx) & (ry >= 0) & (ry < sy) & (rz >= 0) & (rz < sz)
    if not ok.any():
        return 0
    return int(m[rz[ok], ry[ok], rx[ok]].sum())


def best_symmetry_plane(mask_zyx, phi_max=45.0, subsample=7,
                        coarse_step=7.5, fine_step=1.0):
    """
    3D bone mask (z, y, x) 的最佳鏡稱面，一般平面方程 a·x + b·y + c·z = d
    （voxel index 座標；(a,b,c) 為單位法線，方向任意，不限制平行 YZ 面）。
    以「每個 bone voxel 鏡射後四捨五入到最近 voxel 的命中率」為分數，
    參數化 n = R_y(phi)·R_z(theta)·(1,0,0)：
      theta : 法線在 axial (x,y) 平面內的旋轉
      phi   : 法線出平面的傾斜，限制在 [-phi_max, +phi_max] 以確保仍切分左右
    搜尋：theta 由 2D axial 投影粗定位，再 3D 三段式（粗→細→微細）
    （粗/細階段用等距子樣 voxel 加速）。
    回傳 dict：
      plane     : (a, b, c, d)  ->  a*x + b*y + c*z = d
      normal    : (a, b, c)
      offset    : d
      theta_deg / phi_deg : 法線參數
      ratio     : 鏡射命中率（0..1）
      u, v      : 平面內兩個正交方向（供繪製用）
    """
    m = np.asarray(mask_zyx) > 0
    sz, sy, sx = m.shape
    zz, yy, xx = np.nonzero(m)
    if xx.size < 100:
        n = np.array([1.0, 0.0, 0.0])
        d = (sx - 1) / 2.0
        return {'plane': (1.0, 0.0, 0.0, float(d)), 'normal': (1.0, 0.0, 0.0),
                'offset': float(d), 'theta_deg': 0.0, 'phi_deg': 0.0,
                'ratio': 1.0, 'u': (0.0, 1.0, 0.0), 'v': (0.0, 0.0, 1.0)}
    xf = xx.astype(np.float32)
    yf = yy.astype(np.float32)
    zf = zz.astype(np.float32)
    N = xf.size
    c0 = np.array([(sx - 1) / 2.0, (sy - 1) / 2.0, (sz - 1) / 2.0])
    Xs, Ys, Zs = xf[::subsample], yf[::subsample], zf[::subsample]

    def n_of(theta_deg, phi_deg):
        t = np.deg2rad(theta_deg)
        p = np.deg2rad(phi_deg)
        return np.array([np.cos(t) * np.cos(p), np.sin(t), -np.cos(t) * np.sin(p)])

    # theta 初值：2D axial 投影搜尋（垂直面情形）
    proj = m.max(axis=0)
    theta0 = best_symmetry_axis_angle(proj.astype(np.float32))[0] if proj.sum() >= 10 else 0.0

    best = None

    def consider(theta, phi, d, full=False):
        nonlocal best
        n = n_of(theta, phi)
        X, Y, Z = (xf, yf, zf) if full else (Xs, Ys, Zs)
        s = _mirror_hit_count(X, Y, Z, n, d, m, sz, sy, sx)
        if best is None or s > best['score']:
            best = {'score': s, 'theta': float(theta), 'phi': float(phi), 'd': float(d)}

    # stage 1：粗搜尋（子樣）
    for th in np.arange(theta0 - coarse_step * 2, theta0 + coarse_step * 2 + 1e-9, coarse_step):
        for ph in np.arange(-phi_max, phi_max + 1e-9, coarse_step):
            n = n_of(th, ph)
            d0 = float(n @ c0)
            for dd in (-10.0, 0.0, 10.0):
                consider(th, ph, d0 + dd)
    # stage 2：細搜尋（子樣）
    for th in np.arange(best['theta'] - 2 * fine_step * 2, best['theta'] + 2 * fine_step * 2 + 1e-9, fine_step):
        for ph in np.arange(best['phi'] - 2 * fine_step * 2, best['phi'] + 2 * fine_step * 2 + 1e-9, fine_step):
            for dd in (-3.0, -1.0, 0.0, 1.0, 3.0):
                consider(th, ph, best['d'] + dd)
    # stage 3：微細搜尋（全分辨率）
    for th in np.arange(best['theta'] - fine_step, best['theta'] + fine_step + 1e-9, fine_step / 2):
        for ph in np.arange(best['phi'] - fine_step, best['phi'] + fine_step + 1e-9, fine_step / 2):
            for dd in (-1.0, -0.5, 0.0, 0.5, 1.0):
                consider(th, ph, best['d'] + dd, full=True)
    # stage 4：d 微調（全分辨率）
    b = best
    for dd in np.arange(-1.0, 1.0 + 1e-9, 0.25):
        consider(b['theta'], b['phi'], b['d'] + dd, full=True)

    n = n_of(best['theta'], best['phi'])
    theta = best['theta']
    d = best['d']
    if theta > 90.0 or theta < -90.0:     # 正規化到 (-90, 90]，法線翻轉時 d 取反
        theta -= 180.0
        n = -n
        d = -d
    ratio = _mirror_hit_count(xf, yf, zf, n, d, m, sz, sy, sx) / N
    u = np.cross(n, [0.0, 0.0, 1.0])
    if np.linalg.norm(u) < 0.1:
        u = np.cross(n, [1.0, 0.0, 0.0])
    u = u / np.linalg.norm(u)
    v = np.cross(n, u)
    return {
        'plane': (float(n[0]), float(n[1]), float(n[2]), float(d)),
        'normal': (float(n[0]), float(n[1]), float(n[2])),
        'offset': float(d),
        'theta_deg': float(theta),
        'phi_deg': float(best['phi']),
        'ratio': float(ratio),
        'u': (float(u[0]), float(u[1]), float(u[2])),
        'v': (float(v[0]), float(v[1]), float(v[2])),
    }

def segment_spinous_process(mask_zyx, sym, band_frac=0.06, min_band=6.0,
                            min_mass_frac=0.05):
    """
    以 best_symmetry_plane 的結果 sym 從 3D bone mask (z, y, x) 切出棘突。
    棘突是中線後側構造，利用鏡稱面 a·x+b·y+c·z=d 定義：
      1) 中線帶：骨 voxel 到平面的有號距離 |s| <= w，
         w = max(min_band, band_frac * s 全寬)
      2) 前後方向：平面內兩軸 (u, v) 中 |y| 分量大者，
         正規化成 +AP = 後側（本資料系 y 往前遞增，後側 = y 小側）
      3) 中線帶的 AP 分佈呈兩大叢（椎體在前、椎弓/棘突在後），
         以兩叢間的 AP 谷底為界，AP >= 谷底 的中線帶 voxel = 棘突（含中線椎弓）；
         無明顯谷底（如骨橋）fallback 取中線帶後側 15%。
    回傳 (sp_mask (z,y,x) bool, ap_thresh, info dict)；
    資料過少時 sp_mask = None（info['mode'] 說明原因）。
    """
    m = np.asarray(mask_zyx) > 0
    zz, yy, xx = np.nonzero(m)
    info = {'n_bone': int(zz.size), 'n_sp': 0, 'ap_thresh': None,
            'band_w': None, 'mode': 'empty'}
    if zz.size < 50:
        return None, None, info
    a, b, c, d = sym['plane']
    X = xx.astype(np.float64)
    Y = yy.astype(np.float64)
    Z = zz.astype(np.float64)
    s = X * a + Y * b + Z * c - d
    u = np.array(sym['u'])
    v = np.array(sym['v'])
    u_ap = u if abs(u[1]) >= abs(v[1]) else v
    if u_ap[1] > 0:
        u_ap = -u_ap                                  # +AP = 後側（y 小側）
    w = max(float(min_band), float(band_frac) * float(s.max() - s.min()))
    mid = np.abs(s) <= w
    ap = X * u_ap[0] + Y * u_ap[1] + Z * u_ap[2]
    aps = ap[mid]
    info['band_w'] = float(w)
    if aps.size < 50:
        info['mode'] = 'too_few_midline'
        return None, None, info
    lo = int(np.floor(aps.min()))
    hi = int(np.ceil(aps.max()))
    th = None
    mode = 'fallback'
    if hi - lo >= 10:
        hist, edges = np.histogram(aps, bins=range(lo, hi + 1))
        csum = np.concatenate([[0], np.cumsum(hist)])
        total = csum[-1]
        peak = hist.max()
        best_i, best_score = None, -1.0
        for i in range(len(hist)):
            if hist[i] >= 0.05 * peak:
                continue
            if csum[i] < min_mass_frac * total or (total - csum[i + 1]) < min_mass_frac * total:
                continue
            score = min(csum[i], total - csum[i + 1])
            if score > best_score:
                best_score, best_i = score, i
        if best_i is not None:
            th = float(0.5 * (edges[best_i] + edges[best_i + 1]))
            mode = 'gap'
    if th is None:
        th = float(np.quantile(aps, 0.85))
    sp_mask = np.zeros(m.shape, dtype=bool)
    sel = mid & (ap >= th)
    sp_mask[zz[sel], yy[sel], xx[sel]] = True
    info.update(n_sp=int(sel.sum()), ap_thresh=th, mode=mode)
    return sp_mask, th, info

def best_upper_endplate_plane(mask_zyx, angle_max=45.0, thresh=3.0,
                              n_iter=500, seed=42):
    """
    3D bone mask (z, y, x) 的最佳「上終板」近似平面 a·x + b·y + c·z = d
    （voxel index 座標；(a,b,c) 為朝上的單位法線）。
    1) 每個 (y, x) 欄位取最上方 bone voxel 作為頂面點（僅前側半邊，
       y >= COM_y，避開後方元素，與 2D superior endplate 定義一致）
    2) RANSAC 三點擬平面：法線限制在與 +z 軸 ≤ angle_max° 內，
       計數 ±thresh voxel 內的頂面點為 inlier，取 inlier 最多者
    3) SVD 最小二乘微調
    回傳 dict（資料不足時回傳 None）：
      plane : (a, b, c, d)
      normal / offset
      tilt_deg     : 法線與 +z 軸的夾角（上終板傾斜）
      inlier_ratio : 頂面點落在平面 ±thresh 的比例
      n_points / n_inliers
      u, v : 平面內正交方向（供繪製用）
    """
    m = np.asarray(mask_zyx) > 0
    nz, ny, nx = m.shape
    idx = np.where(m, np.arange(nz)[:, None, None], -1)
    ztop = idx.max(axis=0)                                # (y, x) 每欄最上 z
    y_split = int(round(center_of_mass(m)[1])) if m.sum() else 0
    # ztop 是 (y, x)：條件作用在 y 軸（axis 0）
    sel = (ztop >= 0) & (np.arange(ny)[:, None] >= y_split)
    yy, xx = np.where(sel)
    if xx.size < 8:
        return None
    P = np.stack((xx, yy, ztop[yy, xx]), axis=1).astype(np.float64)
    n = P.shape[0]

    rng = np.random.default_rng(seed)
    cos_min = np.cos(np.deg2rad(angle_max))
    best_cnt, best_nv, best_d = -1, None, 0.0
    for _ in range(n_iter):
        i, j, k = rng.choice(n, 3, replace=False)
        cr = np.cross(P[j] - P[i], P[k] - P[i])
        ln = np.linalg.norm(cr)
        if ln < 1e-6:
            continue
        nv = cr / ln
        if nv[2] < 0:
            nv = -nv
        if nv[2] < cos_min:                               # 法線必須朝上
            continue
        d = float(nv @ P[i])
        cnt = int(np.count_nonzero(np.abs(P @ nv - d) <= thresh))
        if cnt > best_cnt:
            best_cnt, best_nv, best_d = cnt, nv, d
    if best_nv is None:
        return None
    # SVD 微調
    inl = P[np.abs(P @ best_nv - best_d) <= thresh]
    if inl.shape[0] < 3:
        return None
    mean = inl.mean(axis=0)
    _, _, Vt = np.linalg.svd(inl - mean, full_matrices=False)
    nv = Vt[2]
    if nv[2] < 0:
        nv = -nv
    d = float(nv @ mean)
    dist = np.abs(P @ nv - d)
    inl = P[dist <= thresh]
    u = np.cross(nv, [1.0, 0.0, 0.0])
    u = u / np.linalg.norm(u)
    v = np.cross(nv, u)
    return {
        'plane': (float(nv[0]), float(nv[1]), float(nv[2]), d),
        'normal': (float(nv[0]), float(nv[1]), float(nv[2])),
        'offset': d,
        'tilt_deg': float(np.degrees(np.arccos(np.clip(nv[2], -1.0, 1.0)))),
        'inlier_ratio': float(inl.shape[0] / n),
        'n_points': int(n),
        'n_inliers': int(inl.shape[0]),
        'u': (float(u[0]), float(u[1]), float(u[2])),
        'v': (float(v[0]), float(v[1]), float(v[2])),
    }

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}°')

        # 中矢狀面：最佳鏡稱面 a·x + b·y + c·z = d（法線方向任意，不平行 YZ 面）
        # 畫該平面在體積中央 z 切片上的截線
        sym3 = best_symmetry_plane(arr_zyx)
        a3, b3, c3, d3 = sym3['plane']
        zc_ = (arr_zyx.shape[0] - 1) / 2.0
        rhs_ = d3 - c3 * zc_                    # a*x + b*y = rhs
        n2d = np.hypot(a3, b3)
        x0_ = a3 * rhs_ / (n2d * n2d)
        y0_ = b3 * rhs_ / (n2d * n2d)
        u2_ = np.array([b3, -a3]) / n2d
        tfg_ = (xs - x0_) * u2_[0] + (ys - y0_) * u2_[1]
        L_ = float(np.abs(tfg_).max())
        plt.plot([x0_ - L_ * u2_[0], x0_ + L_ * u2_[0]],
                 [y0_ - L_ * u2_[1], y0_ + L_ * u2_[1]],
                 color='magenta', lw=2, linestyle='--',
                 label=f'Mirror plane: {a3:+.2f}x {b3:+.2f}y {c3:+.2f}z = {d3:.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