{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"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":"\n# 🫀 ECG Digitization — Grid Detection & Waveform Extraction (Work in Progress)\n\nThis notebook explores the **digitization of ECG paper images** from the *PhysioNet - Digitization of ECG Images* competition.  \nThe aim is to extract ECG waveforms as numerical time-series data through image preprocessing, grid detection, and contour extraction.\n\n---\n\n### 📘 **Notebook Outline**\n1. **Overview & Data Setup**\n2. **ECG Image Visualization**\n3. **Grid Detection using Hough Transform**\n4. **Waveform Extraction via Thresholding**\n5. **Digitization: Pixel → Time-Amplitude Conversion**\n6. **CSV Export and Visualization**\n7. **Next Steps (Rotation, Scale, Multi-lead)**\n\n---\n\n### 💡 **Motivation**\nECG images often exist as scanned printouts, making them hard to analyze digitally.  \nBy detecting the grid and tracing the waveform line, we can reconstruct time-series ECG data for further analysis or ML modeling.\n\n> *Work in progress — shared for community feedback and educational purposes.*\n","metadata":{}},{"cell_type":"markdown","source":"# PhysioNet - Digitization of ECG Images 🫀\n### Starter Notebook (EDA + Grid & Waveform Detection)\n\n本ノートブックは PhysioNet の ECG 画像を対象に、波形データ化（デジタイズ）のための前処理テンプレートです。\n\n**内容:**\n1. データ構造の確認と可視化\n2. ECG画像の表示\n3. グリッド線の検出 (OpenCV Hough変換)\n4. 波形抽出 (二値化 + 輪郭抽出)\n5. ピクセル座標をデジタル波形に変換\n6. CSV出力\n\n---","metadata":{}},{"cell_type":"markdown","source":"## 1️⃣ パス設定と基本情報","metadata":{}},{"cell_type":"code","source":"\nimport os, cv2, numpy as np, pandas as pd, matplotlib.pyplot as plt\n\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","metadata":{"execution":{"iopub.status.busy":"2025-11-06T01:18:30.243810Z","iopub.execute_input":"2025-11-06T01:18:30.244105Z","iopub.status.idle":"2025-11-06T01:18:30.250986Z","shell.execute_reply.started":"2025-11-06T01:18:30.244087Z","shell.execute_reply":"2025-11-06T01:18:30.250015Z"},"trusted":true},"outputs":[],"execution_count":null},{"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-06T01:18:34.297103Z","iopub.execute_input":"2025-11-06T01:18:34.297684Z","iopub.status.idle":"2025-11-06T01:18:35.029478Z","shell.execute_reply.started":"2025-11-06T01:18:34.297658Z","shell.execute_reply":"2025-11-06T01:18:35.028548Z"},"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-06T01:18:55.836057Z","iopub.execute_input":"2025-11-06T01:18:55.836334Z","iopub.status.idle":"2025-11-06T01:18:56.694190Z","shell.execute_reply.started":"2025-11-06T01:18:55.836316Z","shell.execute_reply":"2025-11-06T01:18:56.693162Z"},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:20:56.777479Z","iopub.execute_input":"2025-11-06T01:20:56.777874Z","iopub.status.idle":"2025-11-06T01:20:56.787563Z","shell.execute_reply.started":"2025-11-06T01:20:56.777841Z","shell.execute_reply":"2025-11-06T01:20:56.786704Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"gray_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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:21:01.871191Z","iopub.execute_input":"2025-11-06T01:21:01.871950Z","iopub.status.idle":"2025-11-06T01:21:01.889838Z","shell.execute_reply.started":"2025-11-06T01:21:01.871918Z","shell.execute_reply":"2025-11-06T01:21:01.888590Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:21:05.370718Z","iopub.execute_input":"2025-11-06T01:21:05.371025Z","iopub.status.idle":"2025-11-06T01:21:06.063097Z","shell.execute_reply.started":"2025-11-06T01:21:05.371002Z","shell.execute_reply":"2025-11-06T01:21:06.062081Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:21:13.367043Z","iopub.execute_input":"2025-11-06T01:21:13.367928Z","iopub.status.idle":"2025-11-06T01:21:13.709966Z","shell.execute_reply.started":"2025-11-06T01:21:13.367896Z","shell.execute_reply":"2025-11-06T01:21:13.708935Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:21:18.947713Z","iopub.execute_input":"2025-11-06T01:21:18.948056Z","iopub.status.idle":"2025-11-06T01:21:19.295358Z","shell.execute_reply.started":"2025-11-06T01:21:18.948032Z","shell.execute_reply":"2025-11-06T01:21:19.294433Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:21:24.383958Z","iopub.execute_input":"2025-11-06T01:21:24.384732Z","iopub.status.idle":"2025-11-06T01:21:24.726407Z","shell.execute_reply.started":"2025-11-06T01:21:24.384692Z","shell.execute_reply":"2025-11-06T01:21:24.725455Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:21:29.766171Z","iopub.execute_input":"2025-11-06T01:21:29.766470Z","iopub.status.idle":"2025-11-06T01:21:31.486664Z","shell.execute_reply.started":"2025-11-06T01:21:29.766446Z","shell.execute_reply":"2025-11-06T01:21:31.485573Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:21:40.982550Z","iopub.execute_input":"2025-11-06T01:21:40.982876Z","iopub.status.idle":"2025-11-06T01:21:41.445981Z","shell.execute_reply.started":"2025-11-06T01:21:40.982855Z","shell.execute_reply":"2025-11-06T01:21:41.445057Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#　改善コード","metadata":{}},{"cell_type":"code","source":"for img_path in selected_images:\n    img = cv2.imread(img_path, cv2.IMREAD_GRAYSCALE)\n    thin = preprocess_ecg_image(img, grid_kernel=30)\n    # 以降、thin を使って輪郭検出・波形抽出へ\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:21:47.330368Z","iopub.execute_input":"2025-11-06T01:21:47.331149Z","iopub.status.idle":"2025-11-06T01:21:47.342104Z","shell.execute_reply.started":"2025-11-06T01:21:47.331113Z","shell.execute_reply":"2025-11-06T01:21:47.340956Z"}},"outputs":[],"execution_count":null},{"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:21:50.145293Z","iopub.execute_input":"2025-11-06T01:21:50.145635Z","iopub.status.idle":"2025-11-06T01:21:54.265559Z","shell.execute_reply.started":"2025-11-06T01:21:50.145604Z","shell.execute_reply":"2025-11-06T01:21:54.264577Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:22:06.516913Z","iopub.execute_input":"2025-11-06T01:22:06.517414Z","iopub.status.idle":"2025-11-06T01:22:06.607617Z","shell.execute_reply.started":"2025-11-06T01:22:06.517390Z","shell.execute_reply":"2025-11-06T01:22:06.606660Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:22:13.338200Z","iopub.execute_input":"2025-11-06T01:22:13.338523Z","iopub.status.idle":"2025-11-06T01:22:13.427415Z","shell.execute_reply.started":"2025-11-06T01:22:13.338480Z","shell.execute_reply":"2025-11-06T01:22:13.426526Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:22:20.550544Z","iopub.execute_input":"2025-11-06T01:22:20.550891Z","iopub.status.idle":"2025-11-06T01:22:20.858151Z","shell.execute_reply.started":"2025-11-06T01:22:20.550869Z","shell.execute_reply":"2025-11-06T01:22:20.857241Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:22:29.906898Z","iopub.execute_input":"2025-11-06T01:22:29.907540Z","iopub.status.idle":"2025-11-06T01:22:30.141217Z","shell.execute_reply.started":"2025-11-06T01:22:29.907506Z","shell.execute_reply":"2025-11-06T01:22:30.140561Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:22:34.980187Z","iopub.execute_input":"2025-11-06T01:22:34.980512Z","iopub.status.idle":"2025-11-06T01:22:35.154842Z","shell.execute_reply.started":"2025-11-06T01:22:34.980468Z","shell.execute_reply":"2025-11-06T01:22:35.153924Z"}},"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-06T01:22:44.145518Z","iopub.execute_input":"2025-11-06T01:22:44.145830Z","iopub.status.idle":"2025-11-06T01:22:44.352841Z","shell.execute_reply.started":"2025-11-06T01:22:44.145807Z","shell.execute_reply":"2025-11-06T01:22:44.351980Z"},"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-06T01:22:52.060339Z","iopub.execute_input":"2025-11-06T01:22:52.060678Z","iopub.status.idle":"2025-11-06T01:22:52.093062Z","shell.execute_reply.started":"2025-11-06T01:22:52.060655Z","shell.execute_reply":"2025-11-06T01:22:52.092290Z"},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:23:00.972044Z","iopub.execute_input":"2025-11-06T01:23:00.972344Z","iopub.status.idle":"2025-11-06T01:23:01.269661Z","shell.execute_reply.started":"2025-11-06T01:23:00.972321Z","shell.execute_reply":"2025-11-06T01:23:01.268703Z"}},"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    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\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,#k=wave_mask,          # 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\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\")\nplt.title(f\"Smoothed ECG Trace: {lead_name}\")\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_{lead_name}.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\")\nplt.title(f\"Extracted ECG Trace ({lead_name})\")\nplt.xlabel(\"Time [sec]\")\nplt.ylabel(\"Voltage [mV]\")\nplt.grid(True)\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:55:17.597538Z","iopub.execute_input":"2025-11-06T01:55:17.597848Z","iopub.status.idle":"2025-11-06T01:55:17.616796Z","shell.execute_reply.started":"2025-11-06T01:55:17.597824Z","shell.execute_reply":"2025-11-06T01:55:17.615583Z"}},"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 = 20\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_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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T08:45:50.436241Z","iopub.execute_input":"2025-11-05T08:45:50.436549Z","iopub.status.idle":"2025-11-05T14:31:23.089620Z","shell.execute_reply.started":"2025-11-05T08:45:50.436529Z","shell.execute_reply":"2025-11-05T14:31:23.085739Z"}},"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 = 22   # 画像の横スケールが広めなら少し大きく\ngrid_dy = 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_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    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.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: {lead_name}\")\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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T02:01:29.099609Z","iopub.execute_input":"2025-11-06T02:01:29.099931Z","iopub.status.idle":"2025-11-06T02:02:31.257783Z","shell.execute_reply.started":"2025-11-06T02:01:29.099910Z","shell.execute_reply":"2025-11-06T02:02:31.256488Z"}},"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/1006867983\"  # 対象フォルダを指定\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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-06T01:25:57.362850Z","iopub.status.idle":"2025-11-06T01:25:57.363111Z","shell.execute_reply.started":"2025-11-06T01:25:57.362973Z","shell.execute_reply":"2025-11-06T01:25:57.362984Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T14:58:26.056014Z","iopub.execute_input":"2025-11-05T14:58:26.056306Z","iopub.status.idle":"2025-11-05T14:58:26.072063Z","shell.execute_reply.started":"2025-11-05T14:58:26.056286Z","shell.execute_reply":"2025-11-05T14:58:26.070282Z"}},"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)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-04T08:10:58.116776Z","iopub.execute_input":"2025-11-04T08:10:58.117129Z","iopub.status.idle":"2025-11-04T08:10:58.151766Z","shell.execute_reply.started":"2025-11-04T08:10:58.117105Z","shell.execute_reply":"2025-11-04T08:10:58.150615Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-04T08:36:30.087125Z","iopub.execute_input":"2025-11-04T08:36:30.088610Z","iopub.status.idle":"2025-11-04T08:36:30.108649Z","shell.execute_reply.started":"2025-11-04T08:36:30.088561Z","shell.execute_reply":"2025-11-04T08:36:30.107614Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"truth_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-04T08:12:03.951502Z","iopub.execute_input":"2025-11-04T08:12:03.951822Z","iopub.status.idle":"2025-11-04T08:12:03.972328Z","shell.execute_reply.started":"2025-11-04T08:12:03.951800Z","shell.execute_reply":"2025-11-04T08:12:03.971544Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-04T08:20:56.175560Z","iopub.execute_input":"2025-11-04T08:20:56.176373Z","iopub.status.idle":"2025-11-04T08:20:56.201218Z","shell.execute_reply.started":"2025-11-04T08:20:56.176326Z","shell.execute_reply":"2025-11-04T08:20:56.199942Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-04T08:25:45.995915Z","iopub.execute_input":"2025-11-04T08:25:45.996348Z","iopub.status.idle":"2025-11-04T08:25:46.006834Z","shell.execute_reply.started":"2025-11-04T08:25:45.996321Z","shell.execute_reply":"2025-11-04T08:25:46.005600Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-04T08:26:44.071602Z","iopub.execute_input":"2025-11-04T08:26:44.072003Z","iopub.status.idle":"2025-11-04T08:26:44.077834Z","shell.execute_reply.started":"2025-11-04T08:26:44.071978Z","shell.execute_reply":"2025-11-04T08:26:44.076865Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-04T08:27:27.707427Z","iopub.execute_input":"2025-11-04T08:27:27.708103Z","iopub.status.idle":"2025-11-04T08:27:27.723861Z","shell.execute_reply.started":"2025-11-04T08:27:27.708074Z","shell.execute_reply":"2025-11-04T08:27:27.723103Z"}},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-04T08:37:19.214426Z","iopub.execute_input":"2025-11-04T08:37:19.215817Z","iopub.status.idle":"2025-11-04T08:37:19.226924Z","shell.execute_reply.started":"2025-11-04T08:37:19.215780Z","shell.execute_reply":"2025-11-04T08:37:19.225845Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-04T08:01:49.667478Z","iopub.execute_input":"2025-11-04T08:01:49.668363Z","iopub.status.idle":"2025-11-04T08:01:49.692068Z","shell.execute_reply.started":"2025-11-04T08:01:49.668334Z","shell.execute_reply":"2025-11-04T08:01:49.691213Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"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":{}}]}