Newer
Older
Demo-Maker / stethoscope_analysis.py
import cv2
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd


def calculate_welzl_circle(points: np.ndarray) -> tuple[tuple[float, float], float]:
    """
    Welzlのアルゴリズムを使用して最小包含円を計算
    Args:
        points: shape=(n, 2) のnumpy配列
    Returns:
        ((center_x, center_y), radius)
    """

    def make_circle_3points(
        p1: np.ndarray, p2: np.ndarray, p3: np.ndarray
    ) -> tuple[np.ndarray, float]:
        """3点から円を計算"""
        # 3点から円の中心と半径を計算
        temp = p2[0] * p2[0] + p2[1] * p2[1]
        bc = (p1[0] * p1[0] + p1[1] * p1[1] - temp) / 2.0
        cd = (temp - p3[0] * p3[0] - p3[1] * p3[1]) / 2.0
        det = (p1[0] - p2[0]) * (p2[1] - p3[1]) - (p2[0] - p3[0]) * (p1[1] - p2[1])

        if abs(det) < 1e-10:
            return None

        center_x = (bc * (p2[1] - p3[1]) - cd * (p1[1] - p2[1])) / det
        center_y = ((p1[0] - p2[0]) * cd - (p2[0] - p3[0]) * bc) / det
        center = np.array([center_x, center_y])
        radius = np.sqrt(np.sum((p1 - center) ** 2))

        return center, radius

    def make_circle_2points(p1: np.ndarray, p2: np.ndarray) -> tuple[np.ndarray, float]:
        """2点から円を計算(直径を使用)"""
        center = (p1 + p2) / 2
        radius = np.sqrt(np.sum((p1 - p2) ** 2)) / 2
        return center, radius

    def is_point_in_circle(
        point: np.ndarray, center: np.ndarray, radius: float
    ) -> bool:
        """点が円の中にあるかチェック"""
        return np.sum((point - center) ** 2) <= radius * radius * (1 + 1e-10)

    def minimal_enclosing_circle(points: np.ndarray) -> tuple[np.ndarray, float]:
        """最小包含円を再帰的に計算"""
        if len(points) <= 1:
            if len(points) == 1:
                return points[0], 0
            return None

        # 2点の場合
        if len(points) == 2:
            return make_circle_2points(points[0], points[1])

        # ランダムに2点を選び、それらを通る円を作成
        for i in range(3):
            for j in range(i + 1, len(points)):
                center, radius = make_circle_2points(points[i], points[j])
                all_points_in_circle = all(
                    is_point_in_circle(p, center, radius) for p in points
                )

                if all_points_in_circle:
                    return center, radius

        # 3点を使用して円を作成
        for i in range(len(points) - 2):
            for j in range(i + 1, len(points) - 1):
                for k in range(j + 1, len(points)):
                    result = make_circle_3points(points[i], points[j], points[k])
                    if result is None:
                        continue
                    center, radius = result
                    if all(is_point_in_circle(p, center, radius) for p in points):
                        return center, radius

        raise ValueError("No valid circle found")

    # 最小包含円の計算
    points = points.copy()
    np.random.shuffle(points)  # ランダム化により計算を効率化
    center, radius = minimal_enclosing_circle(points)

    return tuple(center), radius


def calculate_metrics_for_sequence(
    points: np.ndarray,
) -> tuple[float, float, float, tuple[float, float]]:
    """
    時系列データの点列から標準偏差と最小包含円の半径を計算
    Args:
        points: shape=(n, 2) のnumpy配列
    Returns:
        (std_x, std_y, min_circle_radius, (center_x, center_y))
    """
    # 標準偏差の計算
    std_x = np.std(points[:, 0])
    std_y = np.std(points[:, 1])

    # 最小包含円の計算
    center, radius = calculate_welzl_circle(points)

    return std_x, std_y, radius, center


def draw_circles_on_body(body_image: np.ndarray, circles_info: dict) -> np.ndarray:
    """
    体の画像上に最小包含円を描画
    Args:
        body_image: 体の画像
        circles_info: 各手法の円の情報 {method_name: (center_x, center_y, radius)}
    Returns:
        描画された画像
    """
    # 色の定義
    colors = {
        "stethoscope": (0, 255, 255),  # 黄色
        "conv": (0, 255, 0),  # 緑
        "xgboost": (255, 0, 0),  # 青
        "lightgbm": (0, 0, 255),  # 赤
    }

    result_image = body_image.copy()

    # 各手法の円を描画
    for method, (center_x, center_y, radius) in circles_info.items():
        color = colors.get(method, (255, 255, 255))  # デフォルトは白
        # 円を描画
        cv2.circle(result_image, (int(center_x), int(center_y)), int(radius), color, 2)
        # 中心点を描画
        cv2.circle(result_image, (int(center_x), int(center_y)), 3, color, -1)
        # 手法名を描画
        cv2.putText(
            result_image,
            method,
            (int(center_x), int(center_y - radius - 10)),
            cv2.FONT_HERSHEY_SIMPLEX,
            0.5,
            color,
            2,
        )

    return result_image


