diff --git a/src/maldi_tools/extraction.py b/src/maldi_tools/extraction.py
index cab87dd..cb740a0 100644
--- a/src/maldi_tools/extraction.py
+++ b/src/maldi_tools/extraction.py
@@ -16,7 +16,7 @@
 import pandas as pd
 import xarray as xr
 from pyimzml.ImzMLParser import ImzMLParser
-from scipy import signal
+from scipy import integrate, signal
 from tqdm.notebook import tqdm
 
 from maldi_tools import plotting
@@ -232,13 +232,23 @@ def peak_spectra(
     return panel_df
 
 
-def coordinate_integration(peak_df: pd.DataFrame, imz_data: ImzMLParser) -> xr.DataArray:
+def coordinate_integration(
+    peak_df: pd.DataFrame,
+    l_ips_r: np.ndarray,
+    r_ips_r: np.ndarray,
+    peak_widths_height: np.ndarray,
+    imz_data: ImzMLParser,
+) -> xr.DataArray:
     """Integrates the coordinates with the discovered, post-processed peaks and generates an image for
     each of the peaks using the imzML coordinate data.
 
     Args:
     ----
         peak_df (pd.DataFrame): The unique peaks from the data.
+        l_ips_r (np.ndarray): The rounded left (lower) bound.
+        r_ips_r (np.ndarray): The rounded right (upper) bound.
+        peak_widths_height (np.ndarray): The height of the contour lines at which the peak widths were
+        calculated from.
         imz_data (ImzMLParser): The imzML object.
 
     Returns:
@@ -260,10 +270,12 @@ def coordinate_integration(peak_df: pd.DataFrame, imz_data: ImzMLParser) -> xr.D
     for idx, (x, y, _) in tqdm(enumerate(imz_data.coordinates), total=len(imz_data.coordinates)):
         mzs, intensities = imz_data.getspectrum(idx)
 
-        intensity: np.ndarray = intensities[np.isin(mzs, peak_df["m/z"])]
-
         for i_idx, peak in peak_df.loc[peak_df["m/z"].isin(mzs), "peak"].reset_index(drop=True).items():
-            imgs[peak_dict[peak], x - 1, y - 1] += intensity[i_idx]
+            left_idx = abs(intensities.values - l_ips_r[i_idx]).idxmin()
+            right_idx = abs(intensities.values - r_ips_r[i_idx]).idxmin()
+            imgs[peak_dict[peak], x - 1, y - 1] += integrate.simpson(
+                intensities.values[left_idx:right_idx]
+            ) - (peak_widths_height * (right_idx - left_idx))
 
     img_data = xr.DataArray(
         data=imgs,
diff --git a/templates/maldi-pipeline.ipynb b/templates/maldi-pipeline.ipynb
index 3bffb46..e342a77 100644
--- a/templates/maldi-pipeline.ipynb
+++ b/templates/maldi-pipeline.ipynb
@@ -466,7 +466,13 @@
    },
    "outputs": [],
    "source": [
-    "image_data = extraction.coordinate_integration(peak_df=peak_df, imz_data=imz_data)"
+    "image_data = extraction.coordinate_integration(\n",
+    "    peak_df=peak_df,\n",
+    "    l_ips_r=l_ips_r,\n",
+    "    r_ips_r=r_ips_r,\n",
+    "    peak_widths_height=peak_widths_height,\n",
+    "    imz_data=imz_data,\n",
+    ")"
    ]
   },
   {
@@ -598,7 +604,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.1"
+   "version": "3.11.6"
   },
   "vscode": {
    "interpreter": {
diff --git a/tests/extraction_test.py b/tests/extraction_test.py
index c1571cd..7c51d9d 100644
--- a/tests/extraction_test.py
+++ b/tests/extraction_test.py
@@ -97,7 +97,16 @@ def test_peak_spectra(
 
 def test_coordinate_integration(imz_data, peak_widths):
     peak_df, *_ = peak_widths
-    img_data = extraction.coordinate_integration(peak_df=peak_df, imz_data=imz_data)
+    l_ips_r = np.arange(0, peak_df.shape[0], 5)
+    r_ips_r = l_ips_r + 2
+    peak_widths_height = np.repeat(2, peak_df.shape[0])
+    img_data = extraction.coordinate_integration(
+        peak_df=peak_df,
+        l_ips_r=l_ips_r,
+        r_ips_r=r_ips_r,
+        peak_widths_height=peak_widths_height,
+        imz_data=imz_data,
+    )
 
     # Make sure the shape of any given image is correct.
     assert img_data.shape[1:] == (10, 10)