{"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":"gpu","dataSources":[{"sourceId":97984,"databundleVersionId":14096757,"sourceType":"competition"}],"dockerImageVersionId":31153,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Explained baseline for the *PhysioNet - Digitization of ECG Images* competition\n\nThis notebook contributes three 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 with a top-down plane sweep (function `convert_scanned_color()`).\n1. It shows how to extract the time series from grayscale images with a neural network (function `convert_scanned_grayscale()`).\n\nYou can run this notebook with or without GPU, depending on your patience and your GPU quota.\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. → processed by `convert_scanned_grayscale()`\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. → processed by `convert_scanned_grayscale()`\n\nThe sampling frequencies in train are 250, 256, 500, 512, 1000, 1025 per second.\n\n## Hall of fame / acknowledgements\n\nThis notebook was influenced by the work of\n- @sasaleaf ([Edge Fix | small Improvement](https://www.kaggle.com/code/sasaleaf/edge-fix-small-improvement))\n- @guntasdhanjal ([ECG Digitization - Signal Smoothing Enhancement](https://www.kaggle.com/code/guntasdhanjal/ecg-digitization-signal-smoothing-enhancement))\n- @seowoohyeon ([PhysioNet_adjust](https://www.kaggle.com/code/seowoohyeon/physionet-adjust))\n- @antonoof ([large EDA and statistical model]())\n\nI thank them for publishing their ideas.","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\nfrom scipy.signal import medfilt\n\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import r2_score\n\nfrom tensorflow.keras.layers import Input, Dense, Activation, Reshape, GaussianNoise\nfrom tensorflow.keras.initializers import Constant\nfrom tensorflow.keras.models import Model\nfrom tensorflow.keras.optimizers import Adam\nfrom tensorflow.keras.losses import MeanSquaredError\nfrom tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau, TerminateOnNaN","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-02T16:36:22.705152Z","iopub.execute_input":"2025-11-02T16:36:22.705487Z","iopub.status.idle":"2025-11-02T16:36:44.599035Z","shell.execute_reply.started":"2025-11-02T16:36:22.705464Z","shell.execute_reply":"2025-11-02T16:36:44.597658Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Auxiliary functions","metadata":{}},{"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-11-02T16:36:44.600080Z","iopub.execute_input":"2025-11-02T16:36:44.600817Z","iopub.status.idle":"2025-11-02T16:36:44.624986Z","shell.execute_reply.started":"2025-11-02T16:36:44.600791Z","shell.execute_reply":"2025-11-02T16:36:44.623629Z"},"jupyter":{"source_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def plot_training_history(history):\n    \"\"\"Plot a Keras training history.\"\"\"\n    if len(history['loss']) >= 2:\n        _, axs = plt.subplots(1, 1, figsize=(6, 3), squeeze=False)\n        axs = axs.ravel()\n        axs[0].plot(np.arange(len(history['loss'])) + 1, history['loss'], ':', label='train_loss')\n        axs[0].plot(np.arange(len(history['val_loss'])) + 1, history['val_loss'], label='val_loss')\n        axs[0].legend()\n        axs[0].set_title('Training history')\n        plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-02T16:36:44.627583Z","iopub.execute_input":"2025-11-02T16:36:44.628002Z","iopub.status.idle":"2025-11-02T16:36:44.766552Z","shell.execute_reply.started":"2025-11-02T16:36:44.627969Z","shell.execute_reply":"2025-11-02T16:36:44.765328Z"},"_kg_hide-input":true,"jupyter":{"source_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Reading the metadata and the labels","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')\nlabel_dict = {}\nfor idx, row in tqdm(train.iterrows(), total=len(train)):\n    label_dict[idx] = pd.read_csv(f'/kaggle/input/physionet-ecg-image-digitization/train/{row.id}/{row.id}.csv')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-02T16:36:44.767640Z","iopub.execute_input":"2025-11-02T16:36:44.767946Z","iopub.status.idle":"2025-11-02T16:36:58.938682Z","shell.execute_reply.started":"2025-11-02T16:36:44.767918Z","shell.execute_reply":"2025-11-02T16:36:58.937445Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Looping through the training dataset\n\nWe define a generator function `train_images_and_labels`, which loops through a subset of the training images.","metadata":{}},{"cell_type":"code","source":"def train_images_and_labels(start=None, end=None, image_types=None, use_tqdm=True):\n    \"\"\"Generator function which yields a subset of the training images.\n    \n    Parameters\n    start: start index of slice\n    end: end index of slice\n    image_types: list of image types to select\n    \"\"\"\n    t = train.iloc[start:end]\n    iterable = t.iterrows()\n    if use_tqdm:\n        iterable = tqdm(iterable, total=len(t))\n    for idx, row in iterable:\n        png_paths = sorted(glob(f'/kaggle/input/physionet-ecg-image-digitization/train/{row.id}/{row.id}-*.png'))\n        labels = label_dict[idx]\n        for path in png_paths:\n            img_type = int(path[-8:-4])\n            if image_types is None or img_type in image_types:\n                ima = cv2.imread(path)\n\n                # The following lines document the possible shapes for every image type in train\n                # Test files may be different\n                shape = ima.shape\n                assert len(shape) == 3\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) <= (shape in {(3024, 4032, 3), (1344, 1008, 3), (4032, 3024, 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\n                yield idx, ima, img_type, labels\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-02T16:36:58.940018Z","iopub.execute_input":"2025-11-02T16:36:58.940350Z","iopub.status.idle":"2025-11-02T16:36:58.951822Z","shell.execute_reply.started":"2025-11-02T16:36:58.940320Z","shell.execute_reply":"2025-11-02T16:36:58.950049Z"}},"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. \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 = label_dict[idx]\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        # Upsample 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=(6, 1.5))\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.gca().get_xaxis().set_visible(False)\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 = label_dict[idx]\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        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\n# train_test_split_loc = 780\n# mean_dict = fit_mean_model(train.iloc[:train_test_split_loc], verbose=True)\n# validate_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=True)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-02T16:36:58.953335Z","iopub.execute_input":"2025-11-02T16:36:58.953690Z","iopub.status.idle":"2025-11-02T16:37:09.151617Z","shell.execute_reply.started":"2025-11-02T16:36:58.953666Z","shell.execute_reply":"2025-11-02T16:37:09.150460Z"},"_kg_hide-output":true,"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Finding lead endpoints with MarkerFinder\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 inferred 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        # Transform 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 = np.full((17, 2), -1)\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],\n                                    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],\n                                       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-11-02T17:33:37.694041Z","iopub.execute_input":"2025-11-02T17:33:37.694473Z","iopub.status.idle":"2025-11-02T17:33:39.174770Z","shell.execute_reply.started":"2025-11-02T17:33:37.694443Z","shell.execute_reply":"2025-11-02T17:33:39.173481Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Solution 1: A plane sweep for digitizing scanned color 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.\n\nAnd why 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","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    # Extend into columns without black pixels\n    median_top = int(np.median(top))\n    top[top == 0] = median_top\n    top[top > median_top + 250] = median_top\n    \n    # Bound from above\n    strip_width = 64\n    for strip_left in range(0, ima.shape[1], strip_width):\n        median_top_strip = int(np.median(top[strip_left:strip_left+strip_width]))\n        strip = ima[:median_top_strip-100, strip_left:strip_left+strip_width]\n        all_white = strip.all(axis=1)\n        if all_white.size > 0:\n            first_white_row = np.argmax(all_white[::-1])\n            if first_white_row > 0 or all_white[-1]:\n                first_white_row = median_top_strip - 100 - first_white_row\n                ami = np.argmin(ima[first_white_row:, strip_left:strip_left+strip_width], axis=0)\n                top[strip_left:strip_left+strip_width] = np.where(ami != 0, first_white_row + ami, median_top)\n    top[top > median_top + 250] = median_top\n\n    # Bound from below\n    strip_width = 64\n    for strip_left in range(0, ima.shape[1], strip_width):\n        median_top_strip = int(np.median(top[strip_left:strip_left+strip_width]))\n        strip = ima[median_top_strip+80:, strip_left:strip_left+strip_width]\n        all_white = strip.all(axis=1)\n        if all_white.size > 0:\n            first_white_row = np.argmax(all_white)\n            if first_white_row > 0 or all_white[0]:\n                first_white_row += median_top_strip + 80\n                mask = top > first_white_row\n                mask[:strip_left] = False\n                mask[strip_left+strip_width:] = False\n                top[mask] = median_top_strip\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, n_timesteps, 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    n_timesteps: 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    baseline = baseline[:len(pred0)] # in case end is outside the image (only idx 202 and some test)\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 pixels can be obscured by the marker\n        pred[:4] = np.where(pred[:4] > 0.2, pred[4], 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:] = np.where(pred[-5:] > 0.2, pred[-6], pred[-5:])\n    if lead in ['I', 'II-subset', 'III', 'II']:\n        # first two pixels can be obscured by the marker\n        pred[:2] = pred[2]\n\n    # Upsample\n    pred = np.interp(np.linspace(0, 1, n_timesteps),\n                     np.linspace(0, 1, len(pred)),\n                     pred)\n\n    # Fix implausible predictions\n    # From https://www.kaggle.com/code/antonoof/large-eda-and-statistical-model\n    OUTLIER_LOW_THRESHOLD = -1.5\n    OUTLIER_HIGH_THRESHOLD = 0.9\n    outlier_mask = (pred < OUTLIER_LOW_THRESHOLD) | (pred > OUTLIER_HIGH_THRESHOLD)\n    if np.any(outlier_mask):\n        for i in np.where(outlier_mask)[0]:\n            start_idx = max(0, i - 3)\n            end_idx = min(len(pred), i + 4)\n            neighbors = pred[start_idx:end_idx]\n            valid_neighbors = neighbors[(neighbors >= OUTLIER_LOW_THRESHOLD) & (neighbors <= OUTLIER_HIGH_THRESHOLD)]\n            if len(valid_neighbors) > 0:\n                pred[i] = np.median(valid_neighbors)\n            else:\n                pred[i] = 0\n\n    # Median filter proposed by @guntasdhanjal\n    pred = medfilt(pred, kernel_size=5)\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    # Clean the bottom border\n    baseline_II = markers[15:17][0].max() - crop_top\n    all_white = ima[baseline_II:, markers[15][1]:markers[16][1]].all(axis=1)\n    ama = np.argmax(all_white)\n    if ama > 0 or all_white[0]:\n        ima[baseline_II + ama:] = True\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.show()\n        # cv2.imwrite('ima-bw.png', ima.astype(int) * 255)\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    if verbose:\n        left = max(markers[:, 1].min() - 80, 0)\n        right = markers[:, 1].max() + 80\n        _, axs = plt.subplots(4, 1, sharex=True, figsize=(12, 8))\n        for i in range(4):\n            for j in range(5):\n                if 5*i+j < len(markers):\n                    axs[i].axvline(markers[5*i+j][1], color='gray')\n            axs[i].plot(np.arange(left, right), tops[i][left:right], color='b', label='top')\n            axs[i].plot(np.arange(left, right), bottoms[i][left:right], color='m', label='bottom')\n            axs[i].invert_yaxis()\n            axs[i].set_ylabel('pixel')\n            axs[i].legend()\n            axs[i].set_title(f\"Extracted line {i}\")\n        plt.show()\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    # 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    The three equalities are\n    I + III - II = 0\n    aVR + avL + aVF = 0\n    2 * aVR - 2 * aVF + 3 * II = 0\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\n    residual = 2 * preds['aVR'] - 2 * preds['aVF'] + 3 * preds['II'][len(preds['I']):len(preds['I'])+len(preds['aVR'])]\n    correction = residual / 17\n    preds['aVR'] -= 2 * correction\n    preds['aVF'] += 2 * correction\n    preds['II'][len(preds['I']):len(preds['I'])+len(preds['aVR'])] -= 3 * correction\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-02T16:37:10.803620Z","iopub.execute_input":"2025-11-02T16:37:10.803941Z","iopub.status.idle":"2025-11-02T16:37:10.840927Z","shell.execute_reply.started":"2025-11-02T16:37:10.803917Z","shell.execute_reply":"2025-11-02T16:37:10.839596Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Validation of solution 1\n\nWe digitize a few training images, plot the output and compute the signal-to-noise ratio. \n\nIf you look at the diagrams closely, you'll easily get ideas for improvements. The diagrams show:\n1. The original image\n2. The denoised image\n3. Four extracted ECG lines\n4. y_true vs. y_pred\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(convert, start=None, end=None, image_types=None):\n    \"\"\"Convert a few training images, plot the output and compute the signal-to-noise ratio\"\"\"\n    if len(test) != 24:\n        # test file when saving has 2*12=24 rows\n        return # when submitting, nobody will see the validation output\n    snr_list = []\n    index_list = []\n    is_first_ecg = True # we plot only the images of the first ECG\n    for idx, ima, img_type, labels in train_images_and_labels(start=start,\n                                                              end=end,\n                                                              image_types=image_types,\n                                                              use_tqdm=False):\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, 14))\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} {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            is_first_ecg = False\n\n    print(np.array(snr_list).mean(), len(snr_list))\n    snr = (np.array(snr_list).mean() - 1) / 9 * len(image_types) + 1\n    val_score = max(float(10 * np.log10(snr)), -PERFECT_SCORE)\n    print(f\"# Average SNR: {snr:.2f} {val_score=:.2f} {image_types}\")\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(convert_scanned_color, start=400, end=410, image_types=[3, 11])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-02T17:33:52.587606Z","iopub.execute_input":"2025-11-02T17:33:52.587962Z","iopub.status.idle":"2025-11-02T17:36:04.873073Z","shell.execute_reply.started":"2025-11-02T17:33:52.587940Z","shell.execute_reply":"2025-11-02T17:36:04.871785Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Solution 2: A neural network for digitizing grayscale images\n\nPlane sweep is a nice algorithm, but it fails when the image is too noisy. It can't even deal with black gridlines. So how can we transform this ECG image digitization competition into a tractable regression task?\n\nThanks to the `MarkerFinder` class, we know where the ECG lines begin and end, and we can even map the true labels to pixels of the image. We now slice the image between these markers into thin vertical stripes (600 pixels high, 11 pixels wide). The ECG lines intersect the stripes at a certain height, and predicting this height is a regression task.\n\nLook at the diagram below, which shows some random stripes. The regression targets (numbers between 0 and 600) are shown above every stripe. Try to discern the following items:\n1. ECG lines intersecting the stripe (the stripes are high enough so that there is often more than one line, one of them should be at the position indicated at the top of the diagram)\n2. horizontal gridlines\n3. vertical gridlines\n4. markers showing the start of a lead\n5. lead designations (letters 'a' and 'V')\n6. other noise\n\nThe size of the stripes (600\\*11 = 6600 pixels) has been chosen so that the neural network (with 6600 inputs and one output) can be trained easily in Kaggle notebooks.\n\n## Preparing the training data\n\nWe put the training data into arrays X of shape (40000, 6600) and y of shape (40000, ).","metadata":{}},{"cell_type":"code","source":"# Preparing the training data\nn_train_nn, n_training_samples_per_line = 40000, 30\nh0, w0 = 600, 11\n\nX = np.zeros((n_train_nn, h0 * w0), dtype=np.float32)\ny = np.zeros((n_train_nn, ), dtype=np.float32)\n\ndef prepare_dataset(train):\n    i_train_nn = 0\n    rng = np.random.default_rng(1)\n    t = tqdm(total=n_train_nn, position=0)\n    for idx, row in train.iterrows():\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            if img_type in [4, 12]:\n                ima = cv2.imread(path)\n                markers = mf.find_markers(ima)\n                ima = ima.mean(axis=2)\n                ll = np.array([labels['aVR'].fillna(0) + labels['V1'].fillna(0),\n                               labels['aVL'].fillna(0) + labels['V2'].fillna(0),\n                               labels['aVF'].fillna(0) + labels['V3'].fillna(0),])\n                for i in range(3):\n                    for j in range(1, 3):\n                        yl, xl = markers[5*i+j]\n                        yp, xp = markers[5*i+j+1] - markers[5*i+j]\n                        lead = ['', 'aVR', 'V1', '', '', '', 'aVL', 'V2', '', '', '' ,'aVF', 'V3'][5*i+j]\n                        l = labels[lead]\n                        l = l[~l.isna()]\n                        l = np.interp(np.linspace(0, 1, xp),\n                                      np.linspace(0, 1, len(l)),\n                                      l) * 80 # true labels in pixel coordinates\n                        assert xp > 0 # yp may be negative\n                        for k in range(n_training_samples_per_line):\n                            alpha = rng.uniform()\n                            x0, y0 = int(xl + alpha * xp), int(yl + alpha * yp)\n                            X[i_train_nn] = ima[y0-h0//2:y0+h0//2, x0-w0//2:x0+w0//2+1].ravel()\n                            y[i_train_nn] = h0 / 2 - l[x0-xl]\n                            i_train_nn += 1\n                            t.update(1)\n                            if i_train_nn == n_train_nn:\n                                t.close()\n                                print(f\"last idx used: {idx}\")\n                                return\n        \nprepare_dataset(train.iloc[200:])\nprint(X.shape, y.shape, X.size*4) # (n_train_nn, w0*h0) (n_train_nn, ) 40000 -> last idx used: 311\nX_orig = X.copy()\nX = X / 255.0 - 1 # offset so that X.max() = 0\nassert np.isfinite(X).all()\nassert np.isfinite(y).all()\n\n# X is a scaled pixel intensity, 0 = white, -1 = black\n# y is amplitude in pixels, usually between 0 and h0-1, h0//2 is baseline. In rare cases it is below 0 or above h0","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-02T16:39:31.481545Z","iopub.execute_input":"2025-11-02T16:39:31.481990Z","iopub.status.idle":"2025-11-02T16:41:24.965777Z","shell.execute_reply.started":"2025-11-02T16:39:31.481967Z","shell.execute_reply":"2025-11-02T16:41:24.964736Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"l = [17254, 23409, 3713, 6368, 13561, 12331, 21, 9871, 3718, 21152, 9611, 5918, 22371, 33636, 22909, 14500]\n_, axs = plt.subplots(1, len(l), sharey=True, sharex=True, figsize=(9, 12))\nfor i, j in enumerate(l):\n    # print(i, j)\n    axs[i].imshow(X[j].reshape((h0, w0)))\n    axs[i].set_title(f\"{y[j]:.0f}\")\n    axs[i].get_xaxis().set_visible(False)\n# plt.savefig('stripes.png')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-02T16:41:24.966965Z","iopub.execute_input":"2025-11-02T16:41:24.967328Z","iopub.status.idle":"2025-11-02T16:41:26.357796Z","shell.execute_reply.started":"2025-11-02T16:41:24.967295Z","shell.execute_reply":"2025-11-02T16:41:26.356288Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## The network\n\nWe implement a neural network with four hidden layers (7 million parameters) and train it. If everything goes well, the R2 score should be about 0.75, corresponding to a signal-to-noise ratio of 4. Feel free to experiment with other network architectures.","metadata":{}},{"cell_type":"code","source":"def regression_model():\n    \"\"\"A dense feed-forward with 4 hidden layers.\"\"\"\n    initial_lr = 0.005\n    noise = 0.3\n    x_input = Input(shape=(X.shape[-1], ))\n    x = x_input\n    x = GaussianNoise(noise)(x)\n    x = Dense(units=1024, activation='selu')(x)\n    x = Dense(units=512, activation='selu')(x)\n    x = Dense(units=256, activation='selu')(x)\n    x = Dense(units=256, activation='selu')(x)\n    x_output = Dense(units=1, activation='linear', bias_initializer=Constant(h0 / 2))(x)\n\n    model = Model(inputs=x_input, outputs=x_output)\n    return model, initial_lr\n\nmodel, _ = regression_model()\nmodel(X[0:5]).shape\nmodel.summary()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-02T16:41:26.358992Z","iopub.execute_input":"2025-11-02T16:41:26.359286Z","iopub.status.idle":"2025-11-02T16:41:26.674194Z","shell.execute_reply.started":"2025-11-02T16:41:26.359263Z","shell.execute_reply":"2025-11-02T16:41:26.673212Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Train the regression model\n# Inputs for training: X and y\nepochs = 70\nverbose = 2\nbatch_size = 64\n\ngrayscale_model_list = []\noof = np.zeros_like(y)\nkf = KFold(shuffle=True, random_state=1)\nfor fold, (idx_tr, idx_va) in enumerate(kf.split(X)):\n    X_tr = X[idx_tr]\n    X_va = X[idx_va]\n    y_tr = y[idx_tr]\n    y_va = y[idx_va]\n\n    model, initial_lr = regression_model()\n\n    # Train the complete model\n    model.compile(\n        optimizer=Adam(learning_rate=initial_lr),\n        loss=MeanSquaredError(),\n    )\n    print(y_va.var())\n    history = model.fit(\n        X_tr, y_tr,\n        validation_data=(X_va, y_va),\n        epochs=epochs,\n        batch_size=batch_size,\n        callbacks=[EarlyStopping(patience=6, min_delta=0.1),\n                   ReduceLROnPlateau(factor=0.5, patience=2, verbose=1, min_delta=0.1, min_lr=initial_lr/63),\n                   TerminateOnNaN()],\n        verbose=verbose\n    )\n    history = history.history\n    plot_training_history(history)\n\n    y_pred = model.predict(X_va, batch_size=1024, verbose=0).ravel()\n    oof[idx_va] = y_pred\n    mse = np.square(y_va - y_pred).mean()\n    r2 = r2_score(y_va, y_pred)\n    print(f'# Fold {fold} {mse:6.2f} {r2:4.2f} {np.var(y_va):6.2f} {len(X)}*{h0}*{w0} {batch_size}\\n')\n    grayscale_model_list.append(model)\n\n    # Regression\n    plt.figure(figsize=(16, 4))\n    plt.subplot(1, 2, 1)\n    plt.hist(y, bins=50, density=True, label='y_true')\n    plt.hist(y_pred, bins=50, density=True, alpha=0.6, label='y_pred')\n    plt.legend()\n    plt.subplot(1, 2, 2)\n    plt.scatter(y_pred, y_va, s=1)\n    plt.xlabel('y_pred')\n    plt.ylabel('y_true')\n    plt.gca().set_aspect('equal')\n    plt.show()\n    break\n\n# Fold 0  75.41 0.76 317.96 *Model 0 initial_lr=0.005 noise=0.3* 40000*600*11 64","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-02T16:41:26.675207Z","iopub.execute_input":"2025-11-02T16:41:26.675582Z","iopub.status.idle":"2025-11-02T16:54:55.986762Z","shell.execute_reply.started":"2025-11-02T16:41:26.675553Z","shell.execute_reply":"2025-11-02T16:54:55.984707Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"del X, y","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-02T16:54:55.989158Z","iopub.execute_input":"2025-11-02T16:54:55.989533Z","iopub.status.idle":"2025-11-02T16:54:56.002235Z","shell.execute_reply.started":"2025-11-02T16:54:55.989510Z","shell.execute_reply":"2025-11-02T16:54:56.000948Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Validation of solution 2\n\nWe've used train.iloc\\[200:\\] to train the neural network. Consequently, we may use train.iloc\\[:200\\] for validating. We process a few images and evaluate the signal-to-noise ratio. \n\nBy the way, the grayscale model works for color images, too, so that we can ensemble the neural network predictions with the plane sweep predictions.","metadata":{}},{"cell_type":"code","source":"def get_grayscale_lead(ima, lead, n_timesteps, markers):\n    \"\"\"Extract and resample one lead from an ECG image.\n\n    The function returns the time series if possible, or it gives up and returns None.\n    \n    Parameters:\n    ima: single-channel (grayscale image)\n    lead: one of the 12 lead labels (string)\n    n_timesteps: number of samples required (int)\n    markers: 17 markers as list of size-2 integer arrays (row, column)\n\n    Returns:\n    pred: either an array of shape (n_timesteps, ) or None\n    \"\"\"\n    line, begin, end = mf.lead_info(lead)\n\n    X_list = []\n    width = markers[end][1] - markers[begin][1]\n    slope = (markers[end][0] - markers[begin][0]) / width\n\n    for x0 in range(markers[begin][1], markers[end][1]):\n        y0 = int(markers[begin][0] + slope * (x0 - markers[begin][1]))\n        X = ima[y0-h0//2:y0+h0//2, x0-w0//2:x0+w0//2+1]\n        if X.size == 0:\n            X = np.full((h0, w0), 255) # all white\n        else:\n            if X.shape[0] < h0:\n                X = np.vstack([X, np.full((h0 - X.shape[0], X.shape[1]), 255)]) # pad bottom white\n            if X.shape[1] < w0:\n                X = np.hstack([X, np.full((X.shape[0], w0 - X.shape[1]), 255)]) # pad right white\n        X_list.append(X)\n    X = np.stack(X_list).reshape(len(X_list), -1).astype(np.float32)\n    assert X.shape[1] == w0 * h0\n    assert np.isfinite(X).all()\n\n    X = X / 255.0 - 1\n\n    pred = h0 / 2 - np.mean([model.predict(X, batch_size=1024, verbose=0).ravel() for model in grayscale_model_list], axis=0)\n\n    # Scale\n    pred /= 80 # 80 pixels = 1 mV\n\n    # Upsample\n    pred = np.interp(np.linspace(0, 1, n_timesteps),\n                     np.linspace(0, 1, len(pred)),\n                     pred)\n    \n    # Fix implausible predictions\n    pred = np.where(np.abs(pred) <= 2, pred, 0)\n        \n    return pred\n\n\ndef convert_scanned_grayscale(ima, markers, n_timesteps, verbose=False):\n    \"\"\"Convert a scanned grayscale image (type 4 or 12) to 12 time series.\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 or None\n    \"\"\"\n    # Drop the color channels / convert to grayscale\n    ima = ima.mean(axis=2)\n\n    n_timesteps['II-subset'] = n_timesteps['I']\n    preds = {}\n    for i, lead in enumerate(LEADS + ['II-subset']):\n        preds[lead] = get_grayscale_lead(ima, lead, n_timesteps[lead], markers)\n\n    if preds['II'] is not None and preds['II-subset'] is not None:\n        preds['II'][:len(preds['II-subset'])] = (preds['II'][:len(preds['II-subset'])] + preds['II-subset']) / 2\n    elif preds['II-subset'] is not None:\n        preds['II'] = np.zeros(n_timesteps['II'])\n        preds['II'][:len(preds['II-subset'])] = preds['II-subset']\n    del preds['II-subset']\n\n    # Apply Einthoven's law\n    apply_einthoven(preds)\n    \n    return preds\n\ndef convert_scanned_both(ima, markers, n_timesteps, verbose=False):\n    \"\"\"Convert a scanned color image (type 3 or 11) to 12 time series.\n\n    This function ensembles the neural network with the plane sweep algorithm.\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 or None\n    \"\"\"\n    keys = list(n_timesteps.keys())\n    pred1 = convert_scanned_grayscale(ima, markers, n_timesteps, verbose=False)\n    pred2 = convert_scanned_color(ima, markers, n_timesteps, verbose=False)\n\n    def ensemble_two(pred1, pred2):\n        if pred1 is not None and pred2 is not None:\n            return (pred1 + pred2) / 2\n        if pred1 is not None:\n            return pred1\n        return pred2\n        \n    preds = {k: ensemble_two(pred1[k], pred2[k]) for k in keys}\n    return preds\n\n\nvalidate_algorithm(convert_scanned_grayscale, start=400, end=410, image_types=[4, 12])\nvalidate_algorithm(convert_scanned_both, start=400, end=410, image_types=[3, 11])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-02T17:37:03.814748Z","iopub.execute_input":"2025-11-02T17:37:03.815146Z","iopub.status.idle":"2025-11-02T18:01:28.597925Z","shell.execute_reply.started":"2025-11-02T17:37:03.815123Z","shell.execute_reply":"2025-11-02T18:01:28.596265Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Digitizing and submitting the test images\n\nWe digitize all test images of img_types 3, 4, 11 and 12 (the images which were acquired with a scanner). For all other images (i.e., mobile photos), we submit the average training labels.","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-11-02T17:18:41.818853Z","iopub.execute_input":"2025-11-02T17:18:41.819301Z","iopub.status.idle":"2025-11-02T17:18:41.825450Z","shell.execute_reply.started":"2025-11-02T17:18:41.819259Z","shell.execute_reply":"2025-11-02T17:18:41.824093Z"}},"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:\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            if is_color_image(ima):\n                preds = convert_scanned_both(ima, markers, n_timesteps, verbose=False)\n            else:\n                preds = convert_scanned_grayscale(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-11-02T17:18:41.827010Z","iopub.execute_input":"2025-11-02T17:18:41.827315Z","iopub.status.idle":"2025-11-02T17:18:46.288108Z","shell.execute_reply.started":"2025-11-02T17:18:41.827294Z","shell.execute_reply":"2025-11-02T17:18:46.286908Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}