{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.11.13"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":97984,"databundleVersionId":14096757,"sourceType":"competition"}],"dockerImageVersionId":31153,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# ECG Digitization – Robust Centerline Extraction v4.1\n### PhysioNet 2025 – Digitization of ECG Images  \n**Author:** Kei Itano  \n**Version:** v4.1 (Refined COG + Median + Multi-stage smoothing)\n\nThis notebook presents a robust pipeline for extracting clean ECG waveforms\nfrom pre-segmented mask images provided in the PhysioNet ECG Digitization competition.\n\nThe key goal is to create a stable, noise-resistant centerline extraction method\nthat performs consistently across 12-lead ECG images.","metadata":{}},{"cell_type":"markdown","source":"## 1. Motivation\n\nExtracting ECG waveforms from scanned images is challenging due to:\n\n- Variable scan quality  \n- Different line thickness  \n- Jagged / broken line segments  \n- OCR/segmentation noise  \n- Lead-to-lead drift  \n\nThis notebook refines the traditional COG (center of gravity) approach\nusing a hybrid strategy and multi-stage denoising pipeline.\n\n---\n\n## 2. Proposed Pipeline\n\n1. Hybrid centerline extraction  \n   - Weighted COG  \n   - Median peak estimator  \n2. NaN interpolation  \n3. Median filter smoothing  \n4. Gaussian smoothing  \n5. Savitzky–Golay final smoothing  \n6. Pixel → time/mV conversion  \n7. Optional: cross-lead alignment preview  \n\nThis method corresponds to **v4.1**, which currently produces  \nthe most stable waveforms in my experiments.\n","metadata":{}},{"cell_type":"markdown","source":"## 3. Imports & Utility Functions\n","metadata":{}},{"cell_type":"code","source":"import os\nimport cv2\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\nfrom skimage.morphology import skeletonize","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-10T05:03:40.997517Z","iopub.execute_input":"2025-12-10T05:03:40.998346Z","iopub.status.idle":"2025-12-10T05:03:42.042523Z","shell.execute_reply.started":"2025-12-10T05:03:40.998308Z","shell.execute_reply":"2025-12-10T05:03:42.041538Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 4. Centerline Extraction Versions (v1–v4)\n\nIn this section, I keep all historical versions for clarity.\n\n- **v1** – raw COG  \n- **v2** – COG + median peak  \n- **v3** – noise reduction baseline  \n- **v4 (current)** – multi-stage smoothing + hybrid peaks  \n","metadata":{}},{"cell_type":"code","source":"def extract_centerline_from_mask(mask, grid_dx, grid_dy, scale_time, scale_mv):\n    \"\"\"二値マスクから中心線を抽出し、時系列[mV]に変換\"\"\"\n    H, W = mask.shape\n    m = (mask > 0).astype(np.float32)\n    y_idx = np.arange(H, dtype=np.float32)[:, None]\n    col_sum = m.sum(axis=0) + 1e-6\n    y_cog = (m * y_idx).sum(axis=0) / col_sum\n    y_smooth = pd.Series(y_cog).rolling(5, center=True, min_periods=1).mean().values\n\n    # ---- ピクセル → 時間/電圧 ----\n    baseline_row = H / 2.0\n    time_sec   = np.arange(W) / (grid_dx * (scale_time / 25.0))\n    voltage_mV = ((baseline_row - y_smooth) / grid_dy) * (1.0 / (scale_mv / 10.0))\n\n    # ---- DataFrame整形 ----\n    df_wave = pd.DataFrame({\"time_sec\": time_sec, \"voltage_mV\": voltage_mV})\n    if \"voltage_mV\" in df_wave.columns:\n        df_wave[\"voltage_mV\"] = df_wave[\"voltage_mV\"].interpolate(method=\"linear\")\n\n    return df_wave\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 5. Proposed v4.1 Centerline Extraction\n\nThis is the main function used for the final dataset creation.\nImprovements over v4.0:\n- corrected scaling  \n- more stable gaussian sigma  \n- better small-gap interpolation  \n","metadata":{}},{"cell_type":"code","source":"#  v4 (current) – multi-stage smoothing + hybrid peaks\ndef extract_centerline_v41(mask, grid_dx, grid_dy, scale_time, scale_mv):\n    \"\"\"\n    ECG centerline extraction v4.1\n    - v1 の安定性 + v4 の改良をミックス\n    - baseline 補正の強さを適正化\n    - 過平滑を防ぐ\n    \"\"\"\n\n    import numpy as np\n    import pandas as pd\n    from scipy.ndimage import median_filter, gaussian_filter1d\n    from scipy.signal import savgol_filter\n\n    H, W = mask.shape\n    m = (mask > 0).astype(np.float32)\n    y_idx = np.arange(H, dtype=np.float32)[:, None]\n\n    col_sum = m.sum(axis=0)\n    y_cog = np.zeros(W)\n    y_peak = np.zeros(W)\n\n    for x in range(W):\n        col = m[:, x]\n\n        # ---- 1. 有効な点がない列（欠損） ----\n        if col_sum[x] < 1e-3:\n            y_cog[x] = np.nan\n            y_peak[x] = np.nan\n            continue\n\n        # ---- 2. 重心（COG） ----\n        y_cog[x] = (col * y_idx[:, 0]).sum() / col_sum[x]\n\n        # ---- 3. 非ゼロ座標の中央値ピーク ----\n        nonzero_y = np.where(col > 0)[0]\n        if len(nonzero_y) > 0:\n            y_peak[x] = np.median(nonzero_y)\n        else:\n            y_peak[x] = y_cog[x]\n\n    # ---- 4. COG + median peak のハイブリッド ----\n    y_mix = np.nanmean(np.vstack([y_cog, y_peak]), axis=0)\n\n    # ---- 5. 線形補完で欠損を埋める ----\n    y_interp = pd.Series(y_mix).interpolate(limit_direction=\"both\").values\n\n    # ---- 6. 微小スパイク除去（v4 は size=5→過剰なので 3 に戻す）----\n    y_denoised = median_filter(y_interp, size=3)\n\n    # ---- 7. ガウス平滑（適度に）----\n    y_gauss = gaussian_filter1d(y_denoised, sigma=2)\n\n    # ---- 8. Savitzky–Golay（window=81 → 過平滑を防ぐ）----\n    y_smooth = savgol_filter(y_gauss, window_length=81, polyorder=3)\n\n    # ---- 9. ベースライン補正 ----\n    baseline = savgol_filter(y_smooth, window_length=401, polyorder=1)  # v4の301→弱めに401\n    y_corrected = y_smooth - baseline\n\n    # ---- 10. 最終 smoothing（軽め）----\n    y_final = savgol_filter(y_corrected, window_length=21, polyorder=2)\n\n    # ---- 11. ピクセル → 実電圧変換 ----\n    baseline_row = np.nanmedian(y_final)\n    time_sec = np.arange(W) / (grid_dx * (scale_time / 25.0))\n    voltage_mV = ((baseline_row - y_final) / grid_dy) * (1.0 / (scale_mv / 10.0))\n\n    df_wave = pd.DataFrame({\n        \"time_sec\": time_sec,\n        \"voltage_mV\": voltage_mV,\n        \"voltage_mV_smooth\": savgol_filter(voltage_mV, 21, 3)\n    })\n\n    return df_wave","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-10T05:10:23.889934Z","iopub.execute_input":"2025-12-10T05:10:23.891097Z","iopub.status.idle":"2025-12-10T05:10:23.904387Z","shell.execute_reply.started":"2025-12-10T05:10:23.891065Z","shell.execute_reply":"2025-12-10T05:10:23.903185Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 6. Example Output (Single Lead)\n\nHere we apply v4.1 to a representative ECG image.\nThis plot will serve as a visual reference.\n","metadata":{}},{"cell_type":"code","source":"#\n# Kaggle 上のデータ構造に合わせる\nBASE_PATH = \"/kaggle/input/physionet-ecg-image-digitization\"\nTRAIN_DIR = f\"{BASE_PATH}/train\"\nTEST_DIR = f\"{BASE_PATH}/test\"\n\n# サンプルIDを選択\nsample_id = \"1006427285\"\nimg_name = f\"{sample_id}-0001.png\"\nimg_path = f\"{TRAIN_DIR}/{sample_id}/{img_name}\"\nprint(\"Image path:\", img_path)\n\nprint(\"Image path2:\",\"/kaggle/input/physionet-ecg-image-digitization/train/1006427285/1006427285-0001.png\")\n\n# Read Image\nimg = cv2.imread(img_path)\n##print(\"img:\",np.array(img))\ngray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n\nplt.figure(figsize=(12,6))\nplt.imshow(gray, cmap='gray')\nplt.title(f\"Sample ECG: {img_name}\")\nplt.axis('off')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-10T05:17:04.707814Z","iopub.execute_input":"2025-12-10T05:17:04.708172Z","iopub.status.idle":"2025-12-10T05:17:04.716828Z","shell.execute_reply.started":"2025-12-10T05:17:04.708148Z","shell.execute_reply":"2025-12-10T05:17:04.715693Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =====================================================\n# Multi-File ECG Digitization Pipeline/\\/\\\n# =====================================================\n\nimport os\nimport glob\nimport cv2\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom tqdm import tqdm\n\n# --- 既に定義済みの関数を使う想定 ---\n# from previous cell:\n#   extract_centerline_from_mask(mask, grid_dx, grid_dy, scale_time, scale_mv)\n\n# =====================================================\n# 設定\n# =====================================================\ninput_dir = \"/kaggle/input/physionet-ecg-image-digitization/train/1006867983/\"       # ECG画像を入れたフォルダを指定\noutput_dir = \"./ecg_csv_outputs\"  # CSV出力先\n\nos.makedirs(output_dir, exist_ok=True)\n\n# 画像ファイルの一覧取得（png, jpgどちらも対応）\nimage_files = sorted(glob.glob(os.path.join(input_dir, \"*.png\")) +\n                     glob.glob(os.path.join(input_dir, \"*.jpg\")))\n\nprint(f\"Found {len(image_files)} ECG images.\")\n\n# =====================================================\n# パラメータ設定\n# =====================================================\n#grid_dx = 20   # 1mmあたりのピクセル（横）\n#grid_dy = 20   # 1mmあたりのピクセル（縦）\ngrid_dx = 20 #22   # 画像の横スケールが広めなら少し大きく\ngrid_dy = 10#18   # 縦スケールが狭いなら少し小さく\nscale_time = 25.0  # 25mm/s\nscale_mv = 10.0    # 10mm/mV\n\n\ndef preprocess_ecg_image(img_gray, grid_kernel=30, show=False):\n    \"\"\"\n    ECG画像（グレースケール）からグリッド除去＋波形強調を行う。\n    grid_kernel: グリッド線の周期（px単位、25mm/sなら30前後）\n    \"\"\"\n    # --- ① 背景のグリッド線除去 ---\n    # 横線・縦線をそれぞれ削除\n    horizontal = cv2.morphologyEx(img_gray, cv2.MORPH_OPEN,\n                                  cv2.getStructuringElement(cv2.MORPH_RECT, (grid_kernel, 1)))\n    vertical = cv2.morphologyEx(img_gray, cv2.MORPH_OPEN,\n                                cv2.getStructuringElement(cv2.MORPH_RECT, (1, grid_kernel)))\n    grid = cv2.add(horizontal, vertical)\n    no_grid = cv2.subtract(img_gray, grid)\n    \n    # --- ② 二値化 + モルフォロジ（ノイズ除去＋線補強） ---\n    _, binary = cv2.threshold(no_grid, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)\n    kernel = np.ones((3, 3), np.uint8)\n    morph = cv2.morphologyEx(binary, cv2.MORPH_CLOSE, kernel, iterations=2)\n    morph = cv2.morphologyEx(morph, cv2.MORPH_OPEN, kernel, iterations=1)\n\n    # --- ③ 細線化（1px幅の線に） ---\n    # scikit-image の skeletonizeを使用（入力は0,1のfloat画像）\n    thin = skeletonize((morph > 0).astype(np.uint8)).astype(np.uint8) * 255\n\n    if show:\n        plt.figure(figsize=(14,4))\n        plt.subplot(1,4,1); plt.imshow(img_gray, cmap='gray'); plt.title(\"Gray Input\"); plt.axis('off')\n        plt.subplot(1,4,2); plt.imshow(no_grid, cmap='gray'); plt.title(\"No Grid\"); plt.axis('off')\n        plt.subplot(1,4,3); plt.imshow(morph, cmap='gray'); plt.title(\"Morph\"); plt.axis('off')\n        plt.subplot(1,4,4); plt.imshow(thin, cmap='gray'); plt.title(\"Skeletonized\"); plt.axis('off')\n        plt.show()\n\n    return thin\n\n\n# =====================================================\n# メイン処理ループ\n# =====================================================\nfor path in tqdm(image_files, desc=\"Processing ECG images\"):\n\n    # ---- 画像読込 ----\n    img = cv2.imread(path)\n    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n\n    thin = preprocess_ecg_image(gray, grid_kernel=22)\n    \n    blur = cv2.GaussianBlur(thin, (3, 3), 0)\n    _, binary = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)\n    \n    kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (2, 2))\n    binary = cv2.morphologyEx(binary, cv2.MORPH_OPEN, kernel)\n    binary = cv2.bitwise_not(binary)\n\n    # ---- 波形中心線抽出 ----\n    df_wave = extract_centerline_v41(\n        mask=binary,\n        grid_dx=grid_dx,\n        grid_dy=grid_dy,\n        scale_time=scale_time,\n        scale_mv=scale_mv\n    )\n\n    \n    from scipy.signal import savgol_filter\n    from scipy.ndimage import median_filter\n    \n    # --- 波形補間 + 平滑化 ---\n    # 欠損値補完\n    df_wave[\"voltage_mV\"] = df_wave[\"voltage_mV\"].interpolate(method=\"linear\")\n    \n    # 外れ値除去（3点メディアン）\n    df_wave[\"voltage_mV\"] = median_filter(df_wave[\"voltage_mV\"], size=3)\n    \n    # サヴィツキー・ゴレイで滑らか化\n    df_wave[\"voltage_mV_smooth\"] = savgol_filter(df_wave[\"voltage_mV\"], window_length=11, polyorder=3)\n    \n    # ベースライン補正（平均値をゼロに）\n    df_wave[\"voltage_mV_smooth\"] -= df_wave[\"voltage_mV_smooth\"].mean()\n    \n    # --- プロット ---\n    plt.figure(figsize=(8,3))\n    plt.plot(df_wave[\"time_sec\"], df_wave[\"voltage_mV_smooth\"], color=\"red\", lw=1)\n    #plt.title(f\"Smoothed ECG Trace: {lead_name}\")\n    plt.title(f\"Smoothed ECG Trace: \")\n    plt.xlabel(\"Time [sec]\"); plt.ylabel(\"Voltage [mV]\")\n    plt.grid(True); plt.tight_layout(); plt.show()\n\n\n\n\n    \n    # ---- ファイル名設定 ----\n    base = os.path.basename(path)\n    lead_name = os.path.splitext(base)[0]\n    save_csv = os.path.join(output_dir, f\"{lead_name}.csv\")\n    df_wave.to_csv(save_csv, index=False)\n\n    # ---- 確認プロット ----\n    plt.figure(figsize=(8, 3))\n    plt.plot(df_wave[\"time_sec\"], df_wave[\"voltage_mV\"], color=\"red\")\n    plt.title(f\"Extracted ECG Trace: \")\n    plt.xlabel(\"Time [sec]\")\n    plt.ylabel(\"Voltage [mV]\")\n    plt.grid(True)\n    plt.tight_layout()\n    plt.show()\n\nprint(\"✅ Completed ECG extraction for all files.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-10T05:33:29.149606Z","iopub.execute_input":"2025-12-10T05:33:29.150840Z","iopub.status.idle":"2025-12-10T05:37:04.662426Z","shell.execute_reply.started":"2025-12-10T05:33:29.150803Z","shell.execute_reply":"2025-12-10T05:37:04.661412Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 7. Combined 12-Lead Visualization\n\nTo verify temporal consistency and overall stability,\nwe overlay 12 leads after extracting centerlines.\n\nThis allows checking:\n- relative alignment\n- noise level\n- morphological consistency\n","metadata":{}},{"cell_type":"code","source":"# =====================================================\n# Combine 12-lead ECG CSVs into one aligned record ver2\n# =====================================================\n\nimport os\nimport pandas as pd\nimport numpy as np\nimport glob\nfrom scipy.signal import correlate\n\ndef align_and_combine_leads(csv_folder, max_shift=0.2, sample_rate=500):\n    \"\"\"\n    指定フォルダ内の12リードCSVを1つのDataFrameに整形\n    （時間軸・位相をそろえて統合）\n    \n    Parameters:\n      csv_folder: 各リードのCSVが入ったフォルダ\n      max_shift: 時間シフトの許容範囲 [秒]\n      sample_rate: サンプリング周波数 [Hz] （25mm/s → 約500Hz）\n    Returns:\n      pd.DataFrame: time列＋12リード列\n    \"\"\"\n    csv_files = sorted(glob.glob(os.path.join(csv_folder, \"*.csv\")))\n    if not csv_files:\n        print(f\"⚠️ No CSV files in {csv_folder}\")\n        return None\n\n    leads = {}\n    base_time = None\n\n    for f in csv_files:\n        lead_name = os.path.splitext(os.path.basename(f))[0].split(\"-\")[-1]  # ファイル末尾から誘導名推定\n        df = pd.read_csv(f)\n        if base_time is None:\n            base_time = df[\"time_sec\"].values\n        leads[lead_name] = df[\"voltage_mV\"].values\n\n    # 長さ合わせ\n    min_len = min(len(v) for v in leads.values())\n    for k in leads:\n        leads[k] = leads[k][:min_len]\n    base_time = base_time[:min_len]\n\n    # 時間シフト補正（II誘導を基準）\n    ref = leads.get(\"0001\", list(leads.values())[0])  # 最初のリードを基準に\n    aligned = {}\n    for name, sig in leads.items():\n        corr = correlate(ref - np.mean(ref), sig - np.mean(sig), mode=\"full\")\n        lag = np.argmax(corr) - len(sig) + 1\n        lag_time = lag / sample_rate\n        if abs(lag_time) > max_shift:\n            lag = 0  # 大きすぎるシフトは無視\n        aligned[name] = np.roll(sig, -lag)\n\n    # DataFrameへ統合\n    df_combined = pd.DataFrame({\"time_sec\": base_time})\n    for name, sig in aligned.items():\n        df_combined[name] = sig\n\n    return df_combined\n\n# =====================================================\n# Example: Combine one patient's leads\n# =====================================================\nrecord_folder = \"./ecg_csv_outputs/\"  # 対象フォルダを指定\ndf_12lead = align_and_combine_leads(record_folder)\n\nif df_12lead is not None:\n    display(df_12lead.head())\n    plt.figure(figsize=(10,6))\n    for col in df_12lead.columns[1:]:\n        plt.plot(df_12lead[\"time_sec\"], df_12lead[col] + 2*list(df_12lead.columns[1:]).index(col), label=col)\n    plt.title(f\"12-lead ECG Combined View ({os.path.basename(record_folder)})\")\n    plt.xlabel(\"Time [sec]\")\n    plt.ylabel(\"Relative Voltage [mV]\")\n    plt.legend(ncol=4, fontsize=8)\n    plt.grid(True)\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-10T05:38:07.006664Z","iopub.execute_input":"2025-12-10T05:38:07.007018Z","iopub.status.idle":"2025-12-10T05:38:07.551231Z","shell.execute_reply.started":"2025-12-10T05:38:07.006995Z","shell.execute_reply":"2025-12-10T05:38:07.549747Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 8. SNR Analysis (Preliminary)\n\nCurrent SNR values are low (≈ -33 dB).\nThis is expected because:\n\n- segmentation masks contain jitter noise  \n- no baseline wander correction applied yet  \n- SNR is computed per-lead without template matching  \n\nFuture work will include:\n- baseline correction\n- R-peak aligned template averaging\n- improved SNR metrics\n","metadata":{}},{"cell_type":"code","source":"# =====================================================\n# SNR Local Evaluation Function for PhysioNet ECG\n# =====================================================\n\nimport numpy as np\nfrom scipy.signal import correlate\n\ndef align_signals(pred, truth, sample_rate=500, max_shift=0.2):\n    \"\"\"\n    時間軸シフトを補正して2つの信号を整列。\n    最大±max_shift秒の範囲で相互相関を最大化。\n    \"\"\"\n    max_lag = int(max_shift * sample_rate)\n    corr = correlate(truth - np.mean(truth), pred - np.mean(pred), mode=\"full\")\n    lags = np.arange(-len(pred)+1, len(truth))\n    valid_idx = (lags >= -max_lag) & (lags <= max_lag)\n    lag = lags[valid_idx][np.argmax(corr[valid_idx])]\n    aligned_pred = np.roll(pred, lag)\n    return aligned_pred, lag / sample_rate\n\n\ndef compute_snr(pred, truth):\n    \"\"\"\n    縦オフセットを補正し、信号対雑音比（SNR, dB）を算出。\n    \"\"\"\n    pred = pred - np.mean(pred)\n    truth = truth - np.mean(truth)\n    noise = truth - pred\n    snr = 10 * np.log10(np.sum(truth**2) / np.sum(noise**2))\n    return snr\n\n\ndef evaluate_ecg_record(pred_df, truth_df, sample_rate=500):\n    \"\"\"\n    各リードごとにSNRを計算し、平均を返す。\n    \"\"\"\n    leads = [c for c in truth_df.columns if c != \"time_sec\" and c in pred_df.columns]\n    snr_list = []\n\n    for lead in leads:\n        pred = pred_df[lead].values\n        truth = truth_df[lead].values\n\n        # 長さ調整\n        min_len = min(len(pred), len(truth))\n        pred, truth = pred[:min_len], truth[:min_len]\n\n        # 整列・補正\n        aligned_pred, shift = align_signals(pred, truth, sample_rate=sample_rate)\n        snr = compute_snr(aligned_pred, truth)\n        snr_list.append(snr)\n\n        print(f\"{lead:>5s} | Shift: {shift:+.4f}s | SNR: {snr:.2f} dB\")\n\n    avg_snr = np.mean(snr_list)\n    print(f\"\\n✅ Mean SNR across leads: {avg_snr:.2f} dB\")\n    return avg_snr\n\n# truth_df:  本当のECG時系列（PhysioNet提供のground truth）\n# pred_df:   Keiさんの抽出した波形（df_12leadなど）\n\n# 例: Check 同じ/kaggle/input/physionet-ecg-image-digitization/train/1006867983フォルダ内に正解ファイルがある場合\ntruth_path = \"/kaggle/input/physionet-ecg-image-digitization/train/1006867983/1006867983.csv\"\ntruth_df = pd.read_csv(truth_path)\n\n# Keiさんの推定結果（すでに作成済み）\npred_df = df_12lead.copy()\n\n# 典型的なマッピング例（必要に応じて調整）\nalias_to_num = {\n    \"I\":\"0001\",\"II\":\"0003\",\"III\":\"0004\",\n    \"aVR\":\"0005\",\"aVL\":\"0006\",\"aVF\":\"0009\",\n    \"V1\":\"0010\",\"V2\":\"0011\",\"V3\":\"0012\",\n    \"V4\":\"0013\",\"V5\":\"0014\",\"V6\":\"0015\"  # 末尾の番号は実データに合わせて変更する\n}\n\ndef normalize_lead_columns(df):\n    cols = list(df.columns)\n    # 時間列\n    time_col = \"time_sec\" if \"time_sec\" in cols else (\"time\" if \"time\" in cols else None)\n    # リード列の名前を数字系に寄せる\n    new = {}\n    for c in cols:\n        if c == time_col: \n            new[c] = \"time_sec\"\n        elif c in alias_to_num:\n            new[c] = alias_to_num[c]\n        else:\n            new[c] = c  # すでに 0001 等ならそのまま\n    df = df.rename(columns=new)\n    return df\n\ntruth_df = normalize_lead_columns(truth_df)\npred_df  = normalize_lead_columns(df_12lead.copy())\n\nprint(\"truth_df (norm) cols:\", [c for c in truth_df.columns if c!=\"time_sec\"][:15])\nprint(\"pred_df  (norm) cols:\", [c for c in pred_df.columns  if c!=\"time_sec\"][:15])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-10T05:45:58.047427Z","iopub.execute_input":"2025-12-10T05:45:58.047840Z","iopub.status.idle":"2025-12-10T05:45:58.082529Z","shell.execute_reply.started":"2025-12-10T05:45:58.047813Z","shell.execute_reply":"2025-12-10T05:45:58.080555Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =====================================================\n# SNR Local Evaluation Function for PhysioNet ECG\n# =====================================================\n\nimport numpy as np\nfrom scipy.signal import correlate\n\ndef align_signals(pred, truth, sample_rate=500, max_shift=0.2):\n    \"\"\"\n    時間軸シフトを補正して2つの信号を整列。\n    最大±max_shift秒の範囲で相互相関を最大化。\n    \"\"\"\n    max_lag = int(max_shift * sample_rate)\n    corr = correlate(truth - np.mean(truth), pred - np.mean(pred), mode=\"full\")\n    lags = np.arange(-len(pred)+1, len(truth))\n    valid_idx = (lags >= -max_lag) & (lags <= max_lag)\n    lag = lags[valid_idx][np.argmax(corr[valid_idx])]\n    aligned_pred = np.roll(pred, lag)\n    return aligned_pred, lag / sample_rate\n\n\ndef compute_snr(pred, truth):\n    \"\"\"\n    縦オフセットを補正し、信号対雑音比（SNR, dB）を算出。\n    \"\"\"\n    pred = pred - np.mean(pred)\n    truth = truth - np.mean(truth)\n    noise = truth - pred\n    snr = 10 * np.log10(np.sum(truth**2) / np.sum(noise**2))\n    return snr\n\n\ndef evaluate_ecg_record(pred_df, truth_df, sample_rate=500):\n    \"\"\"\n    各リードごとにSNRを計算し、平均を返す。\n    \"\"\"\n    leads = [c for c in truth_df.columns if c != \"time_sec\" and c in pred_df.columns]\n    snr_list = []\n\n    for lead in leads:\n        pred = pred_df[lead].values\n        truth = truth_df[lead].values\n\n        # 長さ調整\n        min_len = min(len(pred), len(truth))\n        pred, truth = pred[:min_len], truth[:min_len]\n\n        # 整列・補正\n        aligned_pred, shift = align_signals(pred, truth, sample_rate=sample_rate)\n        snr = compute_snr(aligned_pred, truth)\n        snr_list.append(snr)\n\n        print(f\"{lead:>5s} | Shift: {shift:+.4f}s | SNR: {snr:.2f} dB\")\n\n    avg_snr = np.mean(snr_list)\n    print(f\"\\n✅ Mean SNR across leads: {avg_snr:.2f} dB\")\n    return avg_snr","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-10T05:48:12.918766Z","iopub.execute_input":"2025-12-10T05:48:12.919084Z","iopub.status.idle":"2025-12-10T05:48:12.930691Z","shell.execute_reply.started":"2025-12-10T05:48:12.919062Z","shell.execute_reply":"2025-12-10T05:48:12.929671Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nfrom scipy.signal import correlate\n\ndef _get_time_and_leads(df):\n    time_col = \"time_sec\" if \"time_sec\" in df.columns else (\"time\" if \"time\" in df.columns else None)\n    if time_col is None:\n        raise ValueError(\"時間列(time/time_sec)が見つかりません。\")\n    leads = [c for c in df.columns if c != time_col]\n    return time_col, leads\n\ndef _align_by_xcorr(pred, truth, sample_rate, max_shift=0.2):\n    max_lag = int(max_shift * sample_rate)\n    corr = correlate(truth - np.mean(truth), pred - np.mean(pred), mode=\"full\")\n    lags = np.arange(-len(pred)+1, len(truth))\n    sel = (lags >= -max_lag) & (lags <= max_lag)\n    lag = lags[sel][np.argmax(corr[sel])] if np.any(sel) else 0\n    return np.roll(pred, lag), lag / sample_rate\n\ndef evaluate_ecg_record_v2(pred_df, truth_df):\n    # 時間列取得\n    p_time, p_leads = _get_time_and_leads(pred_df)\n    t_time, t_leads = _get_time_and_leads(truth_df)\n    # 共通リード\n    common = [l for l in p_leads if l in t_leads]\n    if not common:\n        raise ValueError(f\"共通するリードがありません。\\n pred:{p_leads}\\n truth:{t_leads}\")\n    # サンプルレート推定（時間差の中央値から）\n    dt_pred = np.median(np.diff(pred_df[p_time].values))\n    sample_rate = int(round(1.0 / dt_pred)) if dt_pred > 0 else 500\n\n    snrs = []\n    for lead in common:\n        p = pred_df[lead].values\n        t = truth_df[lead].values\n        n = min(len(p), len(t))\n        p, t = p[:n], t[:n]\n        p_aligned, shift = _align_by_xcorr(p, t, sample_rate=sample_rate, max_shift=0.2)\n        # 縦オフセット除去\n        p_aligned = p_aligned - np.mean(p_aligned)\n        t = t - np.mean(t)\n        noise = t - p_aligned\n        if np.sum(noise**2) == 0:\n            snr = np.inf\n        else:\n            snr = 10*np.log10(np.sum(t**2)/np.sum(noise**2))\n        print(f\"{lead:>5s} | shift={shift:+.4f}s | SNR={snr:.2f} dB\")\n        snrs.append(snr)\n    mean_snr = float(np.mean(snrs))\n    print(f\"\\n✅ Mean SNR across leads: {mean_snr:.2f} dB (n={len(snrs)})\")\n    return mean_snr\n\n# 実行\nmean_snr = evaluate_ecg_record_v2(pred_df, truth_df)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-10T05:48:17.521638Z","iopub.execute_input":"2025-12-10T05:48:17.522094Z","iopub.status.idle":"2025-12-10T05:48:17.568408Z","shell.execute_reply.started":"2025-12-10T05:48:17.522063Z","shell.execute_reply":"2025-12-10T05:48:17.566619Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from scipy.signal import correlate\n\ndef _infer_fs_from_time(df):\n    tcol = \"time_sec\" if \"time_sec\" in df.columns else (\"time\" if \"time\" in df.columns else None)\n    if tcol is None:\n        return 500  # フォールバック\n    dt = np.median(np.diff(df[tcol].values))\n    return int(round(1.0/dt)) if dt and np.isfinite(dt) and dt>0 else 500\n\ndef _align_xcorr_valid(pred, truth, fs, max_shift_s=0.2):\n    \"\"\"NaNを除外しつつ相互相関でシフト最適化（±max_shift_s）\"\"\"\n    max_lag = int(max_shift_s * fs)\n    # NaN除外\n    mask = np.isfinite(pred) & np.isfinite(truth)\n    if mask.sum() < 10:  # データ少なすぎ\n        return pred, 0.0\n    p = pred[mask] - np.nanmean(pred[mask])\n    q = truth[mask] - np.nanmean(truth[mask])\n    # 相関\n    corr = correlate(q, p, mode=\"full\")\n    lags = np.arange(-len(p)+1, len(q))\n    ok = (lags >= -max_lag) & (lags <= max_lag)\n    lag = lags[ok][np.argmax(corr[ok])] if ok.any() else 0\n    # シフト適用（全長に対して）\n    pred_shift = np.roll(pred, lag)\n    return pred_shift, lag / fs\n\ndef evaluate_ecg_record_robust(pred_df, truth_df, fs=None, max_shift_s=0.2, verbose=True):\n    \"\"\"\n    各リードで truth の「有限値が連続している先頭区間」だけを使い、\n    pred を同じ長さに切り出してから整列＆SNR計算。\n    \"\"\"\n    if fs is None:\n        fs = _infer_fs_from_time(pred_df)\n\n    p_time = \"time_sec\" if \"time_sec\" in pred_df.columns else \"time\"\n    t_time = \"time_sec\" if \"time_sec\" in truth_df.columns else \"time\"\n\n    pred_leads = [c for c in pred_df.columns if c != p_time]\n    truth_leads = [c for c in truth_df.columns if c != t_time]\n    common = [l for l in pred_leads if l in truth_leads]\n    if not common:\n        raise ValueError(f\"共通リードがありません: pred={pred_leads}, truth={truth_leads}\")\n\n    snr_list = []\n    for lead in common:\n        p = pred_df[lead].values.astype(float)\n        t = truth_df[lead].values.astype(float)\n\n        # ---- truth の「先頭から連続する有限区間」のみ使う ----\n        finite_truth = np.isfinite(t)\n        if not finite_truth.any():\n            if verbose: print(f\"{lead:>5s} | no finite truth -> skip\")\n            continue\n        # 先頭からの連続区間長\n        valid_len = np.argmax(~finite_truth) if (~finite_truth).any() else len(t)\n        t_seg = t[:valid_len]\n        p_seg = p[:valid_len] if len(p) >= valid_len else p  # predが短ければそのまま\n\n        # 長さ合わせ\n        n = min(len(p_seg), len(t_seg))\n        if n < 10:\n            if verbose: print(f\"{lead:>5s} | too short overlap -> skip\")\n            continue\n        p_seg = p_seg[:n]\n        t_seg = t_seg[:n]\n\n        # 整列（NaNを考慮した相関）\n        p_aligned, shift = _align_xcorr_valid(p_seg, t_seg, fs, max_shift_s=max_shift_s)\n\n        # オフセット除去 & NaN除去\n        mask = np.isfinite(p_aligned) & np.isfinite(t_seg)\n        if mask.sum() < 10:\n            if verbose: print(f\"{lead:>5s} | not enough finite -> skip\")\n            continue\n        p0 = p_aligned[mask] - np.mean(p_aligned[mask])\n        t0 = t_seg[mask] - np.mean(t_seg[mask])\n\n        noise = t0 - p0\n        denom = np.sum(noise**2)\n        if denom == 0:\n            snr = np.inf\n        else:\n            snr = 10*np.log10(np.sum(t0**2) / denom)\n\n        if verbose:\n            print(f\"{lead:>5s} | len={len(t0):4d} | shift={shift:+.4f}s | SNR={snr:.2f} dB\")\n        snr_list.append(snr)\n\n    if len(snr_list)==0:\n        if verbose: print(\"⚠️ 有効なリードがありません（NaNのみ）。\")\n        return np.nan\n\n    mean_snr = float(np.nanmean(snr_list))\n    if verbose:\n        print(f\"\\n✅ Mean SNR across leads: {mean_snr:.2f} dB (n={len(snr_list)})\")\n    return mean_snr","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-10T05:49:09.151769Z","iopub.execute_input":"2025-12-10T05:49:09.152102Z","iopub.status.idle":"2025-12-10T05:49:09.175070Z","shell.execute_reply.started":"2025-12-10T05:49:09.152081Z","shell.execute_reply":"2025-12-10T05:49:09.173472Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =====================================================\n# Example Usage\n# =====================================================\n\n# truth_df:  本当のECG時系列（PhysioNet提供のground truth）\n# pred_df:   Keiさんの抽出した波形（df_12leadなど）\n\n# 例: Check 同じ/kaggle/input/physionet-ecg-image-digitization/train/1006867983フォルダ内に正解ファイルがある場合\ntruth_path = \"/kaggle/input/physionet-ecg-image-digitization/train/1006867983/1006867983.csv\"\ntruth_df = pd.read_csv(truth_path)\n\n# Keiさんの推定結果（すでに作成済み）\npred_df = df_12lead.copy()\n\n# 典型的なマッピング例（必要に応じて調整）\nalias_to_num = {\n    \"I\":\"0001\",\"II\":\"0003\",\"III\":\"0004\",\n    \"aVR\":\"0005\",\"aVL\":\"0006\",\"aVF\":\"0009\",\n    \"V1\":\"0010\",\"V2\":\"0011\",\"V3\":\"0012\",\n    \"V4\":\"0013\",\"V5\":\"0014\",\"V6\":\"0015\"  # 末尾の番号は実データに合わせて変更する\n}\n\ndef normalize_lead_columns(df):\n    cols = list(df.columns)\n    # 時間列\n    time_col = \"time_sec\" if \"time_sec\" in cols else (\"time\" if \"time\" in cols else None)\n    # リード列の名前を数字系に寄せる\n    new = {}\n    for c in cols:\n        if c == time_col: \n            new[c] = \"time_sec\"\n        elif c in alias_to_num:\n            new[c] = alias_to_num[c]\n        else:\n            new[c] = c  # すでに 0001 等ならそのまま\n    df = df.rename(columns=new)\n    return df\n\ntruth_df = normalize_lead_columns(truth_df)\npred_df  = normalize_lead_columns(df_12lead.copy())\n\nprint(\"truth_df (norm) cols:\", [c for c in truth_df.columns if c!=\"time_sec\"][:15])\nprint(\"pred_df  (norm) cols:\", [c for c in pred_df.columns  if c!=\"time_sec\"][:15])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-10T05:49:35.666123Z","iopub.execute_input":"2025-12-10T05:49:35.666471Z","iopub.status.idle":"2025-12-10T05:49:35.685035Z","shell.execute_reply.started":"2025-12-10T05:49:35.666434Z","shell.execute_reply":"2025-12-10T05:49:35.683651Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"mean_snr = evaluate_ecg_record_robust(pred_df, truth_df)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-10T05:50:35.150384Z","iopub.execute_input":"2025-12-10T05:50:35.150736Z","iopub.status.idle":"2025-12-10T05:50:35.165167Z","shell.execute_reply.started":"2025-12-10T05:50:35.150707Z","shell.execute_reply":"2025-12-10T05:50:35.163864Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\n\ntruth_dict = {}\nmax_len = max(len(sig) for sig in truth_signals.values())\ntime_sec = np.arange(max_len) / sample_rate\ntruth_dict[\"time_sec\"] = time_sec\n\nfor lead, sig in truth_signals.items():\n    pad = np.full(max_len, np.nan)\n    pad[:len(sig)] = sig\n    truth_dict[lead] = pad\n\ntruth_df_fixed = pd.DataFrame(truth_dict)\nprint(truth_df_fixed.shape)\ndisplay(truth_df_fixed.head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-10T05:51:27.010815Z","iopub.execute_input":"2025-12-10T05:51:27.011110Z","iopub.status.idle":"2025-12-10T05:51:27.025343Z","shell.execute_reply.started":"2025-12-10T05:51:27.011090Z","shell.execute_reply":"2025-12-10T05:51:27.023766Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"## 2️⃣ ECG画像を読み込み・可視化","metadata":{}},{"cell_type":"code","source":"\n# 画像の読み込み\nimg = cv2.imread(img_path)\n##print(\"img:\",np.array(img))\ngray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n\nplt.figure(figsize=(12,6))\nplt.imshow(gray, cmap='gray')\nplt.title(f\"Sample ECG: {img_name}\")\nplt.axis('off')\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:17:27.306069Z","iopub.execute_input":"2025-11-29T06:17:27.306515Z","iopub.status.idle":"2025-11-29T06:17:28.158024Z","shell.execute_reply.started":"2025-11-29T06:17:27.30649Z","shell.execute_reply":"2025-11-29T06:17:28.156892Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 3️⃣ グリッド線の検出（時間軸と振幅軸の目安）","metadata":{}},{"cell_type":"code","source":"\nblur = cv2.GaussianBlur(gray, (3,3), 0)\nedges = cv2.Canny(blur, 50, 150, apertureSize=3)\n\nlines = cv2.HoughLines(edges, 1, np.pi / 180, 400)\ngrid = img.copy()\n\nif lines is not None:\n    for rho, theta in lines[:,0]:\n        a = np.cos(theta); b = np.sin(theta)\n        x0 = a * rho; y0 = b * rho\n        x1 = int(x0 + 1000 * (-b)); y1 = int(y0 + 1000 * (a))\n        x2 = int(x0 - 1000 * (-b)); y2 = int(y0 - 1000 * (a))\n        cv2.line(grid, (x1, y1), (x2, y2), (0,255,0), 1)\n\nplt.figure(figsize=(12,6))\nplt.imshow(cv2.cvtColor(grid, cv2.COLOR_BGR2RGB))\nplt.title(\"Detected Grid Lines (ECG background)\")\nplt.axis('off')\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:17:33.901484Z","iopub.execute_input":"2025-11-29T06:17:33.901864Z","iopub.status.idle":"2025-11-29T06:17:34.863155Z","shell.execute_reply.started":"2025-11-29T06:17:33.901839Z","shell.execute_reply":"2025-11-29T06:17:34.861927Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 4️⃣ 波形抽出（二値化＋輪郭検出）","metadata":{}},{"cell_type":"code","source":"# =====================================================\n# ① グリッド除去 + ② モルフォロジ + ③ 細線化 自動処理\n# =====================================================\nimport cv2\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom skimage.morphology import skeletonize\n\ndef preprocess_ecg_image(img_gray, grid_kernel=30, show=False):\n    \"\"\"\n    ECG画像（グレースケール）からグリッド除去＋波形強調を行う。\n    grid_kernel: グリッド線の周期（px単位、25mm/sなら30前後）\n    \"\"\"\n    # --- ① 背景のグリッド線除去 ---\n    # 横線・縦線をそれぞれ削除\n    horizontal = cv2.morphologyEx(img_gray, cv2.MORPH_OPEN,\n                                  cv2.getStructuringElement(cv2.MORPH_RECT, (grid_kernel, 1)))\n    vertical = cv2.morphologyEx(img_gray, cv2.MORPH_OPEN,\n                                cv2.getStructuringElement(cv2.MORPH_RECT, (1, grid_kernel)))\n    grid = cv2.add(horizontal, vertical)\n    no_grid = cv2.subtract(img_gray, grid)\n    \n    # --- ② 二値化 + モルフォロジ（ノイズ除去＋線補強） ---\n    _, binary = cv2.threshold(no_grid, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)\n    kernel = np.ones((3, 3), np.uint8)\n    morph = cv2.morphologyEx(binary, cv2.MORPH_CLOSE, kernel, iterations=2)\n    morph = cv2.morphologyEx(morph, cv2.MORPH_OPEN, kernel, iterations=1)\n\n    # --- ③ 細線化（1px幅の線に） ---\n    # scikit-image の skeletonizeを使用（入力は0,1のfloat画像）\n    thin = skeletonize((morph > 0).astype(np.uint8)).astype(np.uint8) * 255\n\n    if show:\n        plt.figure(figsize=(14,4))\n        plt.subplot(1,4,1); plt.imshow(img_gray, cmap='gray'); plt.title(\"Gray Input\"); plt.axis('off')\n        plt.subplot(1,4,2); plt.imshow(no_grid, cmap='gray'); plt.title(\"No Grid\"); plt.axis('off')\n        plt.subplot(1,4,3); plt.imshow(morph, cmap='gray'); plt.title(\"Morph\"); plt.axis('off')\n        plt.subplot(1,4,4); plt.imshow(thin, cmap='gray'); plt.title(\"Skeletonized\"); plt.axis('off')\n        plt.show()\n\n    return thin\n\n# --- テスト例 ---\n# img_gray = cv2.imread('/kaggle/input/physionet-ecg-image-digitization/train/1006867983/1006867983-0001.png', cv2.IMREAD_GRAYSCALE)\n# processed = preprocess_ecg_image(img_gray, grid_kernel=30, show=True)\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:17:40.431005Z","iopub.execute_input":"2025-11-29T06:17:40.431369Z","iopub.status.idle":"2025-11-29T06:17:41.200171Z","shell.execute_reply.started":"2025-11-29T06:17:40.431343Z","shell.execute_reply":"2025-11-29T06:17:41.199057Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# ---- 画像読込 ----\n#img = cv2.imread(path)\n#gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n\nthin = preprocess_ecg_image(gray, grid_kernel=22)\n\ngray_u8 = thin.copy()\nif gray_u8.dtype != np.uint8:\n    gray_u8 = (np.clip(gray_u8, 0, 255)).astype(np.uint8)\n\n# 1) 白い“細い線”を強調（Top-hat）\nkernel = cv2.getStructuringElement(cv2.MORPH_RECT, (7,7))\ntophat = cv2.morphologyEx(gray_u8, cv2.MORPH_TOPHAT, kernel)\n\n# --- Before thresholding ---\nb, g, r = cv2.split(img)\ngray_wave = cv2.subtract(255, b)  # 反転して黒波形を強調\nplt.imshow(gray_wave, cmap='gray')\nplt.title(\"Enhanced waveform channel\")\nplt.axis('off')\nplt.show()\n\n# Otsu’s binarization\nblur = cv2.GaussianBlur(gray_wave, (5,5), 0)\n_, binary = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)\n\nplt.imshow(binary, cmap='gray')\nplt.title(\"Otsu Threshold Result\")\nplt.axis('off')\nplt.show()\n\nkernel = np.ones((3,3), np.uint8)\nbinary_closed = cv2.morphologyEx(binary, cv2.MORPH_CLOSE, kernel, iterations=2)\n\nplt.imshow(binary_closed, cmap='gray')\nplt.title(\"After Morphological Closing\")\nplt.axis('off')\nplt.show()\n\ncontours, _ = cv2.findContours(binary_closed, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\nwave = np.zeros_like(binary_closed)\ncv2.drawContours(wave, contours, -1, 255, 1)\n\nplt.imshow(wave, cmap='gray')\nplt.title(\"Improved Waveform Contours\")\nplt.axis('off')\nplt.show()\n\n\n\n# 2) “非反転”の適応二値化（白=波形）\nbin_wave = cv2.adaptiveThreshold(\n    tophat, 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 55, -5\n)\n\n# 3) グリッド抑圧（水平/垂直の長いカーネルで開演算して引き算）\nhker = cv2.getStructuringElement(cv2.MORPH_RECT, (45,1))\nvker = cv2.getStructuringElement(cv2.MORPH_RECT, (1,45))\ngrid_h = cv2.morphologyEx(bin_wave, cv2.MORPH_OPEN, hker)\ngrid_v = cv2.morphologyEx(bin_wave, cv2.MORPH_OPEN, vker)\ngrid = cv2.bitwise_or(grid_h, grid_v)\n\nbin_wo_grid = cv2.bitwise_and(bin_wave, cv2.bitwise_not(grid))\n\n# 4) 輪郭検出（外枠だけに偏らないよう RETR_LIST）\ncontours, _ = cv2.findContours(bin_wo_grid, cv2.RETR_LIST, cv2.CHAIN_APPROX_NONE)\n\n# 5) 描画（ゼロ初期化の8bitキャンバスに白で描画）\nwave = np.zeros_like(gray_u8)             # ← 同形状・同dtype(uint8)のゼロ配列\ncv2.drawContours(wave, contours, -1, 255, 1)\n\nplt.figure(figsize=(12,6))\nplt.imshow(binary, cmap='gray')\nplt.title(\"Binary ECG Trace\")\nplt.axis('off')\nplt.show()\n\nplt.figure(figsize=(12,6))\nplt.imshow(wave, cmap='gray')\nplt.title(\"Extracted ECG Line Contours\")\nplt.axis('off')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:17:46.254524Z","iopub.execute_input":"2025-11-29T06:17:46.25563Z","iopub.status.idle":"2025-11-29T06:17:55.358612Z","shell.execute_reply.started":"2025-11-29T06:17:46.255595Z","shell.execute_reply":"2025-11-29T06:17:55.35776Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- Step 1. Morphological closing (with larger kernel) ---\nkernel = np.ones((5,5), np.uint8)\nbinary_closed = cv2.morphologyEx(binary, cv2.MORPH_CLOSE, kernel, iterations=2)\n\n# --- Step 2. Ensure binary is uint8 & white-on-black ---\nbinary_inv = cv2.bitwise_not(binary_closed)\nbinary_inv = np.uint8(binary_inv)\n\nplt.imshow(binary_inv, cmap='gray')\nplt.title(\"Binary (Inverted for Contour Detection)\")\nplt.axis('off')\nplt.show()\n\n# --- Step 3. Detect contours ---\ncontours, _ = cv2.findContours(binary_inv, cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)\n\nwave = np.zeros_like(binary_inv)\ncv2.drawContours(wave, contours, -1, 255, 1)\n\nplt.imshow(wave, cmap='gray')\nplt.title(\"Improved Waveform Contours (Fixed)\")\nplt.axis('off')\nplt.show()\n\nprint(\"Detected contour count:\", len(contours))\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:17:57.419671Z","iopub.execute_input":"2025-11-29T06:17:57.420504Z","iopub.status.idle":"2025-11-29T06:17:58.275938Z","shell.execute_reply.started":"2025-11-29T06:17:57.420474Z","shell.execute_reply":"2025-11-29T06:17:58.274981Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"wave = np.zeros_like(binary_inv)\ncv2.drawContours(wave, contours, -1, 255, 3)   # ← thicknessを3に\nplt.imshow(wave, cmap='gray', vmin=0, vmax=255)\nplt.title(\"Contours (Thicker, visible)\")\nplt.axis('off')\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:18:07.559659Z","iopub.execute_input":"2025-11-29T06:18:07.560718Z","iopub.status.idle":"2025-11-29T06:18:08.133073Z","shell.execute_reply.started":"2025-11-29T06:18:07.560685Z","shell.execute_reply":"2025-11-29T06:18:08.131928Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"b, g, r = cv2.split(img)\n# 赤グリッド除去 → 黒波形強調\ngray_wave = cv2.subtract(255, r)  \nblur = cv2.GaussianBlur(gray_wave, (5,5), 0)\n_, binary = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)\n\nplt.imshow(binary, cmap='gray')\nplt.title(\"Red Grid Removed — Black Wave Enhanced\")\nplt.axis('off')\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:18:08.829164Z","iopub.execute_input":"2025-11-29T06:18:08.829509Z","iopub.status.idle":"2025-11-29T06:18:09.265334Z","shell.execute_reply.started":"2025-11-29T06:18:08.829484Z","shell.execute_reply":"2025-11-29T06:18:09.264234Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"filtered_contours = [cnt for cnt in contours if cv2.contourArea(cnt) > 20]\n\nwave = np.zeros_like(binary_inv)\ncv2.drawContours(wave, filtered_contours, -1, 255, 2)\nplt.imshow(wave, cmap='gray')\nplt.title(\"Filtered Contours (Waveform only)\")\nplt.axis('off')\nplt.show()\n\nprint(\"After filtering:\", len(filtered_contours), \"contours kept\")\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:18:15.289469Z","iopub.execute_input":"2025-11-29T06:18:15.289845Z","iopub.status.idle":"2025-11-29T06:18:15.80885Z","shell.execute_reply.started":"2025-11-29T06:18:15.289804Z","shell.execute_reply":"2025-11-29T06:18:15.807895Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- Grid color suppression (remove reddish pixels) ---\nimg_hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)\n# Hue 0〜10 or 160〜180（赤領域）を除去\nmask_red = cv2.inRange(img_hsv, (0,50,50), (10,255,255)) | cv2.inRange(img_hsv, (160,50,50), (180,255,255))\nimg_no_grid = cv2.inpaint(img, mask_red, 3, cv2.INPAINT_TELEA)\n\nplt.figure(figsize=(12,6))\nplt.imshow(cv2.cvtColor(img_no_grid, cv2.COLOR_BGR2RGB))\nplt.title(\"Grid Removed using HSV Mask\")\nplt.axis('off')\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:18:45.688469Z","iopub.execute_input":"2025-11-29T06:18:45.688798Z","iopub.status.idle":"2025-11-29T06:18:47.563383Z","shell.execute_reply.started":"2025-11-29T06:18:45.688776Z","shell.execute_reply":"2025-11-29T06:18:47.562381Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from skimage.morphology import skeletonize\nwave_bin = (wave > 0).astype(np.uint8)\nskeleton = skeletonize(wave_bin).astype(np.uint8) * 255\n\nplt.imshow(skeleton, cmap='gray')\nplt.title(\"Skeletonized ECG waveform (1px width)\")\nplt.axis('off')\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:18:56.216449Z","iopub.execute_input":"2025-11-29T06:18:56.216783Z","iopub.status.idle":"2025-11-29T06:18:56.707592Z","shell.execute_reply.started":"2025-11-29T06:18:56.216758Z","shell.execute_reply":"2025-11-29T06:18:56.706566Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#　改善コード","metadata":{}},{"cell_type":"code","source":"# --- 1️⃣ 赤系＋ピンク系グリッド除去（範囲拡大） ---\nimg_hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)\n\n# 赤～ピンクまで広めに指定\nmask_red = (\n    cv2.inRange(img_hsv, (0, 30, 40), (15, 255, 255)) |  # 赤系\n    cv2.inRange(img_hsv, (150, 30, 40), (180, 255, 255))  # ピンク系\n)\n\nimg_no_grid = cv2.inpaint(img, mask_red, 3, cv2.INPAINT_TELEA)\n\nplt.figure(figsize=(12,6))\nplt.imshow(cv2.cvtColor(img_no_grid, cv2.COLOR_BGR2RGB))\nplt.title(\"Grid Fully Removed (Red + Pink HSV Range)\")\nplt.axis('off')\nplt.show()\n\n\n# --- 2️⃣ 黒波形を強調し、スケルトン化 ---\nb, g, r = cv2.split(img_no_grid)\ngray_wave = cv2.subtract(255, r)\nblur = cv2.GaussianBlur(gray_wave, (5,5), 0)\n_, binary = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)\n\n# 膨張＋収縮で線を整える\nkernel = np.ones((3,3), np.uint8)\nbinary = cv2.morphologyEx(binary, cv2.MORPH_CLOSE, kernel, iterations=2)\n\n# スケルトン化\nfrom skimage.morphology import skeletonize\nskeleton = skeletonize((binary > 0).astype(np.uint8)).astype(np.uint8) * 255\n\nplt.figure(figsize=(12,6))\nplt.imshow(skeleton, cmap='gray')\nplt.title(\"Final Skeletonized ECG waveform (Grid-Free, 1px)\")\nplt.axis('off')\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:19:03.594723Z","iopub.execute_input":"2025-11-29T06:19:03.595676Z","iopub.status.idle":"2025-11-29T06:19:08.044197Z","shell.execute_reply.started":"2025-11-29T06:19:03.595643Z","shell.execute_reply":"2025-11-29T06:19:08.043219Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\n# スケルトン画像から各列の波形座標を抽出\ny_positions = []\nfor x in range(skeleton.shape[1]):\n    ys = np.where(skeleton[:, x] > 0)[0]\n    if len(ys) > 0:\n        y_mean = np.median(ys)\n        y_positions.append(y_mean)\n    else:\n        y_positions.append(np.nan)\n\n# 補間して欠損を埋める\ny_series = pd.Series(y_positions).interpolate().fillna(method='bfill').fillna(method='ffill')\n\n# CSVとして保存\npd.DataFrame({\"x\": np.arange(len(y_series)), \"y\": y_series}).to_csv(\"ecg_waveform_extracted.csv\", index=False)\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:19:14.059943Z","iopub.execute_input":"2025-11-29T06:19:14.060255Z","iopub.status.idle":"2025-11-29T06:19:14.160915Z","shell.execute_reply.started":"2025-11-29T06:19:14.060233Z","shell.execute_reply":"2025-11-29T06:19:14.15957Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#修正\n# 例：Lead II の位置（画像の下1/4あたり）を切り出す\nh, w = skeleton.shape\nlead_y_start = int(h * 0.70)\nlead_y_end = int(h * 0.80)\nlead_ii = skeleton[lead_y_start:lead_y_end, :]\nplt.imshow(lead_ii, cmap='gray')\nplt.title(\"Extracted Lead II Region\")\nplt.axis('off')\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:19:22.11241Z","iopub.execute_input":"2025-11-29T06:19:22.112873Z","iopub.status.idle":"2025-11-29T06:19:22.21935Z","shell.execute_reply.started":"2025-11-29T06:19:22.112806Z","shell.execute_reply":"2025-11-29T06:19:22.217679Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"y_positions = []\nfor x in range(lead_ii.shape[1]):\n    ys = np.where(lead_ii[:, x] > 0)[0]\n    if len(ys) > 0:\n        y_positions.append(np.median(ys))\n    else:\n        y_positions.append(np.nan)\n\n# 欠損補間 & 平滑化\nimport pandas as pd\nimport numpy as np\nfrom scipy.signal import savgol_filter\n\ny_series = pd.Series(y_positions).interpolate().fillna(method='bfill').fillna(method='ffill')\ny_smooth = savgol_filter(y_series, window_length=31, polyorder=3)  # 平滑化\n\nplt.figure(figsize=(12,4))\nplt.plot(y_smooth)\nplt.gca().invert_yaxis()\nplt.title(\"Smoothed Single Lead ECG (Lead II)\")\nplt.xlabel(\"Time (pixels)\")\nplt.ylabel(\"Amplitude (pixels)\")\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:19:23.356766Z","iopub.execute_input":"2025-11-29T06:19:23.357135Z","iopub.status.idle":"2025-11-29T06:19:24.591021Z","shell.execute_reply.started":"2025-11-29T06:19:23.357111Z","shell.execute_reply":"2025-11-29T06:19:24.590025Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom scipy.signal import savgol_filter\nimport matplotlib.pyplot as plt\n\n# === 1. スケルトン波形から1リード領域切り出し ===\nh, w = skeleton.shape\nlead_y_start = int(h * 0.70)   # 前回より少し狭めに設定\nlead_y_end   = int(h * 0.78)\nlead_ii = skeleton[lead_y_start:lead_y_end, :]\n\n# === 2. 各列ごとの白画素の中央値（波形） ===\ny_positions = []\nfor x in range(lead_ii.shape[1]):\n    ys = np.where(lead_ii[:, x] > 0)[0]\n    if len(ys) > 0:\n        y_positions.append(np.median(ys))\n    else:\n        y_positions.append(np.nan)\n\n# 欠損補間\ny_series = pd.Series(y_positions).interpolate().fillna(method='bfill').fillna(method='ffill')\n\n# === 3. ベースライン補正（移動平均差し引き） ===\nbaseline = y_series.rolling(window=200, center=True, min_periods=1).mean()\ny_corrected = y_series - baseline\n\n# === 4. 平滑化（適度なwindowとpoly） ===\ny_smooth = savgol_filter(y_corrected, window_length=51, polyorder=3)\n\n# === 5. 表示（上下反転で心電図らしく） ===\nplt.figure(figsize=(12,4))\nplt.plot(-y_smooth, color='steelblue')\nplt.title(\"Lead II — Baseline Corrected & Smoothed ECG Waveform\")\nplt.xlabel(\"Time (pixels)\")\nplt.ylabel(\"Amplitude (a.u.)\")\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:19:29.520595Z","iopub.execute_input":"2025-11-29T06:19:29.521223Z","iopub.status.idle":"2025-11-29T06:19:29.78208Z","shell.execute_reply.started":"2025-11-29T06:19:29.521194Z","shell.execute_reply":"2025-11-29T06:19:29.780399Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 5️⃣ 波形をピクセル座標からデジタル化","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(12,4))\nplt.plot(y_series.values)\nplt.gca().invert_yaxis()\nplt.title(\"Digitized ECG waveform\")\nplt.xlabel(\"Time (pixels)\")\nplt.ylabel(\"Amplitude (pixels)\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:19:37.276511Z","iopub.execute_input":"2025-11-29T06:19:37.276923Z","iopub.status.idle":"2025-11-29T06:19:37.486923Z","shell.execute_reply.started":"2025-11-29T06:19:37.2769Z","shell.execute_reply":"2025-11-29T06:19:37.485936Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\nys = []\nfor x in range(wave.shape[1]):\n    ys_in_col = np.where(wave[:,x] > 0)[0]\n    if len(ys_in_col) > 0:\n        ys.append(np.mean(ys_in_col))\n    else:\n        ys.append(np.nan)\n\nys = np.array(ys)\nys = pd.Series(ys).interpolate().fillna(method='bfill').fillna(method='ffill')\n\nplt.figure(figsize=(12,4))\nplt.plot(ys.values)\nplt.gca().invert_yaxis()\nplt.title(\"Digitized ECG waveform\")\nplt.xlabel(\"Time (pixels)\")\nplt.ylabel(\"Amplitude (pixels)\")\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:19:42.592421Z","iopub.execute_input":"2025-11-29T06:19:42.593436Z","iopub.status.idle":"2025-11-29T06:19:42.826944Z","shell.execute_reply.started":"2025-11-29T06:19:42.593406Z","shell.execute_reply":"2025-11-29T06:19:42.82592Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 6️⃣ 波形を正規化してCSV保存","metadata":{}},{"cell_type":"code","source":"\ntime = np.arange(len(ys)) / 500  # 仮の500Hzサンプリング\namp = (ys - ys.mean()) / ys.std()\n\ndf_wave = pd.DataFrame({\"time\": time, \"amplitude\": amp})\nout_path = \"/kaggle/working/sample_waveform.csv\"\ndf_wave.to_csv(out_path, index=False)\nprint(f\"Saved → {out_path}\")\ndf_wave.head()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:19:49.935255Z","iopub.execute_input":"2025-11-29T06:19:49.935558Z","iopub.status.idle":"2025-11-29T06:19:49.973656Z","shell.execute_reply.started":"2025-11-29T06:19:49.935534Z","shell.execute_reply":"2025-11-29T06:19:49.972859Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom scipy.signal import find_peaks, savgol_filter\n\n# === 1️⃣ 自動リード帯検出 ===\n# 各行方向の画素数をカウントして最大値の帯域をLead IIと推定\nproj = np.sum(skeleton > 0, axis=1)\npeaks, _ = find_peaks(proj, distance=100, prominence=100)\nlead_index = peaks[1] if len(peaks) > 1 else peaks[0]  # 2番目のピークをLead IIとみなす\nband = 30  # 帯域幅 (±px)\nlead_ii = skeleton[lead_index-band:lead_index+band, :]\n\nplt.figure(figsize=(10,3))\nplt.imshow(lead_ii, cmap='gray')\nplt.title(\"Auto-detected Lead II region\")\nplt.axis('off')\nplt.show()\n\n\n# === 2️⃣ 列ごとに前列との連続性で最も近い点を採用 ===\ny_positions = []\nprev_y = None\nfor x in range(lead_ii.shape[1]):\n    ys = np.where(lead_ii[:, x] > 0)[0]\n    if len(ys) == 0:\n        y_positions.append(np.nan)\n        continue\n    if prev_y is None:\n        y = np.median(ys)\n    else:\n        y = ys[np.argmin(np.abs(ys - prev_y))]  # 最も近い点を選択\n    y_positions.append(y)\n    prev_y = y\n\ny_series = pd.Series(y_positions).interpolate().fillna(method='bfill').fillna(method='ffill')\n\n# === 3️⃣ ベースライン補正 ===\nbaseline = y_series.rolling(window=200, center=True, min_periods=1).median()\ny_corrected = y_series - baseline\n\n# === 4️⃣ 外れ値除去（Hampelフィルタ） ===\ndef hampel(x, k=25, nsigma=3):\n    s = pd.Series(x)\n    med = s.rolling(k, center=True).median()\n    mad = (s - med).abs().rolling(k, center=True).median() * 1.4826\n    outlier = (np.abs(s - med) > nsigma * mad)\n    s[outlier] = med[outlier]\n    return s.to_numpy()\n\ny_denoised = hampel(y_corrected)\n\n# === 5️⃣ 自動平滑化（R–R間隔から窓決定） ===\npeaks, _ = find_peaks(-y_denoised, distance=150)\nif len(peaks) > 1:\n    rr = np.median(np.diff(peaks))\n    win = int(max(31, rr // 3) // 2 * 2 + 1)  # 奇数化\nelse:\n    win = 51\n\ny_smooth = savgol_filter(y_denoised, window_length=win, polyorder=3)\n\n# === 6️⃣ 表示 ===\nplt.figure(figsize=(12,4))\nplt.plot(-y_smooth, color='steelblue')\nplt.title(f\"Lead II — Auto-corrected & Smoothed ECG (window={win})\")\nplt.xlabel(\"Time (pixels)\")\nplt.ylabel(\"Amplitude (a.u.)\")\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:19:50.808703Z","iopub.execute_input":"2025-11-29T06:19:50.809062Z","iopub.status.idle":"2025-11-29T06:19:51.124666Z","shell.execute_reply.started":"2025-11-29T06:19:50.809038Z","shell.execute_reply":"2025-11-29T06:19:51.123366Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =====================================================\n# ECG Waveform Centerline Extraction & CSV Export\n# =====================================================\n\nimport numpy as np\nimport pandas as pd\nimport cv2\nfrom scipy.signal import find_peaks\n\"\"\"\ndef extract_centerline_from_mask(mask, grid_dx\n                                 , grid_dy, scale_time=25.0, scale_mv=10.0):\n\n    # mask: 2D np.ndarray (0–255), 白線が波形のマスク画像\n    # grid_dx, grid_dy: グリッド間隔（ピクセル）\n    # scale_time: mm/s（ECG標準は25mm/s）\n    # scale_mv: mm/mV（標準は10mm/mV）\n\n    # return: DataFrame [time (s), voltage (mV)]\n    \n   \n    grid_dx = 22   # 画像の横スケールが広めなら少し大きく\n    grid_dy = 18   # 縦スケールが狭いなら少し小さく\n    \n    h, w = mask.shape\n    centerline_y = []\n    x_coords = np.arange(w)\n\n    # 列ごとに白ピクセルの平均y座標を求める\n    for x in range(w):\n        ys = np.where(mask[:, x] > 0)[0]\n        if len(ys) > 0:\n            center_y = np.mean(ys)\n        else:\n            center_y = np.nan\n        centerline_y.append(center_y)\n\n    centerline_y = np.array(centerline_y)\n\n    # 欠損を補間（線形補間＋前後補完）\n    nans = np.isnan(centerline_y)\n    if np.any(nans):\n        not_nan = np.logical_not(nans)\n        centerline_y[nans] = np.interp(np.flatnonzero(nans), np.flatnonzero(not_nan), centerline_y[not_nan])\n\n    # ---- ピクセル → 時間 / 電位 換算 ----\n    mm_per_px_x = 1.0 / grid_dx    # 横方向: 1px = (1/grid_dx) mm\n    mm_per_px_y = 1.0 / grid_dy    # 縦方向: 同上\n\n    # 時間軸 [秒]\n    time_sec = x_coords * mm_per_px_x / scale_time\n\n    # 電位 [mV]\n    # y方向は上がマイナスなので、反転して中央基準に\n    baseline = np.nanmedian(centerline_y)\n    mv = -(centerline_y - baseline) * mm_per_px_y / scale_mv\n\n    df = pd.DataFrame({\"time_sec\": time_sec, \"voltage_mV\": mv})\n    # yの欠損補完\n    df_wave[\"voltage_mV\"] = df_wave[\"voltage_mV\"].interpolate(method=\"linear\")\n\n    return df\n\"\"\"\n\ndef extract_centerline_from_mask(mask, grid_dx, grid_dy, scale_time, scale_mv):\n    \"\"\"二値マスクから中心線を抽出し、時系列[mV]に変換\"\"\"\n    H, W = mask.shape\n    m = (mask > 0).astype(np.float32)\n    y_idx = np.arange(H, dtype=np.float32)[:, None]\n    col_sum = m.sum(axis=0) + 1e-6\n    y_cog = (m * y_idx).sum(axis=0) / col_sum\n    y_smooth = pd.Series(y_cog).rolling(5, center=True, min_periods=1).mean().values\n\n    # ---- ピクセル → 時間/電圧 ----\n    baseline_row = H / 2.0\n    time_sec   = np.arange(W) / (grid_dx * (scale_time / 25.0))\n    voltage_mV = ((baseline_row - y_smooth) / grid_dy) * (1.0 / (scale_mv / 10.0))\n\n    # ---- DataFrame整形 ----\n    df_wave = pd.DataFrame({\"time_sec\": time_sec, \"voltage_mV\": voltage_mV})\n    if \"voltage_mV\" in df_wave.columns:\n        df_wave[\"voltage_mV\"] = df_wave[\"voltage_mV\"].interpolate(method=\"linear\")\n\n    return df_wave\n\n\nimport cv2\nimport numpy as np\nimport pandas as pd\nfrom scipy.ndimage import median_filter\nfrom scipy.signal import savgol_filter\nimport matplotlib.pyplot as plt\n\n\nimport numpy as np\nimport pandas as pd\nfrom scipy.signal import savgol_filter\nfrom scipy.ndimage import median_filter, gaussian_filter1d\n\ndef extract_centerline_from_mask_v4(mask, grid_dx, grid_dy,\n                                    scale_time=25.0, scale_mv=10.0):\n    \"\"\"\n    ECG Centerline Extraction v4\n    --------------------------------------------\n    目的：\n      - demo-submission の安定性を保持\n      - 過平滑化しない\n      - QRS の鋭さを壊さない\n      - baseline の崩壊を防ぐ\n    \"\"\"\n\n    H, W = mask.shape\n    m = (mask > 0).astype(np.uint8)\n\n    # ---- 1. median-based centerline (demo方式)\n    y_peak = []\n    for x in range(W):\n        pts = np.where(m[:, x] > 0)[0]\n        if len(pts) == 0:\n            y_peak.append(np.nan)\n        else:\n            y_peak.append(np.median(pts))\n\n    y_peak = pd.Series(y_peak).interpolate(limit_direction=\"both\").values\n\n    # ---- 2. 軽いノイズ除去（QRSを壊さない size=3）\n    y_denoised = median_filter(y_peak, size=3)\n\n    # ---- 3. 軽い平滑化（window=7, poly=2）\n    y_smooth = savgol_filter(y_denoised, window_length=7, polyorder=2)\n\n    # ---- 4. baseline correction（軽め Gaussian）\n    baseline = gaussian_filter1d(y_smooth, sigma=20)\n    y_final = y_smooth - baseline\n\n    # ---- 5. ピクセル → 実電圧変換\n    time_sec = np.arange(W) / (grid_dx * (scale_time / 25.0))\n\n    # 基準は中央値\n    baseline_row = np.nanmedian(y_smooth)\n    voltage_mV = ((baseline_row - y_smooth) / grid_dy) * (10.0 / scale_mv)\n\n    df = pd.DataFrame({\n        \"time_sec\": time_sec,\n        \"voltage_mV\": voltage_mV,\n        \"voltage_mV_smooth\": y_final\n    })\n\n    return df\n\ndef extract_centerline_v41(mask, grid_dx, grid_dy, scale_time, scale_mv):\n    \"\"\"\n    ECG centerline extraction v4.1\n    - v1 の安定性 + v4 の改良をミックス\n    - baseline 補正の強さを適正化\n    - 過平滑を防ぐ\n    \"\"\"\n\n    import numpy as np\n    import pandas as pd\n    from scipy.ndimage import median_filter, gaussian_filter1d\n    from scipy.signal import savgol_filter\n\n    H, W = mask.shape\n    m = (mask > 0).astype(np.float32)\n    y_idx = np.arange(H, dtype=np.float32)[:, None]\n\n    col_sum = m.sum(axis=0)\n    y_cog = np.zeros(W)\n    y_peak = np.zeros(W)\n\n    for x in range(W):\n        col = m[:, x]\n\n        # ---- 1. 有効な点がない列（欠損） ----\n        if col_sum[x] < 1e-3:\n            y_cog[x] = np.nan\n            y_peak[x] = np.nan\n            continue\n\n        # ---- 2. 重心（COG） ----\n        y_cog[x] = (col * y_idx[:, 0]).sum() / col_sum[x]\n\n        # ---- 3. 非ゼロ座標の中央値ピーク ----\n        nonzero_y = np.where(col > 0)[0]\n        if len(nonzero_y) > 0:\n            y_peak[x] = np.median(nonzero_y)\n        else:\n            y_peak[x] = y_cog[x]\n\n    # ---- 4. COG + median peak のハイブリッド ----\n    y_mix = np.nanmean(np.vstack([y_cog, y_peak]), axis=0)\n\n    # ---- 5. 線形補完で欠損を埋める ----\n    y_interp = pd.Series(y_mix).interpolate(limit_direction=\"both\").values\n\n    # ---- 6. 微小スパイク除去（v4 は size=5→過剰なので 3 に戻す）----\n    y_denoised = median_filter(y_interp, size=3)\n\n    # ---- 7. ガウス平滑（適度に）----\n    y_gauss = gaussian_filter1d(y_denoised, sigma=2)\n\n    # ---- 8. Savitzky–Golay（window=81 → 過平滑を防ぐ）----\n    y_smooth = savgol_filter(y_gauss, window_length=81, polyorder=3)\n\n    # ---- 9. ベースライン補正 ----\n    baseline = savgol_filter(y_smooth, window_length=401, polyorder=1)  # v4の301→弱めに401\n    y_corrected = y_smooth - baseline\n\n    # ---- 10. 最終 smoothing（軽め）----\n    y_final = savgol_filter(y_corrected, window_length=21, polyorder=2)\n\n    # ---- 11. ピクセル → 実電圧変換 ----\n    baseline_row = np.nanmedian(y_final)\n    time_sec = np.arange(W) / (grid_dx * (scale_time / 25.0))\n    voltage_mV = ((baseline_row - y_final) / grid_dy) * (1.0 / (scale_mv / 10.0))\n\n    df_wave = pd.DataFrame({\n        \"time_sec\": time_sec,\n        \"voltage_mV\": voltage_mV,\n        \"voltage_mV_smooth\": savgol_filter(voltage_mV, 21, 3)\n    })\n\n    return df_wave\n\n# =====================================================\n# Example usage\n# =====================================================\n\n# 例: wave_mask に波形マスク画像 (0-255) がある場合\n# grid_dx, grid_dy は前処理で推定した平均グリッド間隔を指定\n\"\"\"\ndf_wave = extract_centerline_from_mask(\n    mask=wave,          # waveは白線マスク画像\n    #grid_dx=20,         # 例：1mmあたり20px\n    #grid_dy=20,\n    grid_dx = 22,  # 画像の横スケールが広めなら少し大きく\n    grid_dy = 18,   # 縦スケールが狭いなら少し小さく\n\n    \n    scale_time=25.0,\n    scale_mv=10.0\n)\n\"\"\"\ndf_wave = extract_centerline_v41(\n    mask=binary, \n    grid_dx=20, \n    grid_dy=20, \n    scale_time=25.0, \n    scale_mv=10.0, \n    #debug_plot=True,        # 確認用に波形を表示\n    #overlay_image=gray      # 元画像を重ねて確認可能\n)\n\n\n\nfrom scipy.signal import savgol_filter\n\n# 波形平滑化を適用（ウィンドウ=9, 多項式=3 は経験的に安定）\ndf_wave[\"voltage_mV_smooth\"] = savgol_filter(df_wave[\"voltage_mV\"], window_length=9, polyorder=3)\n\n# プロットを差し替え\nplt.figure(figsize=(8,3))\nplt.plot(df_wave[\"time_sec\"], df_wave[\"voltage_mV_smooth\"], color=\"red\")\n#plt.title(f\"Smoothed ECG Trace: {lead_name}\")\nplt.title(f\"Smoothed ECG Trace: \")\nplt.xlabel(\"Time [sec]\"); plt.ylabel(\"Voltage [mV]\")\nplt.grid(True); plt.tight_layout(); plt.show()\n\n# 保存\nlead_name = \"V1\"  # リード名を適宜変更\nsave_path = f\"./ecg_trace_ver4.csv\"\ndf_wave.to_csv(save_path, index=False)\nprint(f\"Saved: {save_path}\")\ndisplay(df_wave.head())\n\n# 確認用プロット\nplt.figure(figsize=(8, 3))\nplt.plot(df_wave[\"time_sec\"], df_wave[\"voltage_mV\"], color=\"red\")\n#plt.title(f\"Extracted ECG Trace ({lead_name})\")\nplt.title(f\"Extracted ECG Trace\")\nplt.xlabel(\"Time [sec]\")\nplt.ylabel(\"Voltage [mV]\")\nplt.grid(True)\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-29T06:32:20.25323Z","iopub.execute_input":"2025-11-29T06:32:20.253546Z","iopub.status.idle":"2025-11-29T06:32:20.743899Z","shell.execute_reply.started":"2025-11-29T06:32:20.253525Z","shell.execute_reply":"2025-11-29T06:32:20.742798Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =====================================================\n# ECG Waveform Centerline Extraction & CSV Export\n# =====================================================\n\nimport numpy as np\nimport pandas as pd\nimport cv2\nfrom scipy.signal import find_peaks\n\"\"\"\ndef extract_centerline_from_mask(mask, grid_dx\n                                 , grid_dy, scale_time=25.0, scale_mv=10.0):\n\n    # mask: 2D np.ndarray (0–255), 白線が波形のマスク画像\n    # grid_dx, grid_dy: グリッド間隔（ピクセル）\n    # scale_time: mm/s（ECG標準は25mm/s）\n    # scale_mv: mm/mV（標準は10mm/mV）\n\n    # return: DataFrame [time (s), voltage (mV)]\n    \n   \n    grid_dx = 22   # 画像の横スケールが広めなら少し大きく\n    grid_dy = 18   # 縦スケールが狭いなら少し小さく\n    \n    h, w = mask.shape\n    centerline_y = []\n    x_coords = np.arange(w)\n\n    # 列ごとに白ピクセルの平均y座標を求める\n    for x in range(w):\n        ys = np.where(mask[:, x] > 0)[0]\n        if len(ys) > 0:\n            center_y = np.mean(ys)\n        else:\n            center_y = np.nan\n        centerline_y.append(center_y)\n\n    centerline_y = np.array(centerline_y)\n\n    # 欠損を補間（線形補間＋前後補完）\n    nans = np.isnan(centerline_y)\n    if np.any(nans):\n        not_nan = np.logical_not(nans)\n        centerline_y[nans] = np.interp(np.flatnonzero(nans), np.flatnonzero(not_nan), centerline_y[not_nan])\n\n    # ---- ピクセル → 時間 / 電位 換算 ----\n    mm_per_px_x = 1.0 / grid_dx    # 横方向: 1px = (1/grid_dx) mm\n    mm_per_px_y = 1.0 / grid_dy    # 縦方向: 同上\n\n    # 時間軸 [秒]\n    time_sec = x_coords * mm_per_px_x / scale_time\n\n    # 電位 [mV]\n    # y方向は上がマイナスなので、反転して中央基準に\n    baseline = np.nanmedian(centerline_y)\n    mv = -(centerline_y - baseline) * mm_per_px_y / scale_mv\n\n    df = pd.DataFrame({\"time_sec\": time_sec, \"voltage_mV\": mv})\n    # yの欠損補完\n    df_wave[\"voltage_mV\"] = df_wave[\"voltage_mV\"].interpolate(method=\"linear\")\n\n    return df\n\"\"\"\n\ndef extract_centerline_from_mask(mask, grid_dx, grid_dy, scale_time, scale_mv):\n    \"\"\"二値マスクから中心線を抽出し、時系列[mV]に変換\"\"\"\n    H, W = mask.shape\n    m = (mask > 0).astype(np.float32)\n    y_idx = np.arange(H, dtype=np.float32)[:, None]\n    col_sum = m.sum(axis=0) + 1e-6\n    y_cog = (m * y_idx).sum(axis=0) / col_sum\n    y_smooth = pd.Series(y_cog).rolling(5, center=True, min_periods=1).mean().values\n\n    # ---- ピクセル → 時間/電圧 ----\n    baseline_row = H / 2.0\n    time_sec   = np.arange(W) / (grid_dx * (scale_time / 25.0))\n    voltage_mV = ((baseline_row - y_smooth) / grid_dy) * (1.0 / (scale_mv / 10.0))\n\n    # ---- DataFrame整形 ----\n    df_wave = pd.DataFrame({\"time_sec\": time_sec, \"voltage_mV\": voltage_mV})\n    if \"voltage_mV\" in df_wave.columns:\n        df_wave[\"voltage_mV\"] = df_wave[\"voltage_mV\"].interpolate(method=\"linear\")\n\n    return df_wave\n\n\nimport cv2\nimport numpy as np\nimport pandas as pd\nfrom scipy.ndimage import median_filter\nfrom scipy.signal import savgol_filter\nimport matplotlib.pyplot as plt\n\ndef extract_centerline_from_mask_robust(mask, grid_dx=20, grid_dy=20,\n                                        scale_time=25.0, scale_mv=10.0,\n                                        debug_plot=False, overlay_image=None):\n    \"\"\"\n    ECGマスク画像から中心線を抽出してデジタル波形に変換する。\n    太線・ノイズ・欠損に強いロバスト版。\n    \"\"\"\n\n    # ---------- 1. 細線化 + ノイズ除去 ----------\n    mask_uint = (mask > 0).astype(np.uint8)\n    if hasattr(cv2.ximgproc, \"thinning\"):\n        thin = cv2.ximgproc.thinning(mask_uint)\n    else:\n        # ximgproc未インストール環境用 fallback\n        kernel = np.ones((3, 3), np.uint8)\n        thin = cv2.morphologyEx(mask_uint, cv2.MORPH_OPEN, kernel)\n\n    kernel = np.ones((3, 3), np.uint8)\n    thin = cv2.morphologyEx(thin, cv2.MORPH_OPEN, kernel)\n\n    h, w = thin.shape\n    y_positions = []\n\n    # ---------- 2. 列ごとに線の中心を取得 ----------\n    for x in range(w):\n        ys = np.where(thin[:, x] > 0)[0]\n        if len(ys) > 0:\n            y_positions.append(np.mean(ys))  # 中心（重心）\n        else:\n            y_positions.append(np.nan)\n\n    y_positions = np.array(y_positions)\n\n    # ---------- 3. 欠損補間 + スパイク除去 + 平滑化 ----------\n    y_series = pd.Series(y_positions).interpolate(limit_direction=\"both\")\n    y_med = median_filter(y_series, size=5)\n    y_smooth = savgol_filter(y_med, window_length=51, polyorder=3)\n\n    # ---------- 4. 時間軸と電圧スケールに変換 ----------\n    # ピクセル→秒, ピクセル→mV\n    time_sec = np.arange(len(y_smooth)) / grid_dx * (1.0 / scale_time)\n    voltage_mV = (h/2 - y_smooth) / grid_dy * (1.0 / scale_mv)\n\n    df = pd.DataFrame({\"time_sec\": time_sec, \"voltage_mV\": voltage_mV})\n\n    # ---------- 5. デバッグ表示 ----------\n    if debug_plot:\n        plt.figure(figsize=(10, 4))\n        if overlay_image is not None:\n            plt.imshow(overlay_image, cmap=\"gray\", origin=\"upper\")\n        plt.plot(np.arange(len(y_smooth)), y_smooth, color=\"red\", linewidth=1.5)\n        plt.title(\"Extracted ECG Centerline (Robust)\")\n        plt.xlabel(\"Time (px)\")\n        plt.ylabel(\"Y position (px)\")\n        plt.gca().invert_yaxis()\n        plt.tight_layout()\n        plt.show()\n\n        plt.figure(figsize=(8, 3))\n        plt.plot(df[\"time_sec\"], df[\"voltage_mV\"], color=\"red\")\n        plt.title(\"Smoothed ECG Trace (Scaled)\")\n        plt.xlabel(\"Time [sec]\")\n        plt.ylabel(\"Voltage [mV]\")\n        plt.grid(True)\n        plt.tight_layout()\n        plt.show()\n\n    return df\n\n\n# =====================================================\n# Example usage\n# =====================================================\n\n# 例: wave_mask に波形マスク画像 (0-255) がある場合\n# grid_dx, grid_dy は前処理で推定した平均グリッド間隔を指定\n\"\"\"\ndf_wave = extract_centerline_from_mask(\n    mask=wave,          # waveは白線マスク画像\n    #grid_dx=20,         # 例：1mmあたり20px\n    #grid_dy=20,\n    grid_dx = 22,  # 画像の横スケールが広めなら少し大きく\n    grid_dy = 18,   # 縦スケールが狭いなら少し小さく\n\n    \n    scale_time=25.0,\n    scale_mv=10.0\n)\n\"\"\"\ndf_wave = extract_centerline_v41(\n    mask=binary, \n    grid_dx=20, \n    grid_dy=20, \n    scale_time=25.0, \n    scale_mv=10.0, \n    #debug_plot=True,        # 確認用に波形を表示\n    #overlay_image=gray      # 元画像を重ねて確認可能\n)\n\n\n\nfrom scipy.signal import savgol_filter\n\n# 波形平滑化を適用（ウィンドウ=9, 多項式=3 は経験的に安定）\ndf_wave[\"voltage_mV_smooth\"] = savgol_filter(df_wave[\"voltage_mV\"], window_length=9, polyorder=3)\n\n# プロットを差し替え\nplt.figure(figsize=(8,3))\nplt.plot(df_wave[\"time_sec\"], df_wave[\"voltage_mV_smooth\"], color=\"red\")\n#plt.title(f\"Smoothed ECG Trace: {lead_name}\")\nplt.title(f\"Smoothed ECG Trace: \")\nplt.xlabel(\"Time [sec]\"); plt.ylabel(\"Voltage [mV]\")\nplt.grid(True); plt.tight_layout(); plt.show()\n\n# 保存\nlead_name = \"V1\"  # リード名を適宜変更\nsave_path = f\"./ecg_trace_.csv\"\ndf_wave.to_csv(save_path, index=False)\nprint(f\"Saved: {save_path}\")\ndisplay(df_wave.head())\n\n# 確認用プロット\nplt.figure(figsize=(8, 3))\nplt.plot(df_wave[\"time_sec\"], df_wave[\"voltage_mV\"], color=\"red\")\n#plt.title(f\"Extracted ECG Trace ({lead_name})\")\nplt.title(f\"Extracted ECG Trace\")\nplt.xlabel(\"Time [sec]\")\nplt.ylabel(\"Voltage [mV]\")\nplt.grid(True)\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:33:49.781282Z","iopub.execute_input":"2025-11-29T06:33:49.781587Z","iopub.status.idle":"2025-11-29T06:33:50.280238Z","shell.execute_reply.started":"2025-11-29T06:33:49.781567Z","shell.execute_reply":"2025-11-29T06:33:50.279308Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =====================================================\n# Recursive ECG Digitization Pipeline (subfolder aware)\n# =====================================================\n\nimport os\nimport glob\nimport cv2\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\n\n# 既存関数 extract_centerline_from_mask() をそのまま使用\n# from previous cell import extract_centerline_from_mask\n\n# =====================================================\n# 設定\n# =====================================================\nroot_dir = \"/kaggle/input/physionet-ecg-image-digitization/train\"             # trainフォルダをルートに指定\noutput_root = \"/kaggle/working/ecg_csv_outputs\"\n\nos.makedirs(output_root, exist_ok=True)\n\n# サブフォルダも含めて全てのpng/jpgを探索\nimage_files = sorted(glob.glob(os.path.join(root_dir, \"**\", \"*.png\"), recursive=True) +\n                     glob.glob(os.path.join(root_dir, \"**\", \"*.jpg\"), recursive=True))\n\nprint(f\"Found {len(image_files)} ECG images under {root_dir}.\")\n\n# =====================================================\n# パラメータ設定\n# =====================================================\ngrid_dx = 20\ngrid_dy = 10\nscale_time = 25.0\nscale_mv = 10.0\n\n\n\n# =====================================================\n# メインループ\n# =====================================================\nfor path in tqdm(image_files, desc=\"Processing ECG images\"):\n\n    # ---- 画像読み込み ----\n    img = cv2.imread(path)\n    if img is None:\n        print(f\"⚠️ Cannot open: {path}\")\n        continue\n    # ---- 画像読込 ----\n    \n    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n    thin = preprocess_ecg_image(gray, grid_kernel=30)\n\n\n    # ---- 波形抽出（簡易二値化） ----\n    blur = cv2.GaussianBlur(gray, (3, 3), 0)\n    _, binary = cv2.threshold(blur, 180, 255, cv2.THRESH_BINARY_INV)\n    kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (2, 2))\n    binary = cv2.morphologyEx(binary, cv2.MORPH_OPEN, kernel)\n\n    # ---- 波形中心線抽出 ----\n    df_wave = extract_centerline_v41(#extract_centerline_from_mask(\n        mask=binary,\n        grid_dx=grid_dx,\n        grid_dy=grid_dy,\n        scale_time=scale_time,\n        scale_mv=scale_mv\n    )\n\n    # ---- 出力先パスを決定（フォルダ階層を再現） ----\n    rel_path = os.path.relpath(path, root_dir)\n    rel_dir = os.path.dirname(rel_path)\n    out_dir = os.path.join(output_root, rel_dir)\n    os.makedirs(out_dir, exist_ok=True)\n\n    base_name = os.path.splitext(os.path.basename(path))[0]\n    out_csv = os.path.join(out_dir, f\"{base_name}.csv\")\n    df_wave.to_csv(out_csv, index=False)\n\n    # ---- サンプル画像を数件だけプロット ----\n    if np.random.rand() < 0.02:  # 全体の2%だけ可視化（多すぎ防止）\n        plt.figure(figsize=(8, 3))\n        plt.plot(df_wave[\"time_sec\"], df_wave[\"voltage_mV\"], color=\"red\")\n        plt.title(f\"Extracted: {rel_path}\")\n        plt.xlabel(\"Time [sec]\")\n        plt.ylabel(\"Voltage [mV]\")\n        plt.grid(True)\n        plt.tight_layout()\n        plt.show()\n\nprint(\"✅ Completed ECG extraction for all subfolders.\")\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:34:23.177672Z","iopub.execute_input":"2025-11-29T06:34:23.178015Z","iopub.status.idle":"2025-11-29T06:36:26.612588Z","shell.execute_reply.started":"2025-11-29T06:34:23.177992Z","shell.execute_reply":"2025-11-29T06:36:26.611146Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =====================================================\n# Multi-File ECG Digitization Pipeline/\\/\\\n# =====================================================\n\nimport os\nimport glob\nimport cv2\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom tqdm import tqdm\n\n# --- 既に定義済みの関数を使う想定 ---\n# from previous cell:\n#   extract_centerline_from_mask(mask, grid_dx, grid_dy, scale_time, scale_mv)\n\n# =====================================================\n# 設定\n# =====================================================\ninput_dir = \"/kaggle/input/physionet-ecg-image-digitization/train/1006867983/\"       # ECG画像を入れたフォルダを指定\noutput_dir = \"./ecg_csv_outputs\"  # CSV出力先\n\nos.makedirs(output_dir, exist_ok=True)\n\n# 画像ファイルの一覧取得（png, jpgどちらも対応）\nimage_files = sorted(glob.glob(os.path.join(input_dir, \"*.png\")) +\n                     glob.glob(os.path.join(input_dir, \"*.jpg\")))\n\nprint(f\"Found {len(image_files)} ECG images.\")\n\n# =====================================================\n# パラメータ設定\n# =====================================================\n#grid_dx = 20   # 1mmあたりのピクセル（横）\n#grid_dy = 20   # 1mmあたりのピクセル（縦）\ngrid_dx = 20 #22   # 画像の横スケールが広めなら少し大きく\ngrid_dy = 10#18   # 縦スケールが狭いなら少し小さく\nscale_time = 25.0  # 25mm/s\nscale_mv = 10.0    # 10mm/mV\n\n\n\n\n# =====================================================\n# メイン処理ループ\n# =====================================================\nfor path in tqdm(image_files, desc=\"Processing ECG images\"):\n\n    # ---- 画像読込 ----\n    img = cv2.imread(path)\n    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n\n    thin = preprocess_ecg_image(gray, grid_kernel=22)\n    \n    blur = cv2.GaussianBlur(thin, (3, 3), 0)\n    _, binary = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)\n    \n    kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (2, 2))\n    binary = cv2.morphologyEx(binary, cv2.MORPH_OPEN, kernel)\n    binary = cv2.bitwise_not(binary)\n    \"\"\"\n    thin = preprocess_ecg_image(gray, grid_kernel=30)\n\n    # ---- 波形マスク生成（簡易二値化）----\n    # 背景のグリッド除去と波形線抽出\n    blur = cv2.GaussianBlur(thin, (3, 3), 0)\n    _, binary = cv2.threshold(blur, 180, 255, cv2.THRESH_BINARY_INV)\n    kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (2, 2))\n    binary = cv2.morphologyEx(binary, cv2.MORPH_OPEN, kernel)\n    \"\"\"\n    # ---- 波形中心線抽出 ----\n    df_wave = extract_centerline_v41(\n        mask=binary,\n        grid_dx=grid_dx,\n        grid_dy=grid_dy,\n        scale_time=scale_time,\n        scale_mv=scale_mv\n    )\n\n    \n    from scipy.signal import savgol_filter\n    from scipy.ndimage import median_filter\n    \n    # --- 波形補間 + 平滑化 ---\n    # 欠損値補完\n    df_wave[\"voltage_mV\"] = df_wave[\"voltage_mV\"].interpolate(method=\"linear\")\n    \n    # 外れ値除去（3点メディアン）\n    df_wave[\"voltage_mV\"] = median_filter(df_wave[\"voltage_mV\"], size=3)\n    \n    # サヴィツキー・ゴレイで滑らか化\n    df_wave[\"voltage_mV_smooth\"] = savgol_filter(df_wave[\"voltage_mV\"], window_length=11, polyorder=3)\n    \n    # ベースライン補正（平均値をゼロに）\n    df_wave[\"voltage_mV_smooth\"] -= df_wave[\"voltage_mV_smooth\"].mean()\n    \n    # --- プロット ---\n    plt.figure(figsize=(8,3))\n    plt.plot(df_wave[\"time_sec\"], df_wave[\"voltage_mV_smooth\"], color=\"red\", lw=1)\n    #plt.title(f\"Smoothed ECG Trace: {lead_name}\")\n    plt.title(f\"Smoothed ECG Trace: \")\n    plt.xlabel(\"Time [sec]\"); plt.ylabel(\"Voltage [mV]\")\n    plt.grid(True); plt.tight_layout(); plt.show()\n\n\n\n\n    \n    # ---- ファイル名設定 ----\n    base = os.path.basename(path)\n    lead_name = os.path.splitext(base)[0]\n    save_csv = os.path.join(output_dir, f\"{lead_name}.csv\")\n    df_wave.to_csv(save_csv, index=False)\n\n    # ---- 確認プロット ----\n    plt.figure(figsize=(8, 3))\n    plt.plot(df_wave[\"time_sec\"], df_wave[\"voltage_mV\"], color=\"red\")\n    plt.title(f\"Extracted ECG Trace: \")\n    plt.xlabel(\"Time [sec]\")\n    plt.ylabel(\"Voltage [mV]\")\n    plt.grid(True)\n    plt.tight_layout()\n    plt.show()\n\nprint(\"✅ Completed ECG extraction for all files.\")\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:36:35.27569Z","iopub.execute_input":"2025-11-29T06:36:35.276044Z","iopub.status.idle":"2025-11-29T06:40:02.648424Z","shell.execute_reply.started":"2025-11-29T06:36:35.27602Z","shell.execute_reply":"2025-11-29T06:40:02.647256Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =====================================================\n# Combine 12-lead ECG CSVs into one aligned record ver2\n# =====================================================\n\nimport os\nimport pandas as pd\nimport numpy as np\nimport glob\nfrom scipy.signal import correlate\n\ndef align_and_combine_leads(csv_folder, max_shift=0.2, sample_rate=500):\n    \"\"\"\n    指定フォルダ内の12リードCSVを1つのDataFrameに整形\n    （時間軸・位相をそろえて統合）\n    \n    Parameters:\n      csv_folder: 各リードのCSVが入ったフォルダ\n      max_shift: 時間シフトの許容範囲 [秒]\n      sample_rate: サンプリング周波数 [Hz] （25mm/s → 約500Hz）\n    Returns:\n      pd.DataFrame: time列＋12リード列\n    \"\"\"\n    csv_files = sorted(glob.glob(os.path.join(csv_folder, \"*.csv\")))\n    if not csv_files:\n        print(f\"⚠️ No CSV files in {csv_folder}\")\n        return None\n\n    leads = {}\n    base_time = None\n\n    for f in csv_files:\n        lead_name = os.path.splitext(os.path.basename(f))[0].split(\"-\")[-1]  # ファイル末尾から誘導名推定\n        df = pd.read_csv(f)\n        if base_time is None:\n            base_time = df[\"time_sec\"].values\n        leads[lead_name] = df[\"voltage_mV\"].values\n\n    # 長さ合わせ\n    min_len = min(len(v) for v in leads.values())\n    for k in leads:\n        leads[k] = leads[k][:min_len]\n    base_time = base_time[:min_len]\n\n    # 時間シフト補正（II誘導を基準）\n    ref = leads.get(\"0001\", list(leads.values())[0])  # 最初のリードを基準に\n    aligned = {}\n    for name, sig in leads.items():\n        corr = correlate(ref - np.mean(ref), sig - np.mean(sig), mode=\"full\")\n        lag = np.argmax(corr) - len(sig) + 1\n        lag_time = lag / sample_rate\n        if abs(lag_time) > max_shift:\n            lag = 0  # 大きすぎるシフトは無視\n        aligned[name] = np.roll(sig, -lag)\n\n    # DataFrameへ統合\n    df_combined = pd.DataFrame({\"time_sec\": base_time})\n    for name, sig in aligned.items():\n        df_combined[name] = sig\n\n    return df_combined\n\n# =====================================================\n# Example: Combine one patient's leads\n# =====================================================\nrecord_folder = \"./ecg_csv_outputs/\"  # 対象フォルダを指定\ndf_12lead = align_and_combine_leads(record_folder)\n\nif df_12lead is not None:\n    display(df_12lead.head())\n    plt.figure(figsize=(10,6))\n    for col in df_12lead.columns[1:]:\n        plt.plot(df_12lead[\"time_sec\"], df_12lead[col] + 2*list(df_12lead.columns[1:]).index(col), label=col)\n    plt.title(f\"12-lead ECG Combined View ({os.path.basename(record_folder)})\")\n    plt.xlabel(\"Time [sec]\")\n    plt.ylabel(\"Relative Voltage [mV]\")\n    plt.legend(ncol=4, fontsize=8)\n    plt.grid(True)\n    plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:40:02.650163Z","iopub.execute_input":"2025-11-29T06:40:02.650437Z","iopub.status.idle":"2025-11-29T06:40:03.073907Z","shell.execute_reply.started":"2025-11-29T06:40:02.650418Z","shell.execute_reply":"2025-11-29T06:40:03.072751Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =====================================================\n# SNR Local Evaluation Function for PhysioNet ECG\n# =====================================================\n\nimport numpy as np\nfrom scipy.signal import correlate\n\ndef align_signals(pred, truth, sample_rate=500, max_shift=0.2):\n    \"\"\"\n    時間軸シフトを補正して2つの信号を整列。\n    最大±max_shift秒の範囲で相互相関を最大化。\n    \"\"\"\n    max_lag = int(max_shift * sample_rate)\n    corr = correlate(truth - np.mean(truth), pred - np.mean(pred), mode=\"full\")\n    lags = np.arange(-len(pred)+1, len(truth))\n    valid_idx = (lags >= -max_lag) & (lags <= max_lag)\n    lag = lags[valid_idx][np.argmax(corr[valid_idx])]\n    aligned_pred = np.roll(pred, lag)\n    return aligned_pred, lag / sample_rate\n\n\ndef compute_snr(pred, truth):\n    \"\"\"\n    縦オフセットを補正し、信号対雑音比（SNR, dB）を算出。\n    \"\"\"\n    pred = pred - np.mean(pred)\n    truth = truth - np.mean(truth)\n    noise = truth - pred\n    snr = 10 * np.log10(np.sum(truth**2) / np.sum(noise**2))\n    return snr\n\n\ndef evaluate_ecg_record(pred_df, truth_df, sample_rate=500):\n    \"\"\"\n    各リードごとにSNRを計算し、平均を返す。\n    \"\"\"\n    leads = [c for c in truth_df.columns if c != \"time_sec\" and c in pred_df.columns]\n    snr_list = []\n\n    for lead in leads:\n        pred = pred_df[lead].values\n        truth = truth_df[lead].values\n\n        # 長さ調整\n        min_len = min(len(pred), len(truth))\n        pred, truth = pred[:min_len], truth[:min_len]\n\n        # 整列・補正\n        aligned_pred, shift = align_signals(pred, truth, sample_rate=sample_rate)\n        snr = compute_snr(aligned_pred, truth)\n        snr_list.append(snr)\n\n        print(f\"{lead:>5s} | Shift: {shift:+.4f}s | SNR: {snr:.2f} dB\")\n\n    avg_snr = np.mean(snr_list)\n    print(f\"\\n✅ Mean SNR across leads: {avg_snr:.2f} dB\")\n    return avg_snr\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:40:25.954947Z","iopub.execute_input":"2025-11-29T06:40:25.955269Z","iopub.status.idle":"2025-11-29T06:40:25.967167Z","shell.execute_reply.started":"2025-11-29T06:40:25.955246Z","shell.execute_reply":"2025-11-29T06:40:25.965925Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nfrom scipy.signal import correlate\n\ndef _get_time_and_leads(df):\n    time_col = \"time_sec\" if \"time_sec\" in df.columns else (\"time\" if \"time\" in df.columns else None)\n    if time_col is None:\n        raise ValueError(\"時間列(time/time_sec)が見つかりません。\")\n    leads = [c for c in df.columns if c != time_col]\n    return time_col, leads\n\ndef _align_by_xcorr(pred, truth, sample_rate, max_shift=0.2):\n    max_lag = int(max_shift * sample_rate)\n    corr = correlate(truth - np.mean(truth), pred - np.mean(pred), mode=\"full\")\n    lags = np.arange(-len(pred)+1, len(truth))\n    sel = (lags >= -max_lag) & (lags <= max_lag)\n    lag = lags[sel][np.argmax(corr[sel])] if np.any(sel) else 0\n    return np.roll(pred, lag), lag / sample_rate\n\ndef evaluate_ecg_record_v2(pred_df, truth_df):\n    # 時間列取得\n    p_time, p_leads = _get_time_and_leads(pred_df)\n    t_time, t_leads = _get_time_and_leads(truth_df)\n    # 共通リード\n    common = [l for l in p_leads if l in t_leads]\n    if not common:\n        raise ValueError(f\"共通するリードがありません。\\n pred:{p_leads}\\n truth:{t_leads}\")\n    # サンプルレート推定（時間差の中央値から）\n    dt_pred = np.median(np.diff(pred_df[p_time].values))\n    sample_rate = int(round(1.0 / dt_pred)) if dt_pred > 0 else 500\n\n    snrs = []\n    for lead in common:\n        p = pred_df[lead].values\n        t = truth_df[lead].values\n        n = min(len(p), len(t))\n        p, t = p[:n], t[:n]\n        p_aligned, shift = _align_by_xcorr(p, t, sample_rate=sample_rate, max_shift=0.2)\n        # 縦オフセット除去\n        p_aligned = p_aligned - np.mean(p_aligned)\n        t = t - np.mean(t)\n        noise = t - p_aligned\n        if np.sum(noise**2) == 0:\n            snr = np.inf\n        else:\n            snr = 10*np.log10(np.sum(t**2)/np.sum(noise**2))\n        print(f\"{lead:>5s} | shift={shift:+.4f}s | SNR={snr:.2f} dB\")\n        snrs.append(snr)\n    mean_snr = float(np.mean(snrs))\n    print(f\"\\n✅ Mean SNR across leads: {mean_snr:.2f} dB (n={len(snrs)})\")\n    return mean_snr\n\n# 実行\n#mean_snr = evaluate_ecg_record_v2(pred_df, truth_df)\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:40:30.644173Z","iopub.execute_input":"2025-11-29T06:40:30.6445Z","iopub.status.idle":"2025-11-29T06:40:30.657529Z","shell.execute_reply.started":"2025-11-29T06:40:30.644478Z","shell.execute_reply":"2025-11-29T06:40:30.656522Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nfrom scipy.signal import correlate\n\ndef _infer_fs_from_time(df):\n    tcol = \"time_sec\" if \"time_sec\" in df.columns else (\"time\" if \"time\" in df.columns else None)\n    if tcol is None:\n        return 500  # フォールバック\n    dt = np.median(np.diff(df[tcol].values))\n    return int(round(1.0/dt)) if dt and np.isfinite(dt) and dt>0 else 500\n\ndef _align_xcorr_valid(pred, truth, fs, max_shift_s=0.2):\n    \"\"\"NaNを除外しつつ相互相関でシフト最適化（±max_shift_s）\"\"\"\n    max_lag = int(max_shift_s * fs)\n    # NaN除外\n    mask = np.isfinite(pred) & np.isfinite(truth)\n    if mask.sum() < 10:  # データ少なすぎ\n        return pred, 0.0\n    p = pred[mask] - np.nanmean(pred[mask])\n    q = truth[mask] - np.nanmean(truth[mask])\n    # 相関\n    corr = correlate(q, p, mode=\"full\")\n    lags = np.arange(-len(p)+1, len(q))\n    ok = (lags >= -max_lag) & (lags <= max_lag)\n    lag = lags[ok][np.argmax(corr[ok])] if ok.any() else 0\n    # シフト適用（全長に対して）\n    pred_shift = np.roll(pred, lag)\n    return pred_shift, lag / fs\n\ndef evaluate_ecg_record_robust(pred_df, truth_df, fs=None, max_shift_s=0.2, verbose=True):\n    \"\"\"\n    各リードで truth の「有限値が連続している先頭区間」だけを使い、\n    pred を同じ長さに切り出してから整列＆SNR計算。\n    \"\"\"\n    if fs is None:\n        fs = _infer_fs_from_time(pred_df)\n\n    p_time = \"time_sec\" if \"time_sec\" in pred_df.columns else \"time\"\n    t_time = \"time_sec\" if \"time_sec\" in truth_df.columns else \"time\"\n\n    pred_leads = [c for c in pred_df.columns if c != p_time]\n    truth_leads = [c for c in truth_df.columns if c != t_time]\n    common = [l for l in pred_leads if l in truth_leads]\n    if not common:\n        raise ValueError(f\"共通リードがありません: pred={pred_leads}, truth={truth_leads}\")\n\n    snr_list = []\n    for lead in common:\n        p = pred_df[lead].values.astype(float)\n        t = truth_df[lead].values.astype(float)\n\n        # ---- truth の「先頭から連続する有限区間」のみ使う ----\n        finite_truth = np.isfinite(t)\n        if not finite_truth.any():\n            if verbose: print(f\"{lead:>5s} | no finite truth -> skip\")\n            continue\n        # 先頭からの連続区間長\n        valid_len = np.argmax(~finite_truth) if (~finite_truth).any() else len(t)\n        t_seg = t[:valid_len]\n        p_seg = p[:valid_len] if len(p) >= valid_len else p  # predが短ければそのまま\n\n        # 長さ合わせ\n        n = min(len(p_seg), len(t_seg))\n        if n < 10:\n            if verbose: print(f\"{lead:>5s} | too short overlap -> skip\")\n            continue\n        p_seg = p_seg[:n]\n        t_seg = t_seg[:n]\n\n        # 整列（NaNを考慮した相関）\n        p_aligned, shift = _align_xcorr_valid(p_seg, t_seg, fs, max_shift_s=max_shift_s)\n\n        # オフセット除去 & NaN除去\n        mask = np.isfinite(p_aligned) & np.isfinite(t_seg)\n        if mask.sum() < 10:\n            if verbose: print(f\"{lead:>5s} | not enough finite -> skip\")\n            continue\n        p0 = p_aligned[mask] - np.mean(p_aligned[mask])\n        t0 = t_seg[mask] - np.mean(t_seg[mask])\n\n        noise = t0 - p0\n        denom = np.sum(noise**2)\n        if denom == 0:\n            snr = np.inf\n        else:\n            snr = 10*np.log10(np.sum(t0**2) / denom)\n\n        if verbose:\n            print(f\"{lead:>5s} | len={len(t0):4d} | shift={shift:+.4f}s | SNR={snr:.2f} dB\")\n        snr_list.append(snr)\n\n    if len(snr_list)==0:\n        if verbose: print(\"⚠️ 有効なリードがありません（NaNのみ）。\")\n        return np.nan\n\n    mean_snr = float(np.nanmean(snr_list))\n    if verbose:\n        print(f\"\\n✅ Mean SNR across leads: {mean_snr:.2f} dB (n={len(snr_list)})\")\n    return mean_snr\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:40:40.387602Z","iopub.execute_input":"2025-11-29T06:40:40.38795Z","iopub.status.idle":"2025-11-29T06:40:40.406786Z","shell.execute_reply.started":"2025-11-29T06:40:40.387928Z","shell.execute_reply":"2025-11-29T06:40:40.405558Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =====================================================\n# Example Usage\n# =====================================================\n\n# truth_df:  本当のECG時系列（PhysioNet提供のground truth）\n# pred_df:   Keiさんの抽出した波形（df_12leadなど）\n\n# 例: Check 同じ/kaggle/input/physionet-ecg-image-digitization/train/1006867983フォルダ内に正解ファイルがある場合\ntruth_path = \"/kaggle/input/physionet-ecg-image-digitization/train/1006867983/1006867983.csv\"\ntruth_df = pd.read_csv(truth_path)\n\n# Keiさんの推定結果（すでに作成済み）\npred_df = df_12lead.copy()\n\n# 典型的なマッピング例（必要に応じて調整）\nalias_to_num = {\n    \"I\":\"0001\",\"II\":\"0003\",\"III\":\"0004\",\n    \"aVR\":\"0005\",\"aVL\":\"0006\",\"aVF\":\"0009\",\n    \"V1\":\"0010\",\"V2\":\"0011\",\"V3\":\"0012\",\n    \"V4\":\"0013\",\"V5\":\"0014\",\"V6\":\"0015\"  # 末尾の番号は実データに合わせて変更する\n}\n\ndef normalize_lead_columns(df):\n    cols = list(df.columns)\n    # 時間列\n    time_col = \"time_sec\" if \"time_sec\" in cols else (\"time\" if \"time\" in cols else None)\n    # リード列の名前を数字系に寄せる\n    new = {}\n    for c in cols:\n        if c == time_col: \n            new[c] = \"time_sec\"\n        elif c in alias_to_num:\n            new[c] = alias_to_num[c]\n        else:\n            new[c] = c  # すでに 0001 等ならそのまま\n    df = df.rename(columns=new)\n    return df\n\ntruth_df = normalize_lead_columns(truth_df)\npred_df  = normalize_lead_columns(df_12lead.copy())\n\nprint(\"truth_df (norm) cols:\", [c for c in truth_df.columns if c!=\"time_sec\"][:15])\nprint(\"pred_df  (norm) cols:\", [c for c in pred_df.columns  if c!=\"time_sec\"][:15])\n\n#mean_snr = evaluate_ecg_record(pred_df, truth_df)\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:40:47.517926Z","iopub.execute_input":"2025-11-29T06:40:47.518507Z","iopub.status.idle":"2025-11-29T06:40:47.543533Z","shell.execute_reply.started":"2025-11-29T06:40:47.518478Z","shell.execute_reply":"2025-11-29T06:40:47.542626Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 欠損の多いtruth_dfをリードごとに分割\ntruth_signals = {}\nfor col in truth_df.columns:\n    if col.startswith(\"000\") or col in [\"I\", \"II\", \"III\", \"aVR\", \"aVL\", \"aVF\", \"V1\",\"V2\",\"V3\",\"V4\",\"V5\",\"V6\"]:\n        sig = truth_df[col].dropna().values\n        truth_signals[col] = sig\n        print(f\"{col}: {len(sig)} samples\")\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:40:52.808683Z","iopub.execute_input":"2025-11-29T06:40:52.809037Z","iopub.status.idle":"2025-11-29T06:40:52.818125Z","shell.execute_reply.started":"2025-11-29T06:40:52.809015Z","shell.execute_reply":"2025-11-29T06:40:52.817078Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## ステップ2: 仮の時間列を作る","metadata":{}},{"cell_type":"code","source":"sample_rate = 500\ntruth_signals_time = {\n    lead: np.arange(len(sig)) / sample_rate for lead, sig in truth_signals.items()\n}\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:40:57.491086Z","iopub.execute_input":"2025-11-29T06:40:57.491414Z","iopub.status.idle":"2025-11-29T06:40:57.497141Z","shell.execute_reply.started":"2025-11-29T06:40:57.491384Z","shell.execute_reply":"2025-11-29T06:40:57.496173Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## ステップ3: 時間列つき DataFrame 化","metadata":{}},{"cell_type":"code","source":"import pandas as pd\n\ntruth_dict = {}\nmax_len = max(len(sig) for sig in truth_signals.values())\ntime_sec = np.arange(max_len) / sample_rate\ntruth_dict[\"time_sec\"] = time_sec\n\nfor lead, sig in truth_signals.items():\n    pad = np.full(max_len, np.nan)\n    pad[:len(sig)] = sig\n    truth_dict[lead] = pad\n\ntruth_df_fixed = pd.DataFrame(truth_dict)\nprint(truth_df_fixed.shape)\ndisplay(truth_df_fixed.head())\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:40:59.528623Z","iopub.execute_input":"2025-11-29T06:40:59.52924Z","iopub.status.idle":"2025-11-29T06:40:59.547064Z","shell.execute_reply.started":"2025-11-29T06:40:59.529214Z","shell.execute_reply":"2025-11-29T06:40:59.546Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## ステップ4: 評価実行","metadata":{}},{"cell_type":"code","source":"# pred_df はすでに作成済み（12誘導整形したもの）\n# truth_df_fixed は「仮の time_sec を付けたもの」\n\nmean_snr = evaluate_ecg_record_robust(pred_df, truth_df_fixed)\n\n","metadata":{"execution":{"iopub.status.busy":"2025-11-29T06:41:05.26492Z","iopub.execute_input":"2025-11-29T06:41:05.265237Z","iopub.status.idle":"2025-11-29T06:41:05.279268Z","shell.execute_reply.started":"2025-11-29T06:41:05.265217Z","shell.execute_reply":"2025-11-29T06:41:05.278184Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## NaN行を除いて計算するSNR","metadata":{}},{"cell_type":"code","source":"def evaluate_ecg_record_safe(pred_df, truth_df):\n    \"\"\"\n    SNRをリードごとに計算し、NaNを自動的に無視する安全版。\n    \"\"\"\n    snr_results = []\n    for lead in truth_df.columns:\n        if lead not in pred_df.columns:\n            continue  # predに存在しない列はスキップ\n        \n        t = truth_df[lead].values\n        p = pred_df[lead].values\n\n        # NaNを除外（両方NaNでない部分だけを使う）\n        mask = ~np.isnan(t) & ~np.isnan(p)\n        if np.sum(mask) < 10:\n            continue  # サンプル数が少なすぎる場合は除外\n        \n        t_valid = t[mask]\n        p_valid = p[mask]\n\n        # SNR計算\n        signal_power = np.mean(t_valid ** 2)\n        noise_power = np.mean((t_valid - p_valid) ** 2)\n        snr_db = 10 * np.log10(signal_power / noise_power) if noise_power > 0 else np.inf\n\n        snr_results.append((lead, snr_db, len(t_valid)))\n\n    # 結果表示\n    for lead, snr_db, n in snr_results:\n        print(f\"{lead:>4} | len={n:4} | SNR={snr_db:6.2f} dB\")\n\n    mean_snr = np.mean([s for _, s, _ in snr_results]) if snr_results else np.nan\n    print(f\"\\n✅ Mean SNR (valid leads only): {mean_snr:.2f} dB (n={len(snr_results)})\")\n    return mean_snr\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-07T09:02:14.390423Z","iopub.execute_input":"2025-11-07T09:02:14.390775Z","iopub.status.idle":"2025-11-07T09:02:14.399732Z","shell.execute_reply.started":"2025-11-07T09:02:14.390754Z","shell.execute_reply":"2025-11-07T09:02:14.398386Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"mean_snr = evaluate_ecg_record_safe(pred_df,  truth_df_fixed)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-07T09:07:53.252135Z","iopub.execute_input":"2025-11-07T09:07:53.252505Z","iopub.status.idle":"2025-11-07T09:07:53.278239Z","shell.execute_reply.started":"2025-11-07T09:07:53.252482Z","shell.execute_reply":"2025-11-07T09:07:53.276743Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os, glob, pandas as pd\n\nrec_id = \"1006867983\"  # 比較したいレコードID\nfolder = f\"/kaggle/input/physionet-ecg-image-digitization/train/{rec_id}\"\nprint(\"CSV in folder:\", glob.glob(os.path.join(folder, \"*.csv\")))\n\ntruth_path = glob.glob(os.path.join(folder, \"*.csv\"))[0]  # 1つだけ想定\ntruth_df = pd.read_csv(truth_path)\nprint(\"truth_df columns:\", list(truth_df.columns)[:20])\ndisplay(truth_df.head(3))\n","metadata":{"execution":{"iopub.status.busy":"2025-11-07T09:06:05.260214Z","iopub.execute_input":"2025-11-07T09:06:05.260581Z","iopub.status.idle":"2025-11-07T09:06:05.301022Z","shell.execute_reply.started":"2025-11-07T09:06:05.260559Z","shell.execute_reply":"2025-11-07T09:06:05.299876Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"\n---\n## 🚀 Next Improvements\n- [ ] **Auto-rotation correction** — use average Hough angle to align grid horizontally  \n- [ ] **Grid interval estimation** — derive physical scale (ms, mV) from grid spacing  \n- [ ] **Multi-lead segmentation** — detect and separate lead regions (I, II, V1–V6)  \n- [ ] **Waveform denoising and R-peak detection** — prepare for signal reconstruction  \n\n---\n\n🧩 *Author: Keiko Itano (Lumo.Kei)*  \nKaggle Notebook prepared for the PhysioNet Digitization of ECG Images challenge.  \nFeedback and collaboration are welcome!\n","metadata":{}},{"cell_type":"markdown","source":"\n# ✅ v6 Smoothed Pipeline (Subpixel + Resample + Lowpass)\n\nこのセクションは、**サブピクセル中心線抽出**と**等間隔リサンプル＋Butterworth低域フィルタ**を追加した改良版です。  \n以前のパイプラインよりも階段状の段差が減り、滑らかなECG波形が得られます。\n","metadata":{}},{"cell_type":"code","source":"\nimport numpy as np, pandas as pd\nfrom scipy.signal import butter, filtfilt\nfrom scipy.interpolate import interp1d\n\ndef subpixel_centerline(mask, win=2):\n    \"\"\"列方向ごとの白画素重心からサブピクセル中心線を抽出\"\"\"\n    H, W = mask.shape\n    m = (mask > 0).astype(np.float32)\n    y_idx = np.arange(H, dtype=np.float32)[:, None]\n    col_sum = m.sum(axis=0) + 1e-6\n    y_cog = (m * y_idx).sum(axis=0) / col_sum\n    y_smooth = pd.Series(y_cog).rolling(2*win+1, center=True, min_periods=1).mean().values\n    return np.arange(W), y_smooth\n\ndef resample_and_lowpass(df, fs=500, fc=35):\n    \"\"\"等間隔再サンプリング＋Butterworth低域フィルタ（心電図向け）\"\"\"\n    t = df['time_sec'].values.astype(float)\n    v = df['voltage_mV'].values.astype(float)\n    t_new = np.arange(t.min(), t.max(), 1.0/fs)\n    f = interp1d(t, v, kind='cubic', fill_value='extrapolate')\n    v_new = f(t_new)\n    # 欠損があれば線形補間\n    finite_mask = np.isfinite(v_new)\n    if not finite_mask.all():\n        g = interp1d(t_new[finite_mask], v_new[finite_mask], kind='linear', fill_value='extrapolate')\n        v_new = g(t_new)\n    b, a = butter(N=2, Wn=fc/(fs/2), btype='low')\n    from scipy.signal import filtfilt\n    v_lp = filtfilt(b, a, v_new)\n    v_lp -= np.median(v_lp)\n    return pd.DataFrame({'time_sec': t_new, 'voltage_mV': v_lp})\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-07T03:22:38.020476Z","iopub.execute_input":"2025-11-07T03:22:38.021132Z","iopub.status.idle":"2025-11-07T03:22:38.034435Z","shell.execute_reply.started":"2025-11-07T03:22:38.021103Z","shell.execute_reply":"2025-11-07T03:22:38.03303Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# --- 改良版メインループ (v6) ---\nimport os, glob, cv2, matplotlib.pyplot as plt\nfrom tqdm import tqdm\n\ninput_dir = globals().get('input_dir', '/kaggle/input/physionet-ecg-image-digitization/train/1006867983/')\noutput_dir = globals().get('output_dir', './ecg_csv_outputs_v6')\nos.makedirs(output_dir, exist_ok=True)\n\nimage_files = sorted(glob.glob(os.path.join(input_dir, '*.png')) +\n                     glob.glob(os.path.join(input_dir, '*.jpg')))\nprint(f'[v6] Found {len(image_files)} images in', input_dir)\n\ngrid_dx = globals().get('grid_dx', 20)\ngrid_dy = globals().get('grid_dy', 10)#20)\nscale_time = globals().get('scale_time', 25.0)\nscale_mv = globals().get('scale_mv', 10.0)\n\nfor idx, path in enumerate(tqdm(image_files, desc='[v6] Processing')):\n    img = cv2.imread(path)\n    if img is None:\n        print('skip unreadable:', path); continue\n    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n\n    # 既存の関数がNotebookにある前提\n    thin = preprocess_ecg_image(gray, grid_kernel=22)\n\n    blur = cv2.GaussianBlur(thin, (3,3), 0)\n    _, binary = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)\n    kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (2,2))\n    binary = cv2.morphologyEx(binary, cv2.MORPH_OPEN, kernel)\n\n    xs, ys = subpixel_centerline(binary, win=2)\n\n    baseline_row = gray.shape[0] / 2.0\n    time_sec   = xs / (grid_dx * (scale_time / 25.0))\n    voltage_mV = ((baseline_row - ys) / grid_dy) * (1.0 / (scale_mv / 10.0))\n\n    import pandas as pd\n    df_wave = pd.DataFrame({'time_sec': time_sec, 'voltage_mV': voltage_mV})\n    df_wave['voltage_mV'] = df_wave['voltage_mV'].interpolate('linear')\n\n    df_wave = resample_and_lowpass(df_wave, fs=500, fc=35)\n\n    base = os.path.basename(path); lead_name = os.path.splitext(base)[0]\n    save_csv = os.path.join(output_dir, f'{lead_name}.csv')\n    df_wave.to_csv(save_csv, index=False)\n\n    if idx < 5:\n        plt.figure(figsize=(8,3))\n        plt.plot(df_wave['time_sec'], df_wave['voltage_mV'], color='red', lw=1)\n        plt.title(f'Smoothed ECG Trace (v6): {lead_name}')\n        plt.xlabel('Time [sec]'); plt.ylabel('Voltage [mV]')\n        plt.grid(True); plt.tight_layout(); plt.show()\n\nprint('✅ v6 pipeline finished. CSVs ->', output_dir)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-07T04:06:55.430006Z","iopub.execute_input":"2025-11-07T04:06:55.431246Z","iopub.status.idle":"2025-11-07T04:07:13.443674Z","shell.execute_reply.started":"2025-11-07T04:06:55.431186Z","shell.execute_reply":"2025-11-07T04:07:13.442232Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}