diff --git a/src/load_distribution/load_distribution.py b/src/load_distribution/load_distribution.py
index f4877dd..4bfab61 100644
--- a/src/load_distribution/load_distribution.py
+++ b/src/load_distribution/load_distribution.py
@@ -215,14 +215,18 @@ def singularities_to_polygon(los: list[Singularity], xy: bool = False) -> Polygo
     # required shape, even if that means the exact x value is omitted (because we are
     # keeping the value immediately to the left and immediately to the right instead).
     x_acc = sorted(list(set(x_acc)))
+
     x_ord_count = Counter([round(x, 6) for x in x_acc])
     to_filter = []
     for key, count in x_ord_count.items():
         if count == 3:
             to_filter.append(key)
-    for filter_val in to_filter:
-        index = x_acc.index(filter_val)
-        x_acc.pop(index)
+
+    if to_filter:
+        rounded_x_acc = [round(x, 6) for x in x_acc]
+        for filter_val in to_filter:
+            index = rounded_x_acc.index(filter_val)
+            x_acc.pop(index)
 
     # Now, for every x value, compute the corresponding y value
     y_acc = [sum([sing(x) for sing in sorted_sings]) for x in x_acc[:-1]]
diff --git a/tests/test_load_distribution.py b/tests/test_load_distribution.py
index 4417c05..610ec39 100644
--- a/tests/test_load_distribution.py
+++ b/tests/test_load_distribution.py
@@ -173,6 +173,25 @@ def test_singularities_to_polygon():
     square_bow_tie_sings = ld.get_singularity_functions(square_bow_tie)
     square_with_L_hole_sings = ld.get_singularity_functions(square_with_L_hole)
 
+    real_world_singularities = [
+        ld.Singularity(
+            x0=0.034367933446777243,
+            x1=0.051573454339195024,
+            m=29.060439560440315,
+            y0=0.0,
+            precision=6,
+            eps=1e-12,
+        ),
+        ld.Singularity(
+            x0=0.051573454339195024,
+            x1=11.704999999999998,
+            m=0.0,
+            y0=0.5,
+            precision=6,
+            eps=1e-12,
+        ),
+    ]
+
     assert (
         ld.singularities_to_polygon(square_45_sings[0] + square_45_sings[1]).wkt
         == "POLYGON ((0 0, 0 0, 5 10, 5 10, 10 0, 10 0, 0 0))"
@@ -202,13 +221,22 @@ def test_singularities_to_polygon():
         == "POLYGON ((1 0, 1 10, 3 10, 3 0, 1 0))"
     )
     assert (
-        ld.singularities_to_polygon([
-            ld.Singularity(x0=0, x1=4, m=0.0, y0=5, precision=6, eps=1e-12),
-            ld.Singularity(x0=0, x1=2.5, m=0.0, y0=40, precision=6, eps=1e-12)
-        ]).wkt
+        ld.singularities_to_polygon(
+            [
+                ld.Singularity(x0=0, x1=4, m=0.0, y0=5, precision=6, eps=1e-12),
+                ld.Singularity(x0=0, x1=2.5, m=0.0, y0=40, precision=6, eps=1e-12),
+            ]
+        ).wkt
         == "POLYGON ((0 0, 0 45, 2.5 45, 2.5 5, 4 5, 4 0, 0 0))"
     )
 
+    # This test fails due to an error in filtering
+    assert (
+        ld.singularities_to_polygon(real_world_singularities).wkt
+        == "POLYGON ((0 0, 0.034368 0, 0.051573 0.5, 0.051573 0.5, 11.705 0.5, 11.705 0, 0 0))"
+    )
+
+
 def test_overlap_region_to_singularity():
     assert ld.overlap_region_to_singularity(
         ld.Overlap(6.0, 10.0, 0.0, 0.0, 0.0, 10.0), 6