{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.12","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":31234,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Physionet 2025 - ECG images digitization using OpenCV, SciPy and NumPy in python\n\n**See [this](https://github.com/JackTheProgrammer/ECG-digitization-for-physionet-2025-competition) github repo for every step's in-depth knowledge**","metadata":{}},{"cell_type":"markdown","source":"## Step 1 — Imports & Debug","metadata":{}},{"cell_type":"code","source":"import os\nimport random\n\nimport cv2\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom scipy.signal import find_peaks, medfilt\n\n# Reproducibility\nrandom.seed(42)\nnp.random.seed(42)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:31:33.059193Z","iopub.execute_input":"2026-01-16T08:31:33.060179Z","iopub.status.idle":"2026-01-16T08:31:33.064897Z","shell.execute_reply.started":"2026-01-16T08:31:33.060139Z","shell.execute_reply":"2026-01-16T08:31:33.064011Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Step 2 — Image Loading","metadata":{}},{"cell_type":"code","source":"def load_image(path: str) -> np.ndarray:\n    \"\"\"\n    Safely load an ECG image.\n    Ensures BGR format for OpenCV compatibility.\n    \"\"\"\n    img = cv2.imread(path, cv2.IMREAD_UNCHANGED)\n    if img is None:\n        raise FileNotFoundError(f\"Could not read image: {path}\")\n\n    # Convert grayscale → BGR for uniform downstream handling\n    if img.ndim == 2:\n        img = cv2.cvtColor(img, cv2.COLOR_GRAY2BGR)\n    return img","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:31:37.266102Z","iopub.execute_input":"2026-01-16T08:31:37.266446Z","iopub.status.idle":"2026-01-16T08:31:37.271749Z","shell.execute_reply.started":"2026-01-16T08:31:37.266418Z","shell.execute_reply":"2026-01-16T08:31:37.270915Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Step 3 - Splitting the 12 leads","metadata":{}},{"cell_type":"code","source":"def get_precise_boundaries(band_img, num_leads=4, debug=False):\n    \"\"\"\n    Computes precise vertical lead boundaries using ECG grid geometry.\n    Anchors boundaries to bold (5 mm) grid lines.\n    \"\"\"\n\n    # 1. Convert to grayscale\n    gray = cv2.cvtColor(band_img, cv2.COLOR_BGR2GRAY)\n\n    # 2. Emphasize vertical grid lines\n    sobelx = cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=3)\n    sobelx = np.abs(sobelx)\n    sobelx /= (sobelx.max() + 1e-6)\n\n    # 3. Column-wise vertical energy\n    col_energy = sobelx.mean(axis=0)\n\n    # 4. Smooth profile\n    col_energy = cv2.GaussianBlur(\n        col_energy.reshape(1, -1),\n        (1, 31),\n        0\n    ).ravel()\n\n    # 5. Detect bold vertical grid lines\n    peaks, _ = find_peaks(\n        col_energy,\n        distance=30,\n        prominence=np.percentile(col_energy, 80)\n    )\n\n    if len(peaks) < 10:\n        raise RuntimeError(\"Insufficient grid lines detected\")\n\n    # 6. Grid spacing between bold lines (≈5 mm)\n    diffs = np.diff(peaks)\n    bold_grid_spacing_px = np.median(diffs)\n\n    # 7. Convert to pixels per mm\n    pixels_per_mm = bold_grid_spacing_px / 5.0\n\n    # 8. Lead width (ECG standard)\n    lead_width_px = int(62.5 * pixels_per_mm)\n\n    # 9. Anchor to first valid grid line\n    first_x = peaks[0]\n\n    boundaries = []\n    for i in range(num_leads):\n        x_start = int(first_x + i * lead_width_px)\n        x_end   = int(x_start + lead_width_px)\n        boundaries.append((x_start, x_end))\n\n    # 10. Debug visualization\n    if debug:\n        dbg = band_img.copy()\n        for i, (x0, x1) in enumerate(boundaries):\n            cv2.line(dbg, (x0, 0), (x0, dbg.shape[0]), (0, 255, 0), 2)\n            cv2.line(dbg, (x1, 0), (x1, dbg.shape[0]), (0, 0, 255), 2)\n            cv2.putText(\n                dbg, f\"L{i+1}\", (x0 + 5, 30),\n                cv2.FONT_HERSHEY_SIMPLEX, 0.8, (0, 0, 0), 2\n            )\n\n        plt.figure(figsize=(14, 3))\n        plt.imshow(dbg)\n        plt.title(\"Grid-anchored lead boundaries\")\n        plt.axis(\"off\")\n        plt.show()\n\n    return boundaries\n\ndef split_ecg_panels_refined(img_bgr, debug=False):\n    h, w = img_bgr.shape[:2]\n\n    # 1. Crop margins (Your empirical values)\n    top_crop = int(0.32 * h)\n    bottom_crop = int(0.03 * h)\n    core = img_bgr[top_crop : h - bottom_crop, :]\n    core_h, core_w = core.shape[:2]\n\n    # 2. Split into 4 horizontal bands\n\n    h1 = int(core_h)\n    band1 = core[0 : int(h1 * 0.34), :]\n\n    start_crop = int(band1.shape[0])\n    h2 = int(core_h - start_crop) - 55\n    band2 = core[int(start_crop - (start_crop * 0.3)) : h2, :]\n\n    start_crop = int(band2.shape[0]) + 200\n    h3 = int(h1 - int(h1 * 0.14))\n    band3 = core[start_crop : h3, :] # Approximation based on height3\n\n    start_crop = (int(band3.shape[0]) * 2) + 200\n    h4 = int(h1 - int(h1 * 0.03))\n    band4 = core[start_crop : h4, :]\n\n    bands = [band1, band2, band3, band4]\n\n    # 3. Detect precise vertical boundaries using the first band\n    # This identifies the 4 columns by finding the 3 separators\n    lead_boundaries = get_precise_boundaries(band1, debug=debug)\n    lead_layout = [\n        [\"I\",   \"aVR\", \"V1\", \"V4\"],\n        [\"II\",  \"aVL\", \"V2\", \"V5\"],\n        [\"III\", \"aVF\", \"V3\", \"V6\"]\n    ]\n\n    panels = {}\n\n    # 4. Extract 12 leads using dynamic boundaries\n    for r in range(3):\n        for c in range(4):\n            x_start, x_end = lead_boundaries[c]\n            # Sanity check — do NOT silently clip\n            if x_start < 0 or x_end > core_w:\n                raise ValueError(\n                    f\"Lead boundary out of bounds: ({x_start}, {x_end}) \"\n                    f\"for image width {core_w}\"\n                )\n\n            panels[lead_layout[r][c]] = bands[r][:, x_start:x_end]\n\n    # 5. Rhythm strip (Full Lead II)\n    panels[\"II_full\"] = bands[3]\n\n    if debug:\n        dbg = band1.copy()\n        for i, (x0, x1) in enumerate(lead_boundaries):\n            cv2.line(dbg, (x0, 0), (x0, dbg.shape[0]), (0, 255, 0), 2)\n            cv2.line(dbg, (x1, 0), (x1, dbg.shape[0]), (0, 0, 255), 2)\n            cv2.putText(\n                dbg, f\"L{i+1}\", (x0 + 5, 30),\n                cv2.FONT_HERSHEY_SIMPLEX, 0.8, (0, 0, 0), 2\n            )\n\n        plt.figure(figsize=(14, 3))\n        plt.imshow(dbg)\n        plt.title(\"Vertical lead boundaries (grid-anchored)\")\n        plt.axis(\"off\")\n        plt.show()\n\n    return panels","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:31:39.712250Z","iopub.execute_input":"2026-01-16T08:31:39.712606Z","iopub.status.idle":"2026-01-16T08:31:39.731991Z","shell.execute_reply.started":"2026-01-16T08:31:39.712577Z","shell.execute_reply":"2026-01-16T08:31:39.731063Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Step 4 - Mask acquiring by removing grid from splitted ECG Lead","metadata":{}},{"cell_type":"code","source":"def remove_ecg_grid(gray_img, debug=False):\n    \"\"\"\n    Removes ECG background grid using morphology only.\n    Preserves ECG waveform while eliminating horizontal and vertical grid lines.\n    \"\"\"\n\n    # 1. Normalize & binarize (signal becomes white)\n    gray = gray_img.copy()\n\n    # Light smoothing to reduce noise\n    gray = cv2.GaussianBlur(gray, (3, 3), 0)\n\n    # Adaptive threshold works better than fixed for scanned ECGs\n    binary = cv2.adaptiveThreshold(\n        gray,\n        255,\n        cv2.ADAPTIVE_THRESH_MEAN_C,\n        cv2.THRESH_BINARY_INV,\n        blockSize=21,\n        C=10\n    )\n\n    # 2. Detect vertical grid lines (long & straight)\n    vertical_kernel = cv2.getStructuringElement(\n        cv2.MORPH_RECT, (1, 35)\n    )\n    vertical_grid = cv2.morphologyEx(\n        binary,\n        cv2.MORPH_OPEN,\n        vertical_kernel,\n        iterations=1\n    )\n\n    # 3. Detect horizontal grid lines\n    horizontal_kernel = cv2.getStructuringElement(\n        cv2.MORPH_RECT, (35, 1)\n    )\n    horizontal_grid = cv2.morphologyEx(\n        binary,\n        cv2.MORPH_OPEN,\n        horizontal_kernel,\n        iterations=1\n    )\n\n    # 4. Combine grid components\n    grid_mask = cv2.bitwise_or(vertical_grid, horizontal_grid)\n\n    # 5. Remove grid from binary signal\n    signal_only = cv2.subtract(binary, grid_mask)\n\n    # 6. Restore ECG continuity\n    signal_kernel = cv2.getStructuringElement(\n        cv2.MORPH_ELLIPSE, (3, 3)\n    )\n\n    signal_only = cv2.morphologyEx(\n        signal_only,\n        cv2.MORPH_CLOSE,\n        signal_kernel,\n        iterations=1\n    )\n\n    signal_only = cv2.morphologyEx(\n        signal_only,\n        cv2.MORPH_OPEN,\n        signal_kernel,\n        iterations=1\n    )\n\n    # 7. Optional cleanup of tiny dots\n    signal_only = cv2.medianBlur(signal_only, 3)\n\n    if debug:\n        titles = [\n            \"Binary (signal + grid)\",\n            \"Vertical grid\",\n            \"Horizontal grid\",\n            \"Grid mask\",\n            \"Final signal mask\"\n        ]\n        imgs = [\n            binary,\n            vertical_grid,\n            horizontal_grid,\n            grid_mask,\n            signal_only\n        ]\n\n        plt.figure(figsize=(18, 4))\n        for i, (img, title) in enumerate(zip(imgs, titles)):\n            plt.subplot(1, 5, i + 1)\n            plt.imshow(img, cmap=\"gray\")\n            plt.title(title)\n            plt.axis(\"off\")\n        plt.tight_layout()\n        plt.show()\n\n    return signal_only","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:34:04.596093Z","iopub.execute_input":"2026-01-16T08:34:04.596438Z","iopub.status.idle":"2026-01-16T08:34:04.607512Z","shell.execute_reply.started":"2026-01-16T08:34:04.596409Z","shell.execute_reply":"2026-01-16T08:34:04.606637Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Step 5 - ECG trace extraction from binary mask","metadata":{}},{"cell_type":"code","source":"def extract_ecg_trace(signal_mask, debug=False):\n    \"\"\"\n    Extracts a 1D ECG trace from a binary signal mask.\n    Uses column-wise median to remain robust to thickness & noise.\n    \n    Returns:\n        x_coords : array of x pixel indices\n        y_coords : array of y pixel indices (NaN where missing)\n    \"\"\"\n\n    h, w = signal_mask.shape\n    y_trace = np.full(w, np.nan)\n\n    for x in range(w):\n        ys = np.where(signal_mask[:, x] > 0)[0]\n        if len(ys) > 0:\n            y_trace[x] = np.median(ys)\n\n    x_trace = np.arange(w)\n\n    if debug:\n        plt.figure(figsize=(12, 3))\n        plt.imshow(signal_mask, cmap=\"gray\")\n        plt.plot(x_trace, y_trace, color=\"red\", linewidth=1)\n        plt.title(\"Extracted ECG trace (overlay)\")\n        plt.axis(\"off\")\n        plt.show()\n\n    return x_trace, y_trace","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:34:09.312405Z","iopub.execute_input":"2026-01-16T08:34:09.312771Z","iopub.status.idle":"2026-01-16T08:34:09.320090Z","shell.execute_reply.started":"2026-01-16T08:34:09.312741Z","shell.execute_reply":"2026-01-16T08:34:09.319010Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Step 6 – time axis & resampling","metadata":{}},{"cell_type":"code","source":"# Step 6.1 - Extract 1D Trace From Mask\ndef mask_to_trace(mask):\n    \"\"\"\n    Convert binary signal mask into 1D trace (row index per column).\n    Returns float array with NaNs where signal is missing.\n    \"\"\"\n    h, w = mask.shape\n    trace = np.full(w, np.nan, dtype=np.float32)\n\n    for x in range(w):\n        ys = np.where(mask[:, x] > 0)[0]\n        if len(ys) > 0:\n            trace[x] = np.median(ys)\n\n    return trace\n\n# STEP 6.2 — Pixel → mV Conversion\ndef pixels_to_mV(trace_px, pixels_per_mm):\n    \"\"\"\n    Convert pixel-based trace to millivolts.\n    \"\"\"\n    pixels_per_mV = pixels_per_mm * 10.0\n\n    # Baseline = median ignoring NaNs\n    baseline = np.nanmedian(trace_px)\n\n    mV = (baseline - trace_px) / pixels_per_mV\n    return mV\n\n# STEP 6.3 — Resample to Exact `sig_len` (extracted signal has pixel width ≠ `sig_len`)\ndef resample_trace(trace, target_len):\n    \"\"\"\n    Resample ECG trace to target length using linear interpolation.\n    \"\"\"\n    x_old = np.linspace(0, 1, len(trace))\n    x_new = np.linspace(0, 1, target_len)\n\n    valid = np.isfinite(trace)\n\n    #  For the competition data this is always 10 seconds\n    if valid.sum() < 10:\n        return np.zeros(target_len, dtype=np.float32)\n\n    return np.interp(\n        x_new,\n        x_old[valid],\n        trace[valid]\n    ).astype(np.float32)\n\n# Step 6.4 - full wrap up of the entire cell 7\ndef digitize_ecg_lead(\n    panel_bgr,\n    pixels_per_mm,\n    sig_len,\n    debug=False\n):\n    \"\"\"\n    Full ECG lead digitization pipeline:\n    1. Remove grid from ECG panel\n    2. Extract 1D trace from binary mask\n    3. Convert pixels to mV\n    4. Resample to target signal length\n    5. Optional debug plots\n\n    Args:\n        `panel_bgr`: BGR image of ECG lead panel\n        `pixels_per_mm`: pixels per mm (from grid detection)\n        `sig_len`: target signal length in samples\n        `debug`: whether to show debug plots\n    \n    Returns:\n        `trace_resampled`: 1D numpy array of digitized ECG in mV\n    \"\"\"\n    gray = cv2.cvtColor(panel_bgr, cv2.COLOR_BGR2GRAY)\n\n    mask = remove_ecg_grid(gray, debug=False)\n\n    trace_px = mask_to_trace(mask)\n\n    trace_mV = pixels_to_mV(trace_px, pixels_per_mm)\n\n    trace_resampled = resample_trace(trace_mV, sig_len)\n\n    if debug:\n        plt.figure(figsize=(12, 3))\n        plt.plot(trace_resampled)\n        plt.title(\"Final Digitized ECG Lead\")\n        plt.xlabel(\"Samples\")\n        plt.ylabel(\"mV\")\n        plt.grid(True)\n        plt.show()\n\n    return trace_resampled","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:34:14.849940Z","iopub.execute_input":"2026-01-16T08:34:14.850252Z","iopub.status.idle":"2026-01-16T08:34:14.861681Z","shell.execute_reply.started":"2026-01-16T08:34:14.850225Z","shell.execute_reply":"2026-01-16T08:34:14.860677Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Step 7 - Pure calibration function","metadata":{}},{"cell_type":"code","source":"def calibrate_ecg_amplitude(\n    digitized_signal_pixels,\n    debug=False\n):\n    \"\"\"\n    Convert ECG trace from mm to mV.\n    ECG standard: 1 mm = 0.1 mV\n    \"\"\"\n\n    signal_mv = digitized_signal_pixels * 0.1  # Directly using 0.1 mV per pixel for calibration\n\n    if debug:\n        plt.figure(figsize=(12, 3))\n        plt.plot(signal_mv)\n        plt.title(\"Amplitude-Calibrated Lead (mV)\")\n        plt.xlabel(\"Samples\")\n        plt.ylabel(\"mV\")\n        plt.grid(True)\n        plt.show()\n\n    return signal_mv","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:34:20.474057Z","iopub.execute_input":"2026-01-16T08:34:20.474415Z","iopub.status.idle":"2026-01-16T08:34:20.480803Z","shell.execute_reply.started":"2026-01-16T08:34:20.474386Z","shell.execute_reply":"2026-01-16T08:34:20.479849Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Step 8 — Baseline Correction","metadata":{}},{"cell_type":"code","source":"def baseline_correct_ecg(signal, fs, window_sec=0.8):\n    \"\"\"\n    Gentle baseline wander removal.\n    Preserves ECG morphology and amplitude.\n    \"\"\"\n    window = int(fs * window_sec)\n    window = max(3, window | 1)  # odd length\n\n    baseline = np.convolve(signal, np.ones(window)/window, mode=\"same\")\n    corrected = signal - baseline\n\n    return corrected","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:34:24.998923Z","iopub.execute_input":"2026-01-16T08:34:24.999985Z","iopub.status.idle":"2026-01-16T08:34:25.006264Z","shell.execute_reply.started":"2026-01-16T08:34:24.999936Z","shell.execute_reply":"2026-01-16T08:34:25.004584Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Step 9 - Calculating actual image's lead's based pixels per mm","metadata":{}},{"cell_type":"code","source":"def estimate_pixels_per_mm(\n    lead_img: np.ndarray,\n    debug: bool = False\n) -> float:\n    \"\"\"\n    Robust estimation of pixels-per-mm from ECG grid.\n    Assumes standard ECG grid (1 mm small squares, 5 mm bold lines).\n    \"\"\"\n\n    # 1. Convert to grayscale\n    if lead_img.ndim == 3:\n        gray = cv2.cvtColor(lead_img, cv2.COLOR_BGR2GRAY)\n    else:\n        gray = lead_img.copy()\n\n    # 2. Suppress ECG trace (important!)\n    blur = cv2.GaussianBlur(gray, (5, 5), 0) # Light blur removes sharp waveform while preserving grid\n\n    # Adaptive threshold to isolate grid\n    binary = cv2.adaptiveThreshold(\n        blur, 255,\n        cv2.ADAPTIVE_THRESH_MEAN_C,\n        cv2.THRESH_BINARY_INV,\n        31, 10\n    )\n\n    # 3. Isolate vertical grid lines via morphology\n    h, w = binary.shape\n    vertical_kernel = cv2.getStructuringElement(\n        cv2.MORPH_RECT, (1, max(15, h // 40))\n    )\n\n    vertical_lines = cv2.morphologyEx(\n        binary, cv2.MORPH_OPEN, vertical_kernel, iterations=1\n    )\n\n    # 4. Column projection (grid-only)\n    col_sum = np.sum(vertical_lines, axis=0)\n\n    # Normalize for peak detection\n    col_sum = col_sum.astype(np.float32)\n    col_sum /= (np.max(col_sum) + 1e-6)\n\n    # 5. Detect grid peaks\n    peaks, _ = find_peaks(\n        col_sum,\n        distance=8,          # safe lower bound\n        prominence=0.2\n    )\n\n    if len(peaks) < 10:\n        raise RuntimeError(\"Grid detection failed — insufficient peaks\")\n\n    # 6. Estimate spacing\n    spacings = np.diff(peaks)\n\n    # Remove outliers\n    spacings = spacings[(spacings > 3) & (spacings < 40)]\n\n    # Smallest mode ≈ 1 mm\n    ppm = np.median(spacings)\n\n    # Handle bold 5-mm grid domination\n    if ppm > 18:\n        ppm /= 5.0\n\n    # 7. Debug visualization\n    if debug:\n        plt.figure(figsize=(12, 3))\n        plt.plot(col_sum, label=\"Grid projection\")\n        plt.plot(peaks, col_sum[peaks], \"x\", label=\"Grid peaks\")\n        plt.title(f\"Estimated pixels_per_mm = {ppm:.2f}\")\n        plt.legend()\n        plt.grid(True)\n        plt.show()\n\n    return float(ppm)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:34:30.517308Z","iopub.execute_input":"2026-01-16T08:34:30.517730Z","iopub.status.idle":"2026-01-16T08:34:30.530173Z","shell.execute_reply.started":"2026-01-16T08:34:30.517698Z","shell.execute_reply":"2026-01-16T08:34:30.528742Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Step 10 - Signal QC & Auto-Fix","metadata":{}},{"cell_type":"code","source":"# Step 10.1 – QC + Auto-fix Function\ndef ecg_qc_and_fix(\n    signal_mv: np.ndarray,\n    debug: bool = False\n) -> np.ndarray:\n    \"\"\"\n    Final ECG signal quality control and auto-fixing.\n    \"\"\"\n\n    sig = signal_mv.astype(np.float32).copy()\n\n    # 1. NaN / Inf handling\n    if not np.isfinite(sig).all():\n        idx = np.arange(len(sig))\n        valid = np.isfinite(sig)\n        sig[~valid] = np.interp(idx[~valid], idx[valid], sig[valid])\n\n    # 2. Flatline detection\n    if np.std(sig) < 1e-4:\n        raise ValueError(\"QC failed: flatline ECG detected\")\n\n    # 2.5. Impulse artifact suppression (time-aware)\n    # Large slope within very short time = digitization artifact\n    diff = np.abs(np.diff(sig))\n    # max_slope = 5.0  # mV\n    max_allowed = 0.5  # mV per sample (safe for digitized ECG)\n    artifact_idx = np.where(diff > max_allowed)[0]\n\n    if len(artifact_idx) > 0:\n        sig[artifact_idx + 1] = sig[artifact_idx]  # suppress spike\n\n    # 3. Polarity correction (ECG should be upright)\n    if np.abs(np.min(sig)) > np.abs(np.max(sig)):\n        sig = -sig\n\n    # 4. Amplitude sanity clamp (physiological)\n    sig = np.clip(sig, -5.0, 5.0)  # mV\n\n    # 5. Residual DC offset removal\n    sig -= np.median(sig)\n\n    # 6. Final smoothing (VERY light)\n    sig = medfilt(sig, kernel_size=3)\n\n    if debug:\n        plt.figure(figsize=(14, 4))\n        plt.plot(sig, label=\"After QC\", linewidth=1.5)\n        plt.plot(signal_mv, alpha=0.5, label=\"Before QC\")\n        plt.legend()\n        plt.title(\"Step 10 – QC-validated ECG signal\")\n        plt.xlabel(\"Samples\")\n        plt.ylabel(\"mV\")\n        plt.grid(True)\n        plt.show()\n\n    return sig","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:34:34.687426Z","iopub.execute_input":"2026-01-16T08:34:34.687834Z","iopub.status.idle":"2026-01-16T08:34:34.698154Z","shell.execute_reply.started":"2026-01-16T08:34:34.687805Z","shell.execute_reply":"2026-01-16T08:34:34.697132Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Step 11 — Finalization (Submission Guard)","metadata":{}},{"cell_type":"code","source":"import scipy.signal as s_signal\n\n# Local proxy SNR metric for ECG signals\ndef proxy_snr_ecg(ecg, fs) -> float:\n    \"\"\"\n    Metric-aware proxy SNR:\n    - Signal = band-limited ECG energy (0.5–40 Hz)\n    - Noise  = residual (high-frequency + artifacts)\n    \"\"\"\n\n    if len(ecg) < fs:\n        return -np.inf\n\n    # Bandpass filter (ECG band)\n    b, a = s_signal.butter(\n        4,\n        [0.5 / (fs / 2), 40 / (fs / 2)],\n        btype=\"band\"\n    )\n    ecg_band = s_signal.filtfilt(b, a, ecg)\n\n    # Residual = noise\n    noise = ecg - ecg_band\n\n    p_signal = np.mean(ecg_band ** 2)\n    p_noise = np.mean(noise ** 2) + 1e-12\n\n    return 10 * np.log10(p_signal / p_noise)\n\ndef metric_aware_postprocess(signal, lead):\n    \"\"\"\n    Final metric-aware scaling to maximize SNR.\n    \"\"\"\n\n    LEAD_RMS_TARGETS = {\n        \"I\":   0.035,\n        \"II\":  0.060,\n        \"III\": 0.030,\n        \"aVR\": 0.025,\n        \"aVL\": 0.025,\n        \"aVF\": 0.035,\n        \"V1\":  0.030,\n        \"V2\":  0.040,\n        \"V3\":  0.045,\n        \"V4\":  0.050,\n        \"V5\":  0.045,\n        \"V6\":  0.040,\n    }\n\n    target_rms = LEAD_RMS_TARGETS.get(lead, 0.04)\n\n    # Remove DC (metric already does this, but we help)\n    sig = signal - np.mean(signal)\n\n    rms = np.sqrt(np.mean(sig**2))\n    if rms < 1e-6:\n        return sig\n\n    sig = sig * (target_rms / rms)\n    return sig\n\ndef metric_safe_smoother(signal, fs):\n    \"\"\"\n    Removes high-frequency noise that the metric punishes,\n    while preserving QRS morphology.\n    \"\"\"\n\n    # Low-pass @ 25 Hz (metric-safe)\n    b, a = s_signal.butter(\n        4,\n        25 / (fs / 2),\n        btype=\"low\"\n    )\n    return s_signal.filtfilt(b, a, signal)\n\ndef finalize_ecg_signal(signal, fs, duration_sec):\n    expected_len = int(fs * duration_sec)\n\n    if len(signal) > expected_len:\n        return signal[:expected_len]\n    elif len(signal) < expected_len:\n        return np.pad(signal, (0, expected_len - len(signal)), mode=\"edge\")\n    # signal -= np.mean(signal)  # DC removal (safe)\n    return signal","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:34:38.546404Z","iopub.execute_input":"2026-01-16T08:34:38.546797Z","iopub.status.idle":"2026-01-16T08:34:38.559644Z","shell.execute_reply.started":"2026-01-16T08:34:38.546765Z","shell.execute_reply":"2026-01-16T08:34:38.558530Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## STEP 12 — Multi-Lead Batch Processing (FINAL)","metadata":{}},{"cell_type":"code","source":"# Step 12.1 — Canonical Lead Order\nECG_LEADS = {\n    \"I\", \"II\", \"III\",\n    \"aVR\", \"aVL\", \"aVF\",\n    \"V1\", \"V2\", \"V3\", \"V4\", \"V5\", \"V6\"\n}\n\nsnrs = []\n\n# Step 12.2(Per lead oriented) - Full multi-lead pipeline function\ndef process_all_ecg_leads(\n    img,\n    fs,\n    sig_len,\n    duration_sec\n):\n    \"\"\"\n    Process a full ECG image and return all 12 leads as calibrated mV signals.\n    \"\"\"\n\n    panels = split_ecg_panels_refined(img, debug=False)\n\n    all_leads = {}\n    errors = []\n\n    for lead in ECG_LEADS:\n        try:\n            panel = panels[lead]\n\n            # Step 6 — Digitization\n            signal = digitize_ecg_lead(\n                panel,\n                pixels_per_mm=estimate_pixels_per_mm(panel),\n                sig_len=sig_len,\n                debug=False\n            )\n\n            # Step 7 — Amplitude calibration (SNR-safe)\n            signal = calibrate_ecg_amplitude(signal)\n\n            # Step 8 — Baseline correction\n            signal = baseline_correct_ecg(signal, fs=fs)\n\n            # Step 9 — QC & stabilization\n            signal = ecg_qc_and_fix(signal)\n\n            # Metric-aware smoothing\n            signal = metric_safe_smoother(signal, fs)\n\n            # Step 10/11 — Final formatting\n            signal = finalize_ecg_signal(\n                signal,\n                fs=fs,\n                duration_sec=duration_sec\n            )\n\n            # Metric-aware postprocessing\n            signal = metric_aware_postprocess(signal, lead)\n\n            # Computing the signal's SNR locally\n            snr_score = proxy_snr_ecg(signal, fs)\n\n            snrs.append((lead, snr_score))\n\n            all_leads[lead] = signal\n\n        except Exception as e:\n            errors.append((lead, str(e)))\n\n    return all_leads, errors","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:34:44.549499Z","iopub.execute_input":"2026-01-16T08:34:44.549974Z","iopub.status.idle":"2026-01-16T08:34:44.561554Z","shell.execute_reply.started":"2026-01-16T08:34:44.549939Z","shell.execute_reply":"2026-01-16T08:34:44.560369Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Step 13 — Submission File creation","metadata":{}},{"cell_type":"code","source":"# Step 13.1 - Files picking and submission dirs\ncompetition_dataset_dir = '/kaggle/input/physionet-ecg-image-digitization'\ntest_folder_dir = f'{competition_dataset_dir}/test'\ntest_img1_dir = f'{test_folder_dir}/1053922973.png'\ntest_img2_dir = f'{test_folder_dir}/2352854581.png'\ntrain_folder_dir = f'{competition_dataset_dir}/train'\ntest_csv_dir = f'{competition_dataset_dir}/test.csv'\ntrain_csv_dir = f'{competition_dataset_dir}/train.csv'\noutput_dir = '/kaggle/working'\n\n# Step 13.2 - Submission file creation function\ndef build_submission_csv(\n    test_csv_path,\n    test_image_dir,\n    output_csv_path,\n    debug=False\n):\n    \"\"\"\n    Builds Kaggle-compatible submission.csv for PhysioNet ECG Image Digitization.\n    \"\"\"\n\n    test_df = pd.read_csv(test_csv_path)\n    submission_rows = []\n\n    # Group rows by base_id (VERY IMPORTANT)\n    for base_id, group in test_df.groupby(\"id\"):\n\n        base_id = str(base_id)\n        fs = int(group.iloc[0][\"fs\"])\n\n        img_path = os.path.join(test_image_dir, f\"{base_id}.png\")\n        img = cv2.imread(img_path)\n\n        if img is None:\n            raise FileNotFoundError(f\"Missing image: {img_path}\")\n\n        # Process ONCE per image\n        all_leads, errors = process_all_ecg_leads(\n            img,\n            fs=fs,\n            sig_len=fs * 10,\n            duration_sec=10,\n        )\n\n        if errors:\n            raise RuntimeError(f\"Errors for {base_id}: {errors}\")\n\n        # Create submission rows\n        for _, row in group.iterrows():\n            lead = row[\"lead\"]\n            n_rows = int(row[\"number_of_rows\"])\n\n            signal = all_leads[lead]\n\n            # Enforce exact expected length\n            if len(signal) < n_rows:\n                signal = np.pad(signal, (0, n_rows - len(signal)), mode=\"edge\")\n            else:\n                signal = signal[:n_rows]\n\n            for i, value in enumerate(signal):\n                submission_rows.append({\n                    \"id\": f\"{base_id}_{i}_{lead}\",\n                    \"value\": float(value)\n                })\n\n        if debug:\n            print(f\"[OK] {base_id}\")\n\n    submission_df = pd.DataFrame(submission_rows)\n    submission_df.to_csv(output_csv_path, index=False)\n\n    print(f\"\\n submission.csv written to: {output_csv_path}\")\n\n# Step 13.3 - Final submission file\nbuild_submission_csv(\n    test_csv_path=test_csv_dir,\n    test_image_dir=test_folder_dir,\n    output_csv_path=f\"{output_dir}/submission.csv\",\n    debug=True\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:34:49.223601Z","iopub.execute_input":"2026-01-16T08:34:49.224043Z","iopub.status.idle":"2026-01-16T08:34:50.221847Z","shell.execute_reply.started":"2026-01-16T08:34:49.224008Z","shell.execute_reply":"2026-01-16T08:34:50.220845Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"snrs_sorted = sorted(snrs, key=lambda x: x[1], reverse=True)\nfor lead, snr in snrs_sorted:\n    print(f\"Lead {lead}: Proxy SNR = {snr:.2f} dB\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-16T08:34:53.936949Z","iopub.execute_input":"2026-01-16T08:34:53.937266Z","iopub.status.idle":"2026-01-16T08:34:53.943536Z","shell.execute_reply.started":"2026-01-16T08:34:53.937238Z","shell.execute_reply":"2026-01-16T08:34:53.942561Z"}},"outputs":[],"execution_count":null}]}