def main():
    # CSVの読み込み
    df = pd.read_csv("output/results-rtmpose+yolox/results.csv")

    # 体の画像を読み込み
    body_image = cv2.imread("images/body/BodyF.png")
    if body_image is None:
        raise FileNotFoundError("Body image not found")

    # 各手法のカラム名のペア
    method_columns = {
        "conv": ["conv_stethoscope_x", "conv_stethoscope_y"],
    }

    # XGBoostとLightGBMのカラムが存在する場合は追加
    if "Xgboost_stethoscope_x" in df.columns:
        method_columns["xgboost"] = ["Xgboost_stethoscope_x", "Xgboost_stethoscope_y"]
    if "lightGBM_stethoscope_x" in df.columns:
        method_columns["lightgbm"] = [
            "lightGBM_stethoscope_x",
            "lightGBM_stethoscope_y",
        ]

    # 結果を格納する辞書
    results = {}
    circles_info = {}

    # 各手法について計算
    for method_name, (x_col, y_col) in method_columns.items():
        # 有効な点のみを抽出 (0や欠損値を除外)
        valid_points = df[[x_col, y_col]].values
        valid_mask = ~(
            np.isnan(valid_points).any(axis=1) | (valid_points == 0).all(axis=1)
        )
        valid_points = valid_points[valid_mask]

        if len(valid_points) > 0:
            std_x, std_y, radius, center = calculate_metrics_for_sequence(valid_points)
            results[method_name] = {
                "points_count": len(valid_points),
                "std_x": std_x,
                "std_y": std_y,
                "min_circle_radius": radius,
            }
            circles_info[method_name] = (center[0], center[1], radius)

    # 結果をDataFrameに変換
    results_df = pd.DataFrame(results).T

    # 結果の表示
    print("\n=== 手法別分析結果 ===")
    print(results_df.round(3))

    # 結果をCSVに保存
    results_df.to_csv("stethoscope_analysis_results.csv")

    # 体の画像上に円を描画
    result_image = draw_circles_on_body(body_image, circles_info)

    # 結果の画像を保存
    cv2.imwrite("stethoscope_circles_on_body.png", result_image)

    # オリジナルの数値データの可視化も保持
    plt.figure(figsize=(15, 10))

    # サブプロット1: STDの比較
    plt.subplot(2, 1, 1)
    x = np.arange(len(results_df))
    width = 0.35
    plt.bar(
        x - width / 2, results_df["std_x"], width, label="STD X", color="red", alpha=0.7
    )
    plt.bar(
        x + width / 2,
        results_df["std_y"],
        width,
        label="STD Y",
        color="blue",
        alpha=0.7,
    )
    plt.xticks(x, results_df.index, rotation=45)
    plt.title("Standard Deviation by Method")
    plt.ylabel("Value (pixels)")
    plt.legend()
    plt.grid(True)

    # サブプロット2: 最小包含円の半径
    plt.subplot(2, 1, 2)
    plt.bar(results_df.index, results_df["min_circle_radius"], color="green", alpha=0.7)
    plt.title("Minimum Enclosing Circle Radius by Method")
    plt.ylabel("Radius (pixels)")
    plt.xticks(rotation=45)
    plt.grid(True)

    plt.tight_layout()
    plt.savefig("stethoscope_metrics_comparison.png")
    plt.close()

    # 詳細な統計情報の出力
    print("\n=== 詳細統計情報 ===")
    for method_name in results:
        print(f"\n{method_name}:")
        circle_info = circles_info[method_name]
        print(f"データ点数: {results[method_name]['points_count']}")
        print(f"中心座標: ({circle_info[0]:.2f}, {circle_info[1]:.2f})")
        print(f"X方向の標準偏差: {results[method_name]['std_x']:.2f} pixels")
        print(f"Y方向の標準偏差: {results[method_name]['std_y']:.2f} pixels")
        print(
            f"最小包含円の半径: {results[method_name]['min_circle_radius']:.2f} pixels"
        )


if __name__ == "__main__":
    main()