{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"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":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"Forked from:\nhttps://www.kaggle.com/code/ambrosm/ecg-original-explained-baseline  \nMany thanks for the great base notebook.\n\nIn `convert_scanned_color`, I replaced the first and last `n_timesteps / 125` samples with their neighboring values to reduce edge artifacts caused by dark markings near the signal boundaries.","metadata":{}},{"cell_type":"markdown","source":"# Explained baseline for the *PhysioNet - Digitization of ECG Images* competition\n\nThis notebook contributes two innovations to the competition:\n1. It shows how to find the markers in a scanned ECG (object detection in class `MarkerFinder`).\n1. It visualizes how to extract the time series from images of types 3 and 11 (top-down plane sweep function `convert_scanned_color()`).\n\nWhy did I choose types 3 and 11? For two reasons: First, the scanned images have higher quality than the mobile photos, and they always have the same scale (80 pixels per mV). Second, the color makes it easy to distinguish the black ECG lines from the red gridlines.\n\n## Training data overview\n\nThere are 977 electrocardiograms (ids) in train, accounting for 84 GByte. Every electrocardiogram has 9 PNG files and one CSV file.\n\nThere are nine image types per ECG:\n- 0001 Original color ECG image generated by ECG-image-kit.\n- 0003 Image printed in color and scanned in color. → processed by `convert_scanned_color()`\n- 0004 Image printed in color and scanned in black and white.\n- 0005 Mobile photos of color printed images.\n- 0006 Mobile photos of ECGs on the screen of laptop.\n- 0009 Mobile photos of stained and soaked printed ECGs.\n- 0010 Mobile photos of printed ECGs with extensive damage.\n- 0011 Scans of printed ECG images with mold in color. → processed by `convert_scanned_color()`\n- 0012 Scans of printed ECG images with mold in black and white.\n\nThe sampling frequencies in train are 250, 256, 500, 512, 1000, 1025 per second.\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport cv2\nfrom glob import glob\nimport matplotlib.pyplot as plt\nfrom collections import defaultdict\nfrom tqdm import tqdm","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-27T15:49:33.736517Z","iopub.execute_input":"2025-10-27T15:49:33.736881Z","iopub.status.idle":"2025-10-27T15:49:36.764241Z","shell.execute_reply.started":"2025-10-27T15:49:33.736856Z","shell.execute_reply":"2025-10-27T15:49:36.763221Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Competition metric\n# From https://www.kaggle.com/code/metric/physionet-ecg-signal-extraction-metric\nfrom typing import Tuple\n\nimport numpy as np\nimport pandas as pd\n\nimport scipy.optimize\nimport scipy.signal\n\n\nLEADS = ['I', 'II', 'III', 'aVR', 'aVL', 'aVF', 'V1', 'V2', 'V3', 'V4', 'V5', 'V6']\nMAX_TIME_SHIFT = 0.2\nPERFECT_SCORE = 384\n\n\nclass ParticipantVisibleError(Exception):\n    pass\n\n\ndef compute_power(label: np.ndarray, prediction: np.ndarray) -> Tuple[float, float]:\n    if label.ndim != 1 or prediction.ndim != 1:\n        raise ParticipantVisibleError('Inputs must be 1-dimensional arrays.')\n    finite_mask = np.isfinite(prediction)\n    if not np.any(finite_mask):\n        raise ParticipantVisibleError(\"The 'prediction' array contains no finite values (all NaN or inf).\")\n\n    prediction[~np.isfinite(prediction)] = 0\n    noise = label - prediction\n    p_signal = np.sum(label**2)\n    p_noise = np.sum(noise**2)\n    return p_signal, p_noise\n\n\ndef compute_snr(signal: float, noise: float) -> float:\n    if noise == 0:\n        # Perfect reconstruction\n        snr = PERFECT_SCORE\n    elif signal == 0:\n        snr = 0\n    else:\n        snr = min((signal / noise), PERFECT_SCORE)\n    return snr\n\n\ndef align_signals(label: np.ndarray, pred: np.ndarray, max_shift: float = float('inf')) -> np.ndarray:\n    if np.any(~np.isfinite(label)):\n        raise ParticipantVisibleError('values in label should all be finite')\n    if np.sum(np.isfinite(pred)) == 0:\n        raise ParticipantVisibleError('prediction can not all be infinite')\n\n    # Initialize the reference and digitized signals.\n    label_arr = np.asarray(label, dtype=np.float64)\n    pred_arr = np.asarray(pred, dtype=np.float64)\n\n    label_mean = np.mean(label_arr)\n    pred_mean = np.mean(pred_arr)\n\n    label_arr_centered = label_arr - label_mean\n    pred_arr_centered = pred_arr - pred_mean\n\n    # Compute the correlation between the reference and digitized signals and locate the maximum correlation.\n    correlation = scipy.signal.correlate(label_arr_centered, pred_arr_centered, mode='full')\n\n    n_label = np.size(label_arr)\n    n_pred = np.size(pred_arr)\n\n    lags = scipy.signal.correlation_lags(n_label, n_pred, mode='full')\n    valid_lags_mask = (lags >= -max_shift) & (lags <= max_shift)\n\n    max_correlation = np.nanmax(correlation[valid_lags_mask])\n    all_max_indices = np.flatnonzero(correlation == max_correlation)\n    best_idx = min(all_max_indices, key=lambda i: abs(lags[i]))\n    time_shift = lags[best_idx]\n    start_padding_len = max(time_shift, 0)\n    pred_slice_start = max(-time_shift, 0)\n    pred_slice_end = min(n_label - time_shift, n_pred)\n    end_padding_len = max(n_label - n_pred - time_shift, 0)\n    aligned_pred = np.concatenate((np.full(start_padding_len, np.nan), pred_arr[pred_slice_start:pred_slice_end], np.full(end_padding_len, np.nan)))\n\n    def objective_func(v_shift):\n        return np.nansum((label_arr - (aligned_pred - v_shift)) ** 2)\n\n    if np.any(np.isfinite(label_arr) & np.isfinite(aligned_pred)):\n        results = scipy.optimize.minimize_scalar(objective_func, method='Brent')\n        vertical_shift = results.x\n        aligned_pred -= vertical_shift\n    return aligned_pred\n\n\ndef _calculate_image_score(group: pd.DataFrame) -> float:\n    \"\"\"Helper function to calculate the total SNR score for a single image group.\"\"\"\n\n    unique_fs_values = group['fs'].unique()\n    if len(unique_fs_values) != 1:\n        raise ParticipantVisibleError('Sampling frequency should be consistent across each ecg')\n    sampling_frequency = unique_fs_values[0]\n    if sampling_frequency != int(len(group[group['lead'] == 'II']) / 10):\n        raise ParticipantVisibleError('The sequence_length should be sampling frequency * 10s')\n    sum_signal = 0\n    sum_noise = 0\n    for lead in LEADS:\n        sub = group[group['lead'] == lead]\n        label = sub['value_true'].values\n        pred = sub['value_pred'].values\n\n        aligned_pred = align_signals(label, pred, int(sampling_frequency * MAX_TIME_SHIFT))\n        p_signal, p_noise = compute_power(label, aligned_pred)\n        sum_signal += p_signal\n        sum_noise += p_noise\n    return compute_snr(sum_signal, sum_noise)\n\n\ndef score(solution: pd.DataFrame, submission: pd.DataFrame, row_id_column_name: str) -> float:\n    \"\"\"\n    Compute the mean Signal-to-Noise Ratio (SNR) across multiple ECG leads and images for the PhysioNet 2025 competition.\n    The final score is the average of the sum of SNRs over different lines, averaged over all unique images.\n    Args:\n        solution: DataFrame with ground truth values. Expected columns: 'id' and one for each lead.\n        submission: DataFrame with predicted values. Expected columns: 'id' and one for each lead.\n        row_id_column_name: The name of the unique identifier column, typically 'id'.\n    Returns:\n        The final competition score.\n\n    Examples\n    --------\n    >>> import pandas as pd\n    >>> import numpy as np\n    >>> row_id_column_name = \"id\"\n    >>> solution = pd.DataFrame({'id': ['343_0_I', '343_1_I', '343_2_I', '343_0_III', '343_1_III','343_2_III','343_0_aVR', '343_1_aVR','343_2_aVR',\\\n    '343_0_aVL', '343_1_aVL', '343_2_aVL', '343_0_aVF', '343_1_aVF','343_2_aVF','343_0_V1', '343_1_V1', '343_2_V1','343_0_V2', '343_1_V2','343_2_V2',\\\n    '343_0_V3', '343_1_V3', '343_2_V3','343_0_V4', '343_1_V4', '343_2_V4', '343_0_V5', '343_1_V5','343_2_V5','343_0_V6', '343_1_V6','343_2_V6',\\\n    '343_0_II', '343_1_II','343_2_II', '343_3_II', '343_4_II', '343_5_II','343_6_II', '343_7_II','343_8_II','343_9_II','343_10_II','343_11_II'],\\\n    'fs': [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\\\n    'value':[0.1,0.3,0.4,0.6,0.6,0.4,0.2,0.3,0.4,0.5,0.2,0.7,0.2,0.3,0.4,0.8,0.6,0.7, 0.2,0.3,-0.1,0.5,0.6,0.7,0.2,0.9,0.4,0.5,0.6,0.7,0.1,0.3,0.4,\\\n    0.6,0.6,0.4,0.2,0.3,0.4,0.5,0.2,0.7,0.2,0.3,0.4]})\n    >>> submission = solution.copy()\n    >>> round(score(solution, submission, row_id_column_name), 4)\n    25.8433\n    >>> submission.loc[0, 'value'] = 0.9 # Introduce some noise\n    >>> round(score(solution, submission, row_id_column_name), 4)\n    13.6291\n    >>> submission.loc[4, 'value'] = 0.3 # Introduce some noise\n    >>> round(score(solution, submission, row_id_column_name), 4)\n    13.0576\n\n    >>> solution = pd.DataFrame({'id': ['343_0_I', '343_1_I', '343_2_I', '343_0_III', '343_1_III','343_2_III','343_0_aVR', '343_1_aVR','343_2_aVR',\\\n    '343_0_aVL', '343_1_aVL', '343_2_aVL', '343_0_aVF', '343_1_aVF','343_2_aVF','343_0_V1', '343_1_V1', '343_2_V1','343_0_V2', '343_1_V2','343_2_V2',\\\n    '343_0_V3', '343_1_V3', '343_2_V3','343_0_V4', '343_1_V4', '343_2_V4', '343_0_V5', '343_1_V5','343_2_V5','343_0_V6', '343_1_V6','343_2_V6',\\\n    '343_0_II', '343_1_II','343_2_II', '343_3_II', '343_4_II', '343_5_II','343_6_II', '343_7_II','343_8_II','343_9_II','343_10_II','343_11_II'],\\\n    'fs': [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\\\n    'value':[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]})\n    >>> round(score(solution, submission, row_id_column_name), 4)\n    -384\n    >>> submission = solution.copy()\n    >>> round(score(solution, submission, row_id_column_name), 4)\n    25.8433\n\n    >>> # test alignment\n    >>> label = np.array([0, 1, 2, 1, 0])\n    >>> pred = np.array([0, 1, 2, 1, 0])\n    >>> aligned = align_signals(label, pred)\n    >>> expected_array = np.array([0, 1, 2, 1, 0])\n    >>> np.allclose(aligned, expected_array, equal_nan=True)\n    True\n\n    >>> # Test 2: Vertical shift (DC offset) should be removed\n    >>> label = np.array([0, 1, 2, 1, 0])\n    >>> pred = np.array([10, 11, 12, 11, 10])\n    >>> aligned = align_signals(label, pred)\n    >>> expected_array = np.array([0, 1, 2, 1, 0])\n    >>> np.allclose(aligned, expected_array, equal_nan=True)\n    True\n\n    >>> # Test 3: Time shift should be corrected\n    >>> label = np.array([0, 0, 1, 2, 1, 0., 0.])\n    >>> pred = np.array([1, 2, 1, 0, 0, 0, 0])\n    >>> aligned = align_signals(label, pred)\n    >>> expected_array = np.array([np.nan, np.nan, 1, 2, 1, 0, 0])\n    >>> np.allclose(aligned, expected_array, equal_nan=True)\n    True\n    \n    >>> # Test 4: max_shift constraint prevents optimal alignment\n    >>> label = np.array([0, 0, 0, 0, 1, 2, 1]) # Peak is far\n    >>> pred = np.array([1, 2, 1, 0, 0, 0, 0])\n    >>> aligned = align_signals(label, pred, max_shift=10)\n    >>> expected_array = np.array([ np.nan, np.nan, np.nan, np.nan, 1, 2, 1])\n    >>> np.allclose(aligned, expected_array, equal_nan=True)\n    True\n\n    \"\"\"\n    for df in [solution, submission]:\n        if row_id_column_name not in df.columns:\n            raise ParticipantVisibleError(f\"'{row_id_column_name}' column not found in DataFrame.\")\n        if df['value'].isna().any():\n            raise ParticipantVisibleError('NaN exists in solution/submission')\n        if not np.isfinite(df['value']).all():\n            raise ParticipantVisibleError('Infinity exists in solution/submission')\n\n    submission = submission[['id', 'value']]\n    merged_df = pd.merge(solution, submission, on=row_id_column_name, suffixes=('_true', '_pred'))\n    merged_df['image_id'] = merged_df[row_id_column_name].str.split('_').str[0]\n    merged_df['row_id'] = merged_df[row_id_column_name].str.split('_').str[1].astype('int64')\n    merged_df['lead'] = merged_df[row_id_column_name].str.split('_').str[2]\n    merged_df.sort_values(by=['image_id', 'row_id', 'lead'], inplace=True)\n    image_scores = merged_df.groupby('image_id').apply(_calculate_image_score, include_groups=False)\n    return max(float(10 * np.log10(image_scores.mean())), -PERFECT_SCORE)","metadata":{"trusted":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2025-10-27T15:49:40.571032Z","iopub.execute_input":"2025-10-27T15:49:40.571473Z","iopub.status.idle":"2025-10-27T15:49:42.405939Z","shell.execute_reply.started":"2025-10-27T15:49:40.571450Z","shell.execute_reply":"2025-10-27T15:49:42.404556Z"},"jupyter":{"source_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Reading the metadata","metadata":{}},{"cell_type":"code","source":"train = pd.read_csv('/kaggle/input/physionet-ecg-image-digitization/train.csv')\ntest = pd.read_csv('/kaggle/input/physionet-ecg-image-digitization/test.csv')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-27T15:49:42.407375Z","iopub.execute_input":"2025-10-27T15:49:42.407997Z","iopub.status.idle":"2025-10-27T15:49:42.443550Z","shell.execute_reply.started":"2025-10-27T15:49:42.407963Z","shell.execute_reply":"2025-10-27T15:49:42.442440Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# The average ECG\n\nFor regression tasks, the average of the true labels is often a good baseline. Many public notebooks just submit this baseline with a little added noise for leaderboard scores between 0.14 and 0.18. The score difference compared to a zero submission is tiny.\n\nWe compute the mean time series per lead so that we can predict this mean for those image types which our model cannot handle.","metadata":{}},{"cell_type":"code","source":"def fit_mean_model(train, verbose=False):\n    \"\"\"Compute minima, maxima and means of the time series\"\"\"\n    mean_dict = defaultdict(list)\n    for idx, row in tqdm(train.iterrows(), total=len(train)):\n        labels = pd.read_csv(f'/kaggle/input/physionet-ecg-image-digitization/train/{row.id}/{row.id}.csv')\n        for lead in labels.columns:\n            values = labels[lead]\n            values = values[~values.isna()]\n            mean_dict[lead].append(values)\n    \n    for lead in mean_dict.keys():\n        # Resample every time series to 20000 samples\n        mean_dict[lead] = [\n            np.interp(np.linspace(0, len(values)-1, 20000), np.arange(len(values)), values)\n            for values in mean_dict[lead]\n        ]\n\n        # Stack all ECGs\n        mean_dict[lead] = np.stack(mean_dict[lead])\n\n        # Plot the mean ECG\n        if verbose:\n            m = mean_dict[lead].mean(axis=0)\n            # s = mean_dict[lead].std(axis=0)\n            plt.figure(figsize=(12, 2))\n            plt.title(f\"Mean curve for {lead}\")\n            plt.plot(m)\n            # plt.plot(m-s/30)\n            # plt.plot(m+s/30)\n            plt.axhline(0, color='gray')\n            plt.ylabel('mV')\n            plt.show()\n\n    return mean_dict\n\ndef validate_mean_model(val, mean_dict):\n    snr_list = []\n    for idx, row in tqdm(val.iterrows(), total=len(val)):\n        labels = pd.read_csv(f'/kaggle/input/physionet-ecg-image-digitization/train/{row.id}/{row.id}.csv')\n        # Evaluate the signal-to-noise ratio\n        sum_signal = 0\n        sum_noise = 0\n        for lead in labels.columns:\n            label = labels[lead]\n            label = label[~ label.isna()]\n            pred = mean_dict[lead].mean(axis=0)\n            pred = np.interp(np.linspace(0, 1, len(label)), np.linspace(0, 1, len(pred)), pred)\n            assert len(label) == len(pred)\n    \n            aligned_pred = align_signals(label, pred, int(row.fs * MAX_TIME_SHIFT))\n            p_signal, p_noise = compute_power(label, aligned_pred)\n            sum_signal += p_signal\n            sum_noise += p_noise\n    \n        snr = compute_snr(sum_signal, sum_noise)\n        # print(f\"{idx=:4d} id={row.id} SNR: {snr:.2f}\")\n        snr_list.append(snr)\n    \n    snr = np.array(snr_list).mean()\n    val_score = max(float(10 * np.log10(snr)), -PERFECT_SCORE)\n    print(f\"# Validation SNR for mean prediction: {snr:.2f} {val_score=:.2f}\")\n\n# Validate the mean model\ntrain_test_split_loc = 780\nmean_dict = fit_mean_model(train.iloc[:train_test_split_loc], verbose=True)\nvalidate_mean_model(train.iloc[train_test_split_loc:], mean_dict)\n\n# Refit the mean model to the full dataset\nmean_dict = fit_mean_model(train, verbose=False)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-27T15:49:44.765007Z","iopub.execute_input":"2025-10-27T15:49:44.765362Z","iopub.status.idle":"2025-10-27T15:51:24.807112Z","shell.execute_reply.started":"2025-10-27T15:49:44.765336Z","shell.execute_reply":"2025-10-27T15:51:24.805705Z"},"_kg_hide-output":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Finding the lead endpoints\n\nBefore decoding an image, it's good to know the coordinates of the 17 lead endpoints in the ECG. The following cell defines the class `MarkerFinder`, which determines these points. 13 points are found by the pattern matching function `cv2.matchTemplate()`; the right endpoints of the four lines are found as linear combinations of the other vectors.","metadata":{}},{"cell_type":"code","source":"class MarkerFinder:\n    \"\"\"This class finds the 13 markers in scanned ecg images and guesses the 4 line ends.\"\"\"\n    # From https://www.kaggle.com/code/ambrosm/ecg-original-explained-baseline\n    \n    def __init__(self, show_templates=False):\n        # Derive the templates from type 1 images\n        # np.max keeps the gridlines and markers and removes the ecg lines\n        ima = np.max([\n            cv2.imread('/kaggle/input/physionet-ecg-image-digitization/train/4292118763/4292118763-0001.png'),\n            cv2.imread('/kaggle/input/physionet-ecg-image-digitization/train/4289880010/4289880010-0001.png'),\n            cv2.imread('/kaggle/input/physionet-ecg-image-digitization/train/4284351157/4284351157-0001.png'),\n        ], axis=0)\n\n        # Template points in global coordinates of type 1 images\n        absolute_points = np.zeros((17, 2), dtype=int)\n        for i in range(3):\n            absolute_points[5 * i] = np.array([707 + 284 * i, 118]) # y, x\n            for j in range(1, 5):\n                absolute_points[5 * i + j] = np.array([707 + 284 * i, 118 + 492 * j])\n        absolute_points[5 * 3] = np.array([1535, 118])\n        absolute_points[5 * 3 + 1] = np.array([1535, 118 + 492 * 4])\n\n        # Top left corner of template rectangle\n        template_positions = [None] * 17\n        for i in range(len(absolute_points)):\n            if absolute_points[i][1] < 118 + 492 * 4:\n                if i % 5 == 0:\n                    template_positions[i] = (absolute_points[i][0] - 87, absolute_points[i][1] - 50) # y, x\n                else:\n                    template_positions[i] = (absolute_points[i][0] - 37, absolute_points[i][1] - 13)\n\n        # Height and width of the templates\n        template_sizes = np.array([(105, 60)] * 17) # height, width\n\n        # Convert the points to relative coordinates (inside the template)\n        template_points = [np.array([absolute_points[i][0] - template_positions[i][0],\n                                     absolute_points[i][1] - template_positions[i][1]])\n                           if template_positions[i] is not None\n                           else None\n                           for i in range(len(absolute_points))]\n\n        # Save the template matrices\n        templates = [None] * 17\n        for i in range(len(template_positions)):\n            if template_points[i] is not None:\n                template = (ima[template_positions[i][0]:template_positions[i][0]+template_sizes[i][0],\n                            template_positions[i][1]:template_positions[i][1]+template_sizes[i][1]])\n                templates[i] = template\n\n        # Plot the template matrices\n        if show_templates:\n            _, axs = plt.subplots(4, 4, figsize=(5, 7))\n            for i in range(len(template_positions)):\n                if template_points[i] is not None:\n                    template = templates[i].copy()\n                    cv2.rectangle(template,\n                                  (template_points[i][1]-1, template_points[i][0]-1),\n                                  (template_points[i][1]+1, template_points[i][0]+1), \n                                  [255, 0, 0], 2)\n                    axs[i // 5, i % 5].imshow(template)\n            for i in range(13, len(axs.ravel())):\n                axs.ravel()[i].axis('off')\n            plt.tight_layout()\n            plt.suptitle('The templates for the 13 markers', y=1.01)\n            plt.show()\n\n        self._absolute_points = absolute_points\n        self._template_positions = template_positions\n        self._template_sizes = template_sizes\n        self._template_points = template_points\n        self._templates = templates\n        \n    def find_markers(self, ima, warn=False, plot=False, title=''):\n        \"\"\"Return 17 markers as list of size-2 integer arrays (row, column)\n\n        Parameters:\n        ima: array of shape (1652, height, 3)\n        \"\"\"\n        \n        if ima.shape[0] != 1652:\n            raise ValueError(\"Implemented only for scanned images (image types 3, 4, 11, 12)\")\n\n        markers = [None] * 17\n\n        # Find 13 template-based markers\n        for j in range(len(self._templates)):\n            if self._template_points[j] is not None:\n                t = self._template_positions[j][0]-100\n                l = max(self._template_positions[j][1]-100, 0)\n                search_range = (ima[t:self._template_positions[j][0]+100+self._template_sizes[j][0],\n                                l:self._template_positions[j][1]+250+self._template_sizes[j][0]])\n                res = cv2.matchTemplate(search_range, self._templates[j], cv2.TM_CCOEFF)\n                min_val, max_val, min_loc, max_loc = cv2.minMaxLoc(res)\n    \n                top_left = max_loc\n                if warn and max_val < 3e7:\n                    bottom_right = (top_left[0] + self._templates[j].shape[1], top_left[1] + self._templates[j].shape[0])\n                    print(j, top_left, max_val)\n                    search_range = search_range.copy()\n                    cv2.rectangle(search_range, top_left, bottom_right, 0, 2)\n                    plt.imshow(search_range)\n                    plt.show()\n                markers[j] = np.array((t + top_left[1] + self._template_points[j][0], l + top_left[0] + self._template_points[j][1]))\n\n        # Guess the ends of the first three lines (can be outside the bounding box of the image)\n        for i in range(3):\n            m = markers[5 * i + 3] * 2 - markers[5 * i + 2]\n            markers[5 * i + 4] = m\n\n        # Guess the end of the fourth line (can be outside the bounding box of the image)\n        markers[16] = ((markers[14] * (284 + 260) - markers[9] * 260) / 284).astype(int)\n\n        if plot:\n            ima = ima.copy()\n            for m in markers:\n                if m is not None:\n                    cv2.rectangle(ima, (m[1]-40, m[0]-40), (m[1]+40, m[0]+40), (255, 0, 0), 2)\n            # plt.figure(figsize=(12, 8))\n            plt.imshow(ima)\n            plt.title(title)\n            plt.show()\n\n        return markers\n\n    # def baseline(self, i):\n    #     \"\"\"y coordinate of ith baseline in type 1 images\"\"\"\n    #     if i not in [0, 1, 2, 3]:\n    #         raise ValueError(\"i must be in [0, 1, 2, 3]\")\n    #     return self._absolute_points[5 * i][0]\n        \n    @staticmethod\n    def lead_info(lead):\n        \"\"\"Specify which markers mark the begin and the end of a lead.\"\"\"\n        begin, end = {\n            'I': (0, 1),\n            'II-subset': (5, 6),\n            'III': (10, 11),\n            'aVR': (1, 2),\n            'aVL': (6, 7),\n            'aVF': (11, 12),\n            'V1': (2, 3),\n            'V2': (7, 8),\n            'V3': (12, 13),\n            'V4': (3, 4),\n            'V5': (8, 9),\n            'V6': (13, 14),\n            'II': (15, 16),\n        }[lead]\n        return begin // 5, begin, end\n\n    def demo(self, ima, warn=False, title=''):\n        \"\"\"Plot the image with red markers\"\"\"\n        markers = self.find_markers(ima, warn, plot=True, title=title)\n\nmf = MarkerFinder(show_templates=False)\n\nima = cv2.imread('/kaggle/input/physionet-ecg-image-digitization/train/1026034238/1026034238-0011.png') # correct\nmf.demo(ima, warn=False, title='Scanned ECG with 17 line endpoints')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-27T15:51:26.638651Z","iopub.execute_input":"2025-10-27T15:51:26.639074Z","iopub.status.idle":"2025-10-27T15:51:28.125056Z","shell.execute_reply.started":"2025-10-27T15:51:26.639033Z","shell.execute_reply":"2025-10-27T15:51:28.123673Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Converting the images\n\nWe define the function `convert_scanned_color()`, which converts an image to twelve time series. This function fulfills the main task of the competition, but it works only for images of types 3 and 11. It does not yet generalize to images with a black grid in the background or to mobile photos.\n\nThe algorithm sweeps the image from top to bottom. The first black pixels detected during the sweep define the first line. The sweep continues over the white pixels below the line, and the next black pixels define the second line. An so on for the third and fourth line.\n\nAfter we have the four lines, we use the markers found by `MarkerFinder` to select the segments which form the 12 leads.","metadata":{}},{"cell_type":"code","source":"def find_line_by_topdown_sweep(ima):\n    \"\"\"Find the topmost black line in an image and remove it.\n\n    Parameters:\n    ima: 2d boolean image array (False = black, True = white), will be updated\n\n    Return values:\n    top: topmost black pixel in every column of the matrix\n    bottom: topmost white pixel in every column of the matrix below the topmost black pixel\n    \"\"\"\n    # Find the topmost black pixels\n    top = np.argmin(ima, axis=0) # topmost black (False) pixel per column; 0 if there are no black pixels\n\n    # Paint black everything above\n    mask = np.tile(np.arange(len(ima)).reshape(-1, 1), reps=(1, ima.shape[1]))\n    mask = mask >= top # True for lower part of image\n    ima &= mask # paint black whatever is above the line\n\n    # Find the topmost white pixels\n    bottom = np.argmax(ima, axis=0) # topmost white (True) pixel per column; 0 if there are no white pixels\n\n    # Paint white everything above\n    bottomx = np.maximum(bottom, np.median(top) + 100) # overpaints the letters\n    mask = np.tile(np.arange(len(ima)).reshape(-1, 1), reps=(1, ima.shape[1]))\n    mask = mask < bottomx # True for upper part of image\n    ima |= mask # paint white whatever is above the line\n    ima[:,:-1] |= mask[:,1:]\n    ima[:,1:] |= mask[:,:-1]\n\n    return top, bottom\n\ndef get_lead_from_top_bottom(tops, bottoms, lead, number_of_rows, markers):\n    \"\"\"Extract and resample one lead from an ECG line.\n    \n    Parameters:\n    tops: list of 4 arrays of shape (image_width, )\n    bottoms: list of 4 arrays of shape (image_width, )\n    lead: one of the 12 lead labels (string)\n    number_of_rows: number of samples required (int)\n    markers: 17 markers as list of size-2 integer arrays (row, column)\n    \"\"\"\n\n    # Select the markers and determine the baseline\n    line, begin, end = mf.lead_info(lead)\n    top = tops[line]\n    bottom = bottoms[line]\n    begin, end = markers[begin], markers[end]\n    baseline = np.linspace(begin[0], end[0], end[1] - begin[1])\n\n    pred0 = (top[begin[1]:end[1]] + bottom[begin[1]:end[1]]) / 2\n    if len(pred0) < len(baseline):\n        print(f\"smaller: {len(pred0)} < {len(baseline)}\")\n    baseline = baseline[:len(pred0)] # in case end is outside the image\n    pred = baseline - pred0\n\n    # Scale\n    pred /= 80 # 80 pixels = 1 mV\n\n    # Fix pixels obscured by the markers\n    if lead in ['aVR', 'aVL', 'aVF', 'V1', 'V2', 'V3', 'V4', 'V5', 'V6']:\n        # first four values can be obscured by the marker\n        pred[:4][pred[:4] > 0.2] = pred[4]\n    if lead in ['I', 'II-subset', 'III', 'aVR', 'aVL', 'aVF', 'V1', 'V2', 'V3']:\n        # last five pixels can be obscured by the marker\n        pred[-5:][pred[-5:] > 0.2] = pred[-6]\n    if lead in ['I', 'II-subset', 'III', 'II']:\n        # first pixel can be obscured by the marker\n        if 0.9 < pred[0] and pred[0] < 1.1 and pred[1] < 0.5:\n            pred[0] = pred[1]\n\n    # Upsample\n    pred = np.interp(np.linspace(0, 1, number_of_rows),\n                     np.linspace(0, 1, len(pred)),\n                     pred)\n    \n    # Fix implausible predictions\n    pred = np.where(np.abs(pred) <= 0.9, pred, 0)\n        \n    return pred\n\ndef convert_scanned_color(ima, markers, n_timesteps, verbose=False):\n    \"\"\"Convert a scanned color image (type 3 or 11) to 12 leads.\n\n    The function first extracts the four lines from the image. As\n    the four lines have nonnegligible width, we construct two lists:\n    - tops = y coordinates of the topmost black pixels in the lines\n    - bottoms = y coordinates of the topmost white pixels below the lines\n    Either list is a list of 4 arrays of shape (image_width, )\n\n    Parameters:\n    ima: 3-channel BGR image with height 1652 and width ≈2200.\n    markers: 17 markers as list of size-2 integer arrays (row, column)\n    n_timesteps: number of samples required per lead (dict)\n\n    Returns:\n    preds: dict with 12 time series\n    \"\"\"\n    # Crop the image and convert to black and white\n    # We use only the red channel (channel 2) so that the red gridlines disappear\n    # False = black, True = white\n    # The text at the top of the image is discarded.\n    crop_top = 400\n    ima = ima[crop_top:, :, 2] > 160\n\n    # Denoise single and double black pixels\n    iima = ima.astype(np.uint8)\n    ima = (iima[:-2, :-2] + iima[:-2, 1:-1] + iima[:-2, 2:]\n           + iima[1:-1, :-2] + iima[1:-1, 1:-1] + iima[1:-1, 2:]\n           + iima[2:, :-2] + iima[2:, 1:-1] + iima[2:, 2:]) >= 7\n\n    # Plot the denoised black-and-white image\n    if verbose:\n        plt.figure(figsize=(6, 4))\n        plt.imshow(ima)\n        plt.title('Denoised black-and-white')\n        # plt.savefig('ima.png')\n        plt.show()\n    \n    # Find the four lines\n    tops, bottoms = [], []\n    for i in range(4):\n        top, bottom = find_line_by_topdown_sweep(ima)\n        tops.append(top)\n        bottoms.append(bottom)\n\n    # Transform to global coordinates\n    tops = [t + crop_top for t in tops]\n    bottoms = [b + crop_top for b in bottoms]\n\n    # Extract the twelve leads from the four lines\n    # (as the first part of II is duplicated, we extract it twice\n    # and take the average)\n    n_timesteps['II-subset'] = n_timesteps['I']\n    preds = {}\n    for i, lead in enumerate(LEADS + ['II-subset']):\n        pred = get_lead_from_top_bottom(tops, bottoms, lead, n_timesteps[lead], markers)\n        preds[lead] = pred\n\n    preds['II'][:len(preds['II-subset'])] = (preds['II'][:len(preds['II-subset'])] + preds['II-subset']) / 2\n    del preds['II-subset']\n\n    # === ★ ===\n    for lead, sig in preds.items():\n        sig = np.asarray(sig)\n        n = len(sig)\n        if lead not in n_timesteps:\n            continue\n        replace_len = int(n_timesteps[lead] / 125)\n\n        if replace_len > 0 and replace_len * 2 < n:\n            # 始端部分を置換\n            sig[:replace_len] = sig[replace_len]\n            # 終端部分を置換\n            sig[-replace_len:] = sig[-replace_len - 1]\n\n        preds[lead] = sig\n    # === ★ ===\n\n    # Apply Einthoven's law\n    apply_einthoven(preds)\n\n    return preds\n\ndef apply_einthoven(preds):\n    \"\"\"Apply Einthoven's law to improve the predictions.\n    \n    Parameters:\n    pred: dict of time series, will be updated\n    \"\"\"\n    residual = preds['I'] + preds['III'] - preds['II'][:len(preds['III'])]\n    correction = residual / 3\n    preds['I'] -= correction\n    preds['III'] -= correction\n    preds['II'][:len(preds['III'])] += correction\n    \n    residual = preds['aVR'] + preds['aVL'] + preds['aVF']\n    correction = residual / 3\n    preds['aVR'] -= correction\n    preds['aVL'] -= correction\n    preds['aVF'] -= correction\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-27T16:02:02.236678Z","iopub.execute_input":"2025-10-27T16:02:02.238151Z","iopub.status.idle":"2025-10-27T16:02:02.271101Z","shell.execute_reply.started":"2025-10-27T16:02:02.238097Z","shell.execute_reply":"2025-10-27T16:02:02.270030Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Validation\n\nWe convert a few training images, plot the output and compute the signal-to-noise ratio. If you look at the diagrams closely, you'll easily get ideas for improvements.\n\nValidation is important: If you only get feedback from five submissions per day, your progress will be much too slow.","metadata":{}},{"cell_type":"code","source":"def validate_algorithm(train, image_types, convert):\n    \"\"\"Convert a few training images, plot the output and compute the signal-to-noise ratio\"\"\"\n    snr_list = []\n    index_list = []\n    is_first_ecg = True\n    for idx, row in train.iterrows():\n        # print(idx, row.id)\n        labels = pd.read_csv(f'/kaggle/input/physionet-ecg-image-digitization/train/{row.id}/{row.id}.csv')\n        png_paths = sorted(glob(f'/kaggle/input/physionet-ecg-image-digitization/train/{row.id}/{row.id}-*.png'))\n        for path in png_paths:\n            img_type = int(path[-8:-4])\n            # ima = cv2.imread(path)\n            # shape = ima.shape\n            # assert (img_type == 1) <= (shape == (1700, 2200, 3)) # 200 pixels per inch on Letter paper\n            # assert (img_type == 3) <= (shape[0] == 1652)\n            # assert (img_type == 4) <= (shape[0] == 1652)\n            # assert (img_type == 5) <= (len(shape) == 3 and shape[2] == 3)\n            # assert (img_type == 6) <= ((shape == (4000, 3000, 3) or (shape == (3000, 4000, 3))))\n            # assert (img_type == 9) <= ((shape == (3024, 4032, 3) or (shape == (4032, 3024, 3))))\n            # assert (img_type == 10) <= ((shape == (3024, 4032, 3) or (shape == (4032, 3024, 3))))\n            # assert (img_type == 11) <= (shape[0] == 1652)\n            # assert (img_type == 12) <= (shape[0] == 1652)\n            if img_type in image_types:\n                ima = cv2.imread(path)\n\n                # Find the 17 line endpoints\n                markers = mf.find_markers(ima, plot=is_first_ecg, title='Image with 17 markers')\n\n                # Convert the image to 12 leads\n                n_timesteps = {lead: (~ labels[lead].isna()).sum() for lead in LEADS}\n                preds = convert(ima, markers, n_timesteps, verbose=is_first_ecg)\n                \n                # Evaluate the signal-to-noise ratio, plot y_true vs. y_pred\n                if is_first_ecg:\n                    _, axs = plt.subplots(6, 2, figsize=(12, 18))\n                sum_signal = 0\n                sum_noise = 0\n                for i, lead in enumerate(LEADS):\n                    label = labels[lead]\n                    label = label[~ label.isna()]\n                    pred = preds[lead]\n            \n                    aligned_pred = align_signals(label, pred, int(row.fs * MAX_TIME_SHIFT))\n                    p_signal, p_noise = compute_power(label, aligned_pred)\n                    sum_signal += p_signal\n                    sum_noise += p_noise\n    \n                    if is_first_ecg:\n                        ax = axs.T.ravel()[i]\n                        ax.set_title(lead)\n                        ax.plot(label.values, label='y_true')\n                        ax.plot(pred, label='y_pred')\n                        ax.set_xlabel('timestep')\n                        ax.set_ylabel('mV')\n                        ax.legend()\n                if is_first_ecg:\n                    plt.tight_layout()\n                    plt.suptitle('y_true vs. y_pred', y=1.01)\n                    plt.show()\n                snr = compute_snr(sum_signal, sum_noise)\n                print(f\"{idx=:4d} id={row.id:10d} {img_type=:2d} SNR: {snr:5.2f}\")\n                snr_list.append(snr)\n                index_list.append([idx, img_type])\n    \n                if is_first_ecg:\n                    print('\\n')\n            else:\n                snr_list.append(1)\n                index_list.append([idx, img_type])\n        is_first_ecg = False\n    \n    snr = np.array(snr_list).mean()\n    val_score = max(float(10 * np.log10(snr)), -PERFECT_SCORE)\n    print(f\"# Average SNR: {snr:.2f} {val_score=:.2f}\")\n    snr_df = pd.DataFrame(index_list, columns=['idx', 'type'])\n    snr_df['snr'] = snr_list\n    snr_df.to_csv('~snr.csv', index=False)\n\nvalidate_algorithm(train.iloc[100:110], image_types=[3, 11], convert=convert_scanned_color)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-27T16:02:06.468209Z","iopub.execute_input":"2025-10-27T16:02:06.468983Z","iopub.status.idle":"2025-10-27T16:02:27.531614Z","shell.execute_reply.started":"2025-10-27T16:02:06.468956Z","shell.execute_reply":"2025-10-27T16:02:27.529784Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Convert the test images\n\nWe convert all test images of img_types 3 and 11 (the images which were scanned in color). For all other images (i.e., grayscale images and mobile photos), we submit the average training label.","metadata":{}},{"cell_type":"code","source":"def is_color_image(ima):\n    \"\"\" Test if a 3-channel image has colors.\"\"\"\n    return ima.std(axis=2).mean() != 0","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-27T16:02:31.282912Z","iopub.execute_input":"2025-10-27T16:02:31.283250Z","iopub.status.idle":"2025-10-27T16:02:31.289562Z","shell.execute_reply.started":"2025-10-27T16:02:31.283228Z","shell.execute_reply":"2025-10-27T16:02:31.287772Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission_data = []\nold_id = None\nleads = None\nfor idx, row in test.iterrows():\n    if row.id != old_id:\n        path = f\"/kaggle/input/physionet-ecg-image-digitization/test/{row.id}.png\"\n        # path = '/kaggle/input/physionet-ecg-image-digitization/train/1006427285/1006427285-0004.png'\n        # path = '/kaggle/input/physionet-ecg-image-digitization/train/1006427285/1006427285-0011.png'\n        ima = cv2.imread(path)\n        shape = ima.shape\n        good_shape = shape[0] == 1652 # scanned images have 1652 rows\n        \n        if good_shape and is_color_image(ima): # the image has red gridlines -> interpret the image\n            # Find the 17 line endpoints\n            markers = mf.find_markers(ima)\n\n            # Convert the image to 12 time series\n            n_timesteps = {lead: row.fs * 10 if lead == 'II' else row.fs * 10 // 4 for lead in LEADS}\n            preds = convert_scanned_color(ima, markers, n_timesteps, verbose=False)\n\n        else: # we cannot interpret the image -> predict the mean\n            preds = None\n            \n        old_id = row.id\n\n    if row.lead == 'II':\n        assert row.number_of_rows == row.fs * 10\n    else:\n        assert row.number_of_rows == row.fs * 10 // 4\n\n    if preds is not None:\n        pred = preds[row.lead]\n    else:\n        pred = mean_dict[row.lead].mean(axis=0)\n        pred = np.interp(np.linspace(0, 1, row.number_of_rows),\n                         np.linspace(0, 1, len(pred)),\n                         pred)\n    assert len(pred) == row.number_of_rows\n\n    for timestep in range(row.number_of_rows):\n        signal_id = f\"{row.id}_{timestep}_{row.lead}\"\n        submission_data.append({\n            'id': signal_id,\n            'value': pred[timestep]\n        })\n\nsubmission_df = pd.DataFrame(submission_data)\nprint(f\"Length: {len(submission_df)}\")\nsubmission_df.to_csv('submission.csv', index=False)\n!head submission.csv","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-27T16:02:38.395138Z","iopub.execute_input":"2025-10-27T16:02:38.395655Z","iopub.status.idle":"2025-10-27T16:02:40.651743Z","shell.execute_reply.started":"2025-10-27T16:02:38.395632Z","shell.execute_reply":"2025-10-27T16:02:40.649903Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}