diff --git a/src/pe_coupling/geometry/PeIntersectionRatio.cpp b/src/pe_coupling/geometry/PeIntersectionRatio.cpp
index 8d9fde471086ae01aed4fd554f44d1feca496031..dc4f5e7608ba33fe3fe84848144d3c27a92f0742 100644
--- a/src/pe_coupling/geometry/PeIntersectionRatio.cpp
+++ b/src/pe_coupling/geometry/PeIntersectionRatio.cpp
@@ -22,8 +22,11 @@
 
 #include "pe_coupling/geometry/PeIntersectionRatio.h"
 
+#include "core/math/Limits.h"
 #include "core/math/Matrix3.h"
 
+#include <algorithm>
+
 namespace walberla {
 namespace lbm {
 
@@ -48,9 +51,9 @@ real_t intersectionRatioSpherePe( const pe::Sphere & sphere,
    delta /= dirLength;
 
    WALBERLA_ASSERT_GREATER_EQUAL( delta, real_t( 0 ) );
-   WALBERLA_ASSERT_LESS_EQUAL( delta, real_t( 1 ) );
+   WALBERLA_ASSERT_LESS_EQUAL( delta - real_t(1), math::Limits<real_t>::accuracy(), delta );
 
-   return delta;
+   return std::min(std::max(delta, real_c(0)), real_c(1));
 
 }
 
@@ -74,9 +77,9 @@ real_t intersectionRatioPlanePe( const pe::Plane & plane,
    real_t delta = diff * planeNormal / denom;
 
    WALBERLA_ASSERT_GREATER_EQUAL( delta, real_t( 0 ) );
-   WALBERLA_ASSERT_LESS_EQUAL( delta, real_t( 1 ) );
+   WALBERLA_ASSERT_LESS_EQUAL( delta - real_t(1), math::Limits<real_t>::accuracy(), delta );
 
-   return delta;
+   return std::min(std::max(delta, real_c(0)), real_c(1));
 
 }
 
@@ -110,9 +113,9 @@ real_t intersectionRatioEllipsoidPe( const pe::Ellipsoid & ellipsoid,
    real_t delta = (-b - root) / (real_t(2) * a);
 
    WALBERLA_ASSERT_GREATER_EQUAL( delta, real_t( 0 ) );
-   WALBERLA_ASSERT_LESS_EQUAL( delta, real_t( 1 ) );
+   WALBERLA_ASSERT_LESS_EQUAL( delta - real_t(1), math::Limits<real_t>::accuracy(), delta );
 
-   return delta;
+   return std::min(std::max(delta, real_c(0)), real_c(1));
 
 }