diff --git a/.vscode/launch.json b/.vscode/launch.json
new file mode 100644
index 0000000..497c1b8
--- /dev/null
+++ b/.vscode/launch.json
@@ -0,0 +1,217 @@
+{
+    // Use IntelliSense to learn about possible attributes.
+    // Hover to view descriptions of existing attributes.
+    // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
+    "version": "0.2.0",
+    "configurations": [
+        {
+            "type": "lldb",
+            "request": "launch",
+            "name": "Debug unit tests in library 'collision'",
+            "cargo": {
+                "args": [
+                    "test",
+                    "--no-run",
+                    "--lib",
+                    "--package=collision"
+                ],
+                "filter": {
+                    "name": "collision",
+                    "kind": "lib"
+                }
+            },
+            "args": [],
+            "cwd": "${workspaceFolder}"
+        },
+        {
+            "type": "lldb",
+            "request": "launch",
+            "name": "Debug integration test 'serde'",
+            "cargo": {
+                "args": [
+                    "test",
+                    "--no-run",
+                    "--test=serde",
+                    "--package=collision"
+                ],
+                "filter": {
+                    "name": "serde",
+                    "kind": "test"
+                }
+            },
+            "args": [],
+            "cwd": "${workspaceFolder}"
+        },
+        {
+            "type": "lldb",
+            "request": "launch",
+            "name": "Debug integration test 'point'",
+            "cargo": {
+                "args": [
+                    "test",
+                    "--no-run",
+                    "--test=point",
+                    "--package=collision"
+                ],
+                "filter": {
+                    "name": "point",
+                    "kind": "test"
+                }
+            },
+            "args": [],
+            "cwd": "${workspaceFolder}"
+        },
+        {
+            "type": "lldb",
+            "request": "launch",
+            "name": "Debug integration test 'bound'",
+            "cargo": {
+                "args": [
+                    "test",
+                    "--no-run",
+                    "--test=bound",
+                    "--package=collision"
+                ],
+                "filter": {
+                    "name": "bound",
+                    "kind": "test"
+                }
+            },
+            "args": [],
+            "cwd": "${workspaceFolder}"
+        },
+        {
+            "type": "lldb",
+            "request": "launch",
+            "name": "Debug integration test 'dbvt'",
+            "cargo": {
+                "args": [
+                    "test",
+                    "--no-run",
+                    "--test=dbvt",
+                    "--package=collision"
+                ],
+                "filter": {
+                    "name": "dbvt",
+                    "kind": "test"
+                }
+            },
+            "args": [],
+            "cwd": "${workspaceFolder}"
+        },
+        {
+            "type": "lldb",
+            "request": "launch",
+            "name": "Debug integration test 'plane'",
+            "cargo": {
+                "args": [
+                    "test",
+                    "--no-run",
+                    "--test=plane",
+                    "--package=collision"
+                ],
+                "filter": {
+                    "name": "plane",
+                    "kind": "test"
+                }
+            },
+            "args": [],
+            "cwd": "${workspaceFolder}"
+        },
+        {
+            "type": "lldb",
+            "request": "launch",
+            "name": "Debug integration test 'frustum'",
+            "cargo": {
+                "args": [
+                    "test",
+                    "--no-run",
+                    "--test=frustum",
+                    "--package=collision"
+                ],
+                "filter": {
+                    "name": "frustum",
+                    "kind": "test"
+                }
+            },
+            "args": [],
+            "cwd": "${workspaceFolder}"
+        },
+        {
+            "type": "lldb",
+            "request": "launch",
+            "name": "Debug integration test 'sphere'",
+            "cargo": {
+                "args": [
+                    "test",
+                    "--no-run",
+                    "--test=sphere",
+                    "--package=collision"
+                ],
+                "filter": {
+                    "name": "sphere",
+                    "kind": "test"
+                }
+            },
+            "args": [],
+            "cwd": "${workspaceFolder}"
+        },
+        {
+            "type": "lldb",
+            "request": "launch",
+            "name": "Debug integration test 'aabb'",
+            "cargo": {
+                "args": [
+                    "test",
+                    "--no-run",
+                    "--test=aabb",
+                    "--package=collision"
+                ],
+                "filter": {
+                    "name": "aabb",
+                    "kind": "test"
+                }
+            },
+            "args": [],
+            "cwd": "${workspaceFolder}"
+        },
+        {
+            "type": "lldb",
+            "request": "launch",
+            "name": "Debug benchmark 'dbvt'",
+            "cargo": {
+                "args": [
+                    "test",
+                    "--no-run",
+                    "--bench=dbvt",
+                    "--package=collision"
+                ],
+                "filter": {
+                    "name": "dbvt",
+                    "kind": "bench"
+                }
+            },
+            "args": [],
+            "cwd": "${workspaceFolder}"
+        },
+        {
+            "type": "lldb",
+            "request": "launch",
+            "name": "Debug benchmark 'polyhedron'",
+            "cargo": {
+                "args": [
+                    "test",
+                    "--no-run",
+                    "--bench=polyhedron",
+                    "--package=collision"
+                ],
+                "filter": {
+                    "name": "polyhedron",
+                    "kind": "bench"
+                }
+            },
+            "args": [],
+            "cwd": "${workspaceFolder}"
+        }
+    ]
+}
\ No newline at end of file
diff --git a/Cargo.toml b/Cargo.toml
index 9ca2d0d..9724146 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,7 +1,7 @@
 [package]
 
 name = "collision"
-version = "0.20.1"
+version = "0.20.2"
 authors = ["Brendan Zabarauskas <bjzaba@yahoo.com.au>",
         "Brian Heylin",
         "Colin Sherratt",
diff --git a/benches/dbvt.rs b/benches/dbvt.rs
index 3697aaa..b9dabee 100644
--- a/benches/dbvt.rs
+++ b/benches/dbvt.rs
@@ -5,11 +5,11 @@ extern crate collision;
 extern crate rand;
 extern crate test;
 
-use cgmath::{Point2, Vector2};
 use cgmath::prelude::*;
-use collision::{Aabb2, Ray2};
+use cgmath::{Point2, Vector2};
 use collision::dbvt::*;
 use collision::prelude::*;
+use collision::{Aabb2, Ray2};
 use rand::Rng;
 use test::Bencher;
 
@@ -100,7 +100,8 @@ fn benchmark_query(b: &mut Bencher) {
         })
         .collect();
 
-    let mut visitors: Vec<DiscreteVisitor<Ray2<f32>, Value>> = rays.iter()
+    let mut visitors: Vec<DiscreteVisitor<Ray2<f32>, Value>> = rays
+        .iter()
         .map(|ray| DiscreteVisitor::<Ray2<f32>, Value>::new(ray))
         .collect();
 
diff --git a/benches/polyhedron.rs b/benches/polyhedron.rs
index e0f5d86..dde54bb 100644
--- a/benches/polyhedron.rs
+++ b/benches/polyhedron.rs
@@ -6,12 +6,12 @@ extern crate genmesh;
 extern crate rand;
 extern crate test;
 
-use cgmath::{Decomposed, Point3, Quaternion, Vector3};
 use cgmath::prelude::*;
+use cgmath::{Decomposed, Point3, Quaternion, Vector3};
 use collision::prelude::*;
 use collision::primitive::ConvexPolyhedron;
-use genmesh::Triangulate;
 use genmesh::generators::{IndexedPolygon, SharedVertex, SphereUV};
+use genmesh::Triangulate;
 use rand::Rng;
 use test::{black_box, Bencher};
 
@@ -93,12 +93,14 @@ fn dirs(n: usize) -> Vec<Vector3<f32>> {
 fn sphere(n: usize, with_faces: bool) -> ConvexPolyhedron<f32> {
     let gen = SphereUV::new(n, n);
 
-    let vertices = gen.shared_vertex_iter()
+    let vertices = gen
+        .shared_vertex_iter()
         .map(|v| Point3::from(v.pos))
         .collect::<Vec<_>>();
 
     if with_faces {
-        let faces = gen.indexed_polygon_iter()
+        let faces = gen
+            .indexed_polygon_iter()
             .triangulate()
             .map(|f| (f.x, f.y, f.z))
             .collect::<Vec<_>>();
diff --git a/src/algorithm/broad_phase/mod.rs b/src/algorithm/broad_phase/mod.rs
index 191f5b6..6f2081d 100644
--- a/src/algorithm/broad_phase/mod.rs
+++ b/src/algorithm/broad_phase/mod.rs
@@ -5,5 +5,5 @@ pub use self::dbvt::DbvtBroadPhase;
 pub use self::sweep_prune::{SweepAndPrune, SweepAndPrune2, SweepAndPrune3, Variance};
 
 mod brute_force;
-mod sweep_prune;
 mod dbvt;
+mod sweep_prune;
diff --git a/src/algorithm/broad_phase/sweep_prune.rs b/src/algorithm/broad_phase/sweep_prune.rs
index dddd342..1438c2a 100644
--- a/src/algorithm/broad_phase/sweep_prune.rs
+++ b/src/algorithm/broad_phase/sweep_prune.rs
@@ -126,7 +126,8 @@ where
         }
 
         // compute sweep axis for the next iteration
-        let (axis, _) = self.variance
+        let (axis, _) = self
+            .variance
             .compute_axis(NumCast::from(shapes.len()).unwrap());
         self.sweep_axis = axis;
 
@@ -138,8 +139,8 @@ mod variance {
     use std::marker;
 
     use crate::Bound;
-    use cgmath::{BaseFloat, Point2, Point3, Vector2, Vector3};
     use cgmath::prelude::*;
+    use cgmath::{BaseFloat, Point2, Point3, Vector2, Vector3};
 
     /// Trait for variance calculation in sweep and prune algorithm
     pub trait Variance {
diff --git a/src/algorithm/minkowski/epa/epa2d.rs b/src/algorithm/minkowski/epa/epa2d.rs
index 3500fec..18be46d 100644
--- a/src/algorithm/minkowski/epa/epa2d.rs
+++ b/src/algorithm/minkowski/epa/epa2d.rs
@@ -1,44 +1,45 @@
 use std::marker;
 
-use cgmath::{BaseFloat, Point2, Vector2};
-use cgmath::prelude::*;
-use cgmath::num_traits::NumCast;
 use approx::assert_ulps_ne;
+use approx::ulps_ne;
+use cgmath::num_traits::NumCast;
+use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point2, Vector2};
 
 use super::*;
-use crate::{CollisionStrategy, Contact};
 use crate::prelude::*;
-use crate::primitive::util::triple_product;
+use crate::primitive::util::vector_normal;
+use crate::{CollisionStrategy, Contact};
 
 /// EPA algorithm implementation for 2D. Only to be used in [`GJK`](struct.GJK.html).
 #[derive(Debug)]
 pub struct EPA2<S> {
-    m: marker::PhantomData<S>,
     tolerance: S,
     max_iterations: u32,
 }
 
-impl<S> EPA for EPA2<S>
+impl<S> EPA2<S>
 where
     S: BaseFloat,
 {
-    type Point = Point2<S>;
-
-    fn process<SL, SR, TL, TR>(
+    fn closest_edge<SL, SR, TL, TR, F>(
         &self,
         simplex: &mut Vec<SupportPoint<Point2<S>>>,
         left: &SL,
         left_transform: &TL,
         right: &SR,
         right_transform: &TR,
-    ) -> Option<Contact<Point2<S>>>
+        closest_edge_fn: F,
+    ) -> Option<Edge<S>>
     where
-        SL: Primitive<Point = Self::Point>,
-        SR: Primitive<Point = Self::Point>,
-        TL: Transform<Self::Point>,
-        TR: Transform<Self::Point>,
+        SL: Primitive<Point = <Self as EPA>::Point>,
+        SR: Primitive<Point = <Self as EPA>::Point>,
+        TL: Transform<<Self as EPA>::Point>,
+        TR: Transform<<Self as EPA>::Point>,
+        F: Fn(&[SupportPoint<Point2<S>>]) -> Option<Edge<S>>,
     {
-        let mut e = match closest_edge(simplex) {
+        let e = closest_edge_fn(simplex);
+        let mut e = match e {
             None => return None,
             Some(e) => e,
         };
@@ -57,9 +58,41 @@ where
             } else {
                 simplex.insert(e.index, p);
             }
-            e = closest_edge(simplex).unwrap();
+            e = closest_edge_fn(simplex)?;
         }
 
+        Some(e)
+    }
+}
+
+impl<S> EPA for EPA2<S>
+where
+    S: BaseFloat,
+{
+    type Point = Point2<S>;
+
+    fn process<SL, SR, TL, TR>(
+        &self,
+        simplex: &mut Vec<SupportPoint<Point2<S>>>,
+        left: &SL,
+        left_transform: &TL,
+        right: &SR,
+        right_transform: &TR,
+    ) -> Option<Contact<Point2<S>>>
+    where
+        SL: Primitive<Point = Self::Point>,
+        SR: Primitive<Point = Self::Point>,
+        TL: Transform<Self::Point>,
+        TR: Transform<Self::Point>,
+    {
+        let e = self.closest_edge(
+            simplex,
+            left,
+            left_transform,
+            right,
+            right_transform,
+            closest_edge,
+        )?;
         Some(Contact::new_with_point(
             CollisionStrategy::FullResolution,
             e.normal,
@@ -77,13 +110,92 @@ where
         max_iterations: u32,
     ) -> Self {
         Self {
-            m: marker::PhantomData,
             tolerance,
             max_iterations,
         }
     }
 }
 
+/// Alternate EPA algorithm implementation for 2D. Only to be used in [`GJK`](struct.GJK.html).
+/// The difference is that the normal returned is guaranteed to be a normal
+/// that exists on the left collider.
+#[derive(Debug)]
+pub struct EPALeft2<S>(EPA2<S>);
+
+impl<S> EPA for EPALeft2<S>
+where
+    S: BaseFloat,
+{
+    type Point = Point2<S>;
+
+    fn process<SL, SR, TL, TR>(
+        &self,
+        simplex: &mut Vec<SupportPoint<Point2<S>>>,
+        left: &SL,
+        left_transform: &TL,
+        right: &SR,
+        right_transform: &TR,
+    ) -> Option<Contact<Point2<S>>>
+    where
+        SL: Primitive<Point = Self::Point>,
+        SR: Primitive<Point = Self::Point>,
+        TL: Transform<Self::Point>,
+        TR: Transform<Self::Point>,
+    {
+        let mut e = self.0.closest_edge(
+            simplex,
+            left,
+            left_transform,
+            right,
+            right_transform,
+            closest_edge,
+        )?;
+        let n = left.closest_valid_normal_local(
+            &left_transform
+                .inverse_transform_vector(e.normal)
+                .unwrap_or(e.normal)
+                .normalize(),
+        );
+
+        if ulps_ne!(e.normal, n) {
+            e = self.0.closest_edge(
+                simplex,
+                left,
+                left_transform,
+                right,
+                right_transform,
+                |s| closest_edge_along_direction(s, n),
+            )?;
+
+            return Some(Contact::new_with_point(
+                CollisionStrategy::FullResolution,
+                n,
+                e.distance / n.dot(e.normal),
+                // TODO: Fix collision point
+                point(simplex, &e),
+            ));
+        }
+
+        Some(Contact::new_with_point(
+            CollisionStrategy::FullResolution,
+            e.normal,
+            e.distance,
+            point(simplex, &e),
+        ))
+    }
+
+    fn new() -> Self {
+        Self::new_with_tolerance(NumCast::from(EPA_TOLERANCE).unwrap(), MAX_ITERATIONS)
+    }
+
+    fn new_with_tolerance(
+        tolerance: <Self::Point as EuclideanSpace>::Scalar,
+        max_iterations: u32,
+    ) -> Self {
+        Self(EPA2::new_with_tolerance(tolerance, max_iterations))
+    }
+}
+
 /// This function returns the contact point in world space coordinates on shape A.
 ///
 /// Compute the closest point to the origin on the given simplex edge, then use that to interpolate
@@ -143,8 +255,7 @@ where
             let a = simplex[i].v;
             let b = simplex[j].v;
             let e = b - a;
-            let oa = a;
-            let n = triple_product(&e, &oa, &e).normalize();
+            let n = vector_normal(&e).normalize();
             let d = n.dot(a);
             if d < edge.distance {
                 edge = Edge::new(n, d, j);
@@ -155,10 +266,49 @@ where
     }
 }
 
+// Assumes direction is normalized
+fn closest_edge_along_direction<S>(
+    simplex: &[SupportPoint<Point2<S>>],
+    direction: Vector2<S>,
+) -> Option<Edge<S>>
+where
+    S: BaseFloat,
+{
+    if simplex.len() < 3 {
+        None
+    } else {
+        let mut edge = Edge::new(Vector2::zero(), S::infinity(), 0);
+        let mut distance = S::infinity();
+
+        for i in 0..simplex.len() {
+            let j = if i + 1 == simplex.len() { 0 } else { i + 1 };
+            let a = simplex[i].v;
+            let b = simplex[j].v;
+            let e = b - a;
+            let n = vector_normal(&e).normalize();
+            let cos = direction.dot(n);
+            if cos <= S::zero() {
+                continue; // Edge is in the opposite direction
+            }
+            let d = n.dot(a);
+            let dist = d / cos;
+            if dist < distance {
+                edge = Edge::new(n, d, j);
+                distance = dist;
+            }
+        }
+        if ulps_ne!(S::infinity(), distance) {
+            Some(edge)
+        } else {
+            None
+        }
+    }
+}
+
 #[cfg(test)]
 mod tests {
-    use cgmath::{Basis2, Decomposed, Point2, Rad, Rotation2, Vector2};
     use approx::assert_ulps_eq;
+    use cgmath::{Basis2, Decomposed, Point2, Rad, Rotation2, Vector2};
 
     use super::*;
     use crate::algorithm::minkowski::SupportPoint;
@@ -196,17 +346,15 @@ mod tests {
         let left_transform = transform(15., 0., 0.);
         let right = Rectangle::new(10., 10.);
         let right_transform = transform(7., 2., 0.);
-        assert!(
-            EPA2::new()
-                .process(
-                    &mut vec![],
-                    &left,
-                    &left_transform,
-                    &right,
-                    &right_transform
-                )
-                .is_none()
-        );
+        assert!(EPA2::new()
+            .process(
+                &mut vec![],
+                &left,
+                &left_transform,
+                &right,
+                &right_transform
+            )
+            .is_none());
     }
 
     #[test]
@@ -216,17 +364,15 @@ mod tests {
         let right = Rectangle::new(10., 10.);
         let right_transform = transform(7., 2., 0.);
         let mut simplex = vec![sup(-2., 8.)];
-        assert!(
-            EPA2::new()
-                .process(
-                    &mut simplex,
-                    &left,
-                    &left_transform,
-                    &right,
-                    &right_transform
-                )
-                .is_none()
-        );
+        assert!(EPA2::new()
+            .process(
+                &mut simplex,
+                &left,
+                &left_transform,
+                &right,
+                &right_transform
+            )
+            .is_none());
     }
 
     #[test]
@@ -236,17 +382,15 @@ mod tests {
         let right = Rectangle::new(10., 10.);
         let right_transform = transform(7., 2., 0.);
         let mut simplex = vec![sup(-2., 8.), sup(18., -12.)];
-        assert!(
-            EPA2::new()
-                .process(
-                    &mut simplex,
-                    &left,
-                    &left_transform,
-                    &right,
-                    &right_transform
-                )
-                .is_none()
-        );
+        assert!(EPA2::new()
+            .process(
+                &mut simplex,
+                &left,
+                &left_transform,
+                &right,
+                &right_transform
+            )
+            .is_none());
     }
 
     #[test]
diff --git a/src/algorithm/minkowski/epa/epa3d.rs b/src/algorithm/minkowski/epa/epa3d.rs
index 4201380..f7d7a40 100644
--- a/src/algorithm/minkowski/epa/epa3d.rs
+++ b/src/algorithm/minkowski/epa/epa3d.rs
@@ -1,14 +1,14 @@
 use std::marker;
 
-use cgmath::{BaseFloat, Point3, Vector3};
-use cgmath::prelude::*;
 use cgmath::num_traits::NumCast;
+use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point3, Vector3};
 
-use super::*;
 use super::SupportPoint;
-use crate::{CollisionStrategy, Contact};
+use super::*;
 use crate::prelude::*;
 use crate::primitive::util::barycentric_vector;
+use crate::{CollisionStrategy, Contact};
 
 /// EPA algorithm implementation for 3D. Only to be used in [`GJK`](struct.GJK.html).
 #[derive(Debug)]
@@ -220,8 +220,8 @@ fn remove_or_add_edge(edges: &mut Vec<(usize, usize)>, edge: (usize, usize)) {
 
 #[cfg(test)]
 mod tests {
-    use cgmath::{Decomposed, Quaternion, Rad, Vector3};
     use approx::assert_ulps_eq;
+    use cgmath::{Decomposed, Quaternion, Rad, Vector3};
 
     use super::*;
     use crate::primitive::*;
@@ -253,25 +253,11 @@ mod tests {
         let faces = Face::new(&simplex);
         assert_eq!(4, faces.len());
         assert_face(
-            &faces[0],
-            3,
-            2,
-            1,
-            -0.8728715,
-            0.43643576,
-            0.21821788,
-            1.0910894,
+            &faces[0], 3, 2, 1, -0.8728715, 0.43643576, 0.21821788, 1.0910894,
         );
         assert_face(&faces[1], 3, 1, 0, 0., -0.89442724, 0.44721362, 2.236068);
         assert_face(
-            &faces[2],
-            3,
-            0,
-            2,
-            0.8728715,
-            0.43643576,
-            0.21821788,
-            1.0910894,
+            &faces[2], 3, 0, 2, 0.8728715, 0.43643576, 0.21821788, 1.0910894,
         );
         assert_face(&faces[3], 2, 0, 1, 0., 0., -1., 1.0);
     }
diff --git a/src/algorithm/minkowski/epa/mod.rs b/src/algorithm/minkowski/epa/mod.rs
index 724ab32..a77861a 100644
--- a/src/algorithm/minkowski/epa/mod.rs
+++ b/src/algorithm/minkowski/epa/mod.rs
@@ -1,5 +1,6 @@
 //! Expanding Polytope Algorithm
 
+pub use self::epa2d::EPALeft2;
 pub use self::epa2d::EPA2;
 pub use self::epa3d::EPA3;
 
@@ -9,8 +10,8 @@ mod epa3d;
 use cgmath::prelude::*;
 
 use super::SupportPoint;
-use crate::Contact;
 use crate::prelude::*;
+use crate::Contact;
 
 pub const EPA_TOLERANCE: f32 = 0.00001;
 pub const MAX_ITERATIONS: u32 = 100;
diff --git a/src/algorithm/minkowski/gjk/mod.rs b/src/algorithm/minkowski/gjk/mod.rs
index 5069efc..67e6895 100644
--- a/src/algorithm/minkowski/gjk/mod.rs
+++ b/src/algorithm/minkowski/gjk/mod.rs
@@ -6,15 +6,15 @@ pub use self::simplex::SimplexProcessor;
 use std::cmp::Ordering;
 use std::ops::{Neg, Range};
 
-use cgmath::BaseFloat;
-use cgmath::prelude::*;
 use cgmath::num_traits::NumCast;
+use cgmath::prelude::*;
+use cgmath::BaseFloat;
 use cgmath::UlpsEq;
 
 use self::simplex::{Simplex, SimplexProcessor2, SimplexProcessor3};
-use crate::{CollisionStrategy, Contact};
-use crate::algorithm::minkowski::{EPA2, EPA3, SupportPoint, EPA};
+use crate::algorithm::minkowski::{EPALeft2, SupportPoint, EPA, EPA2, EPA3};
 use crate::prelude::*;
+use crate::{CollisionStrategy, Contact};
 use approx::ulps_eq;
 
 mod simplex;
@@ -26,6 +26,11 @@ const GJK_CONTINUOUS_TOLERANCE: f32 = 0.000001;
 /// GJK algorithm for 2D, see [GJK](struct.GJK.html) for more information.
 pub type GJK2<S> = GJK<SimplexProcessor2<S>, EPA2<S>, S>;
 
+/// GJK algorithm for 2D, see [GJK](struct.GJK.html) for more information.
+/// This one guarantees that the normal returned in the case of full resolution
+/// is from the left collider.
+pub type GJKLeft2<S> = GJK<SimplexProcessor2<S>, EPALeft2<S>, S>;
+
 /// GJK algorithm for 3D, see [GJK](struct.GJK.html) for more information.
 pub type GJK3<S> = GJK<SimplexProcessor3<S>, EPA3<S>, S>;
 
@@ -133,7 +138,8 @@ where
                 return None;
             } else {
                 simplex.push(a);
-                if self.simplex_processor
+                if self
+                    .simplex_processor
                     .reduce_to_closest_feature(&mut simplex, &mut d)
                 {
                     return Some(simplex);
@@ -233,7 +239,8 @@ where
             }
             // we construct the simplex around the current ray origin (if we can)
             simplex.push(p - ray_origin);
-            v = self.simplex_processor
+            v = self
+                .simplex_processor
                 .get_closest_point_to_origin(&mut simplex);
         }
         if v.magnitude2() <= self.continuous_tolerance {
@@ -300,7 +307,8 @@ where
             ));
         }
         for _ in 0..self.max_iterations {
-            let d = self.simplex_processor
+            let d = self
+                .simplex_processor
                 .get_closest_point_to_origin(&mut simplex);
             if ulps_eq!(d, zero) {
                 return None;
@@ -587,9 +595,10 @@ where
 
 #[cfg(test)]
 mod tests {
-    use cgmath::{Basis2, Decomposed, Point2, Point3, Quaternion, Rad, Rotation2, Rotation3,
-                 Vector2, Vector3};
     use approx::assert_ulps_eq;
+    use cgmath::{
+        Basis2, Decomposed, Point2, Point3, Quaternion, Rad, Rotation2, Rotation3, Vector2, Vector3,
+    };
 
     use super::*;
     use crate::primitive::*;
@@ -655,17 +664,18 @@ mod tests {
         let right = Rectangle::new(10., 10.);
         let right_transform = transform(-15., 0., 0.);
         let gjk = GJK2::new();
-        assert!(
-            gjk.intersect(&left, &left_transform, &right, &right_transform)
-                .is_none()
-        );
-        assert!(gjk.intersection(
-            &CollisionStrategy::FullResolution,
-            &left,
-            &left_transform,
-            &right,
-            &right_transform
-        ).is_none())
+        assert!(gjk
+            .intersect(&left, &left_transform, &right, &right_transform)
+            .is_none());
+        assert!(gjk
+            .intersection(
+                &CollisionStrategy::FullResolution,
+                &left,
+                &left_transform,
+                &right,
+                &right_transform
+            )
+            .is_none())
     }
 
     #[test]
@@ -691,6 +701,83 @@ mod tests {
         assert_eq!(Point2::new(10., 1.), contact.contact_point);
     }
 
+    #[test]
+    fn test_gjk_either_2d_hit() {
+        // Make sure only GJKLeft is changed
+        let left = Rectangle::new(2.0, 2.0);
+        let left_transform = transform(0., 0., 0.);
+        let right = ConvexPolygon::new(vec![
+            Point2::new(-2.0, 0.0),
+            Point2::new(0.0, -1.0),
+            Point2::new(2.0, 0.0),
+        ]);
+        let right_transform = transform(-1.5, 1.5, 0.);
+        let gjk = GJK2::new();
+        let contact = gjk.intersection(
+            &CollisionStrategy::FullResolution,
+            &left,
+            &left_transform,
+            &right,
+            &right_transform,
+        );
+        assert!(contact.is_some());
+        let contact = contact.unwrap();
+        // This is a normal GJK. The normal should point slightly to the left.
+        assert!(contact.normal.x < -0.4);
+    }
+
+    #[test]
+    fn test_gjk_left_2d_hit() {
+        let left = Rectangle::new(2.0, 2.0);
+        let left_transform = transform(0., 0., 0.);
+        let right = ConvexPolygon::new(vec![
+            Point2::new(-2.0, 0.0),
+            Point2::new(0.0, -1.0),
+            Point2::new(2.0, 0.0),
+        ]);
+        let right_transform = transform(-1.5, 1.5, 0.);
+        let gjk = GJKLeft2::new();
+        let contact = gjk.intersection(
+            &CollisionStrategy::FullResolution,
+            &left,
+            &left_transform,
+            &right,
+            &right_transform,
+        );
+        assert!(contact.is_some());
+        let contact = contact.unwrap();
+        // If this were a normal GJK2, the normal would point slightly to the left
+        assert_eq!(Vector2::new(0., 1.), contact.normal);
+        assert_ulps_eq!(0.25, contact.penetration_depth);
+    }
+
+    #[test]
+    fn test_gjk_left_2d_hit_vertex() {
+        // Case where the penetration causes a vertex of the triangle
+        // to hit the rectangle instead of an edge.
+        // Penetration depth is especially important to check.
+        let left = Rectangle::new(2.0, 2.0);
+        let left_transform = transform(0., 0., 0.);
+        let right = ConvexPolygon::new(vec![
+            Point2::new(-2.0, 0.0),
+            Point2::new(0.0, -1.0),
+            Point2::new(2.0, 0.0),
+        ]);
+        let right_transform = transform(-0.99, 1.5, 0.);
+        let gjk = GJKLeft2::new();
+        let contact = gjk.intersection(
+            &CollisionStrategy::FullResolution,
+            &left,
+            &left_transform,
+            &right,
+            &right_transform,
+        );
+        assert!(contact.is_some());
+        let contact = contact.unwrap();
+        assert_eq!(Vector2::new(0., 1.), contact.normal);
+        assert_ulps_eq!(0.5, contact.penetration_depth);
+    }
+
     #[test]
     fn test_gjk_3d_hit() {
         let left = Cuboid::new(10., 10., 10.);
@@ -714,6 +801,127 @@ mod tests {
         assert_ulps_eq!(Point3::new(10., 1., 5.), contact.contact_point);
     }
 
+    /// Test for [issue #115](https://github.com/rustgd/collision-rs/issues/115)
+    #[test]
+    fn test_gjk_hit_rounding_edge_case() {
+        use cgmath::Matrix3;
+
+        let left = Rectangle::new(0.75, 1.);
+        // We must use the full power of the f64 to replicate the bug.
+        // Not using `transform()` because it takes `f32`s
+        let left_transform =
+            Matrix3::new(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.858333333333344, 0.0, 1.0);
+        let right = ConvexPolygon::new(vec![
+            Point2::new(0.5, 0.5),
+            Point2::new(-0.5, 0.5),
+            Point2::new(-0.5, -0.5),
+        ]);
+        let right_transform = Matrix3::new(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0);
+        let gjk = GJK2::new();
+        let contact = gjk.intersection(
+            &CollisionStrategy::FullResolution,
+            &left,
+            &left_transform,
+            &right,
+            &right_transform,
+        );
+        assert!(contact.is_some());
+        let contact = contact.unwrap();
+        assert!(
+            contact.penetration_depth < 0.5,
+            format!(
+                "Penetration depth is {}, which is too big",
+                contact.penetration_depth
+            )
+        );
+    }
+    #[test]
+    fn test_gjk_hit_regression_1() {
+        use cgmath::Matrix3;
+
+        let left = ConvexPolygon::new(vec![
+            Point2::new(0.5, -0.5),
+            Point2::new(0.5, 0.5),
+            Point2::new(-0.5, 0.5),
+        ]);
+        let left_transform = Matrix3::new(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 12.5, 6.5, 1.0);
+        let right = ConvexPolygon::new(vec![
+            Point2::new(0.375, -1.0),
+            Point2::new(0.375, -0.01),
+            Point2::new(0.0, 0.0),
+            Point2::new(-0.375, -0.01),
+            Point2::new(-0.375, -1.0),
+        ]);
+        let right_transform = Matrix3::new(
+            1.0,
+            0.0,
+            0.0,
+            0.0,
+            1.0,
+            0.0,
+            12.404999999999987,
+            6.230000000000013,
+            1.0,
+        );
+        let gjk = GJK2::new();
+        let contact = gjk.intersection(
+            &CollisionStrategy::FullResolution,
+            &left,
+            &left_transform,
+            &right,
+            &right_transform,
+        );
+        assert!(contact.is_some());
+        let contact = contact.unwrap();
+        assert_ulps_eq!(contact.normal, Vector2::new(-0.5f64.sqrt(), -0.5f64.sqrt()));
+        assert!(
+            contact.penetration_depth < 0.1,
+            format!(
+                "Penetration depth is {}, which is too big",
+                contact.penetration_depth
+            )
+        );
+    }
+
+    #[test]
+    fn test_gjk_hit_regression_2() {
+        use cgmath::Matrix3;
+
+        let left = ConvexPolygon::new(vec![
+            Point2::new(0.5, 0.5),
+            Point2::new(-0.5, 0.5),
+            Point2::new(-0.5, -0.5),
+        ]);
+        let left_transform = Matrix3::new(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 13.5, 6.5, 1.0);
+        let right = ConvexPolygon::new(vec![
+            Point2::new(0.375, -1.0),
+            Point2::new(0.375, -0.01),
+            Point2::new(0.0, 0.0),
+            Point2::new(-0.375, -0.01),
+            Point2::new(-0.375, -1.0),
+        ]);
+        let right_transform = Matrix3::new(
+            1.0,
+            0.0,
+            0.0,
+            0.0,
+            1.0,
+            0.0,
+            13.585769040132476,
+            6.154102373465809,
+            1.0,
+        );
+        let gjk = GJK2::new();
+        let contact = gjk.intersection(
+            &CollisionStrategy::FullResolution,
+            &left,
+            &left_transform,
+            &right,
+            &right_transform,
+        );
+        assert!(contact.is_none());
+    }
+
     #[test]
     fn test_gjk_distance_2d() {
         let left = Rectangle::new(10., 10.);
@@ -762,24 +970,28 @@ mod tests {
         let right_transform = transform(15., 0., 0.);
         let gjk = GJK2::new();
 
-        let contact = gjk.intersection_time_of_impact(
-            &left,
-            &left_start_transform..&left_end_transform,
-            &right,
-            &right_transform..&right_transform,
-        ).unwrap();
+        let contact = gjk
+            .intersection_time_of_impact(
+                &left,
+                &left_start_transform..&left_end_transform,
+                &right,
+                &right_transform..&right_transform,
+            )
+            .unwrap();
 
         assert_ulps_eq!(0.1666667, contact.time_of_impact);
         assert_eq!(Vector2::new(-1., 0.), contact.normal);
         assert_eq!(0., contact.penetration_depth);
         assert_eq!(Point2::new(10., 0.), contact.contact_point);
 
-        assert!(gjk.intersection_time_of_impact(
-            &left,
-            &left_start_transform..&left_start_transform,
-            &right,
-            &right_transform..&right_transform
-        ).is_none());
+        assert!(gjk
+            .intersection_time_of_impact(
+                &left,
+                &left_start_transform..&left_start_transform,
+                &right,
+                &right_transform..&right_transform
+            )
+            .is_none());
     }
 
     #[test]
@@ -791,23 +1003,27 @@ mod tests {
         let right_transform = transform_3d(15., 0., 0., 0.);
         let gjk = GJK3::new();
 
-        let contact = gjk.intersection_time_of_impact(
-            &left,
-            &left_start_transform..&left_end_transform,
-            &right,
-            &right_transform..&right_transform,
-        ).unwrap();
+        let contact = gjk
+            .intersection_time_of_impact(
+                &left,
+                &left_start_transform..&left_end_transform,
+                &right,
+                &right_transform..&right_transform,
+            )
+            .unwrap();
 
         assert_ulps_eq!(0.1666667, contact.time_of_impact);
         assert_eq!(Vector3::new(-1., 0., 0.), contact.normal);
         assert_eq!(0., contact.penetration_depth);
         assert_eq!(Point3::new(10., 0., 0.), contact.contact_point);
 
-        assert!(gjk.intersection_time_of_impact(
-            &left,
-            &left_start_transform..&left_start_transform,
-            &right,
-            &right_transform..&right_transform
-        ).is_none());
+        assert!(gjk
+            .intersection_time_of_impact(
+                &left,
+                &left_start_transform..&left_start_transform,
+                &right,
+                &right_transform..&right_transform
+            )
+            .is_none());
     }
 }
diff --git a/src/algorithm/minkowski/gjk/simplex/simplex2d.rs b/src/algorithm/minkowski/gjk/simplex/simplex2d.rs
index f37dd57..cd0144c 100644
--- a/src/algorithm/minkowski/gjk/simplex/simplex2d.rs
+++ b/src/algorithm/minkowski/gjk/simplex/simplex2d.rs
@@ -1,9 +1,9 @@
 use std::marker;
 use std::ops::Neg;
 
-use cgmath::{BaseFloat, Point2, Vector2};
-use cgmath::prelude::*;
 use approx::ulps_eq;
+use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point2, Vector2};
 
 use super::{Simplex, SimplexProcessor};
 use crate::primitive::util::{get_closest_point_on_edge, triple_product};
@@ -58,7 +58,18 @@ where
             let ab = b - a;
 
             *d = triple_product(&ab, &ao, &ab);
-            if ulps_eq!(*d, Vector2::zero()) {
+            // Product can be dangerously close to 0 and fail
+            // the ulps_eq! check if the origin is on the
+            // line segment connecting a and b
+            //
+            // The maximum error is determined to be around
+            // 2 * ab.magnitude() * ao.magnitude() * ab.magnitude() * S::EPSILON
+            // at least for normal floating-point numbers
+            //
+            // We approximate ao.magnitude() as 1.5 * <max coord of ao>
+            let error =
+                S::from(3.0).unwrap() * ab.magnitude2() * ao.x.abs().max(ao.y.abs()) * S::epsilon();
+            if d.x.abs() < error && d.y.abs() < error {
                 *d = Vector2::new(-ab.y, ab.x);
             }
         }
diff --git a/src/algorithm/minkowski/gjk/simplex/simplex3d.rs b/src/algorithm/minkowski/gjk/simplex/simplex3d.rs
index f6c4090..9bedd03 100644
--- a/src/algorithm/minkowski/gjk/simplex/simplex3d.rs
+++ b/src/algorithm/minkowski/gjk/simplex/simplex3d.rs
@@ -1,10 +1,10 @@
 use std::marker;
 use std::ops::Neg;
 
-use cgmath::{BaseFloat, Point3, Vector3};
-use cgmath::prelude::*;
-use cgmath::num_traits::cast;
 use approx::ulps_eq;
+use cgmath::num_traits::cast;
+use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point3, Vector3};
 
 use super::{Simplex, SimplexProcessor};
 use crate::primitive::util::{barycentric_vector, get_closest_point_on_edge};
@@ -223,8 +223,8 @@ fn check_side<S>(
 mod tests {
     use std::ops::Neg;
 
-    use cgmath::{Point3, Vector3};
     use approx::assert_ulps_eq;
+    use cgmath::{Point3, Vector3};
 
     use super::*;
     use crate::algorithm::minkowski::SupportPoint;
diff --git a/src/algorithm/minkowski/mod.rs b/src/algorithm/minkowski/mod.rs
index eef3aad..23adaeb 100644
--- a/src/algorithm/minkowski/mod.rs
+++ b/src/algorithm/minkowski/mod.rs
@@ -1,12 +1,12 @@
 //! Algorithms using the Minkowski Sum/Difference
 
-pub use self::epa::{EPA2, EPA3, EPA};
-pub use self::gjk::{GJK2, GJK3, SimplexProcessor, GJK};
+pub use self::epa::{EPALeft2, EPA, EPA2, EPA3};
+pub use self::gjk::{GJKLeft2, SimplexProcessor, GJK, GJK2, GJK3};
 
 use std::ops::{Neg, Sub};
 
-use cgmath::prelude::*;
 use crate::prelude::*;
+use cgmath::prelude::*;
 
 mod epa;
 mod gjk;
@@ -105,7 +105,8 @@ mod tests {
                 &right,
                 &right_transform,
                 &direction,
-            ).v
+            )
+            .v
         );
     }
 }
diff --git a/src/algorithm/mod.rs b/src/algorithm/mod.rs
index 030bcbd..57a2720 100644
--- a/src/algorithm/mod.rs
+++ b/src/algorithm/mod.rs
@@ -1,4 +1,4 @@
 //! Collision detection algorithms
 
-pub mod minkowski;
 pub mod broad_phase;
+pub mod minkowski;
diff --git a/src/bound.rs b/src/bound.rs
index ef73c68..fe95e98 100644
--- a/src/bound.rs
+++ b/src/bound.rs
@@ -1,8 +1,8 @@
 //! Generic spatial bounds.
 
-use cgmath::{EuclideanSpace, Point3};
 use cgmath::BaseFloat;
 use cgmath::Matrix4;
+use cgmath::{EuclideanSpace, Point3};
 use std::{cmp, fmt};
 
 use crate::frustum::Frustum;
@@ -37,14 +37,15 @@ pub trait PlaneBound<S: BaseFloat>: fmt::Debug {
             frustum.bottom,
             frustum.near,
             frustum.far,
-        ].iter()
-            .fold(Relation::In, |cur, p| {
-                let r = self.relate_plane(*p);
-                // If any of the planes are `Out`, the bound is outside.
-                // Otherwise, if any are `Cross`, the bound is crossing.
-                // Otherwise, the bound is fully inside.
-                cmp::max(cur, r)
-            })
+        ]
+        .iter()
+        .fold(Relation::In, |cur, p| {
+            let r = self.relate_plane(*p);
+            // If any of the planes are `Out`, the bound is outside.
+            // Otherwise, if any are `Cross`, the bound is crossing.
+            // Otherwise, the bound is fully inside.
+            cmp::max(cur, r)
+        })
     }
 }
 
diff --git a/src/dbvt/mod.rs b/src/dbvt/mod.rs
index f9c8167..3bee1a8 100644
--- a/src/dbvt/mod.rs
+++ b/src/dbvt/mod.rs
@@ -100,9 +100,9 @@ use rand::Rng;
 
 use crate::prelude::*;
 
-mod wrapped;
-mod visitor;
 mod util;
+mod visitor;
+mod wrapped;
 
 const SURFACE_AREA_IMPROVEMENT_FOR_ROTATION: f32 = 0.3;
 const PERFORM_ROTATION_PERCENTAGE: u32 = 10;
@@ -505,11 +505,13 @@ where
 
                 // if we encounter a branch, do intersection test, and push the children if the
                 // branch intersected
-                Node::Branch(ref branch) => if visitor.accept(&branch.bound, false).is_some() {
-                    stack[stack_pointer] = branch.left;
-                    stack[stack_pointer + 1] = branch.right;
-                    stack_pointer += 2;
-                },
+                Node::Branch(ref branch) => {
+                    if visitor.accept(&branch.bound, false).is_some() {
+                        stack[stack_pointer] = branch.left;
+                        stack[stack_pointer + 1] = branch.right;
+                        stack_pointer += 2;
+                    }
+                }
                 Node::Nil => (),
             }
         }
@@ -561,7 +563,8 @@ where
     /// all insert/remove/updates have been performed this frame.
     ///
     pub fn update(&mut self) {
-        let nodes = self.updated_list
+        let nodes = self
+            .updated_list
             .iter()
             .filter_map(|&index| {
                 if let Node::Leaf(ref l) = self.nodes[index] {
diff --git a/src/dbvt/util.rs b/src/dbvt/util.rs
index e2c9fc4..77cb562 100644
--- a/src/dbvt/util.rs
+++ b/src/dbvt/util.rs
@@ -4,12 +4,12 @@
 
 use std::marker::PhantomData;
 
-use cgmath::BaseFloat;
 use cgmath::prelude::*;
+use cgmath::BaseFloat;
 
 use super::{ContinuousVisitor, DynamicBoundingVolumeTree, TreeValue, Visitor};
-use crate::Ray;
 use crate::prelude::*;
+use crate::Ray;
 
 struct RayClosestVisitor<S, P, T>
 where
diff --git a/src/dbvt/visitor.rs b/src/dbvt/visitor.rs
index 9c1121d..a22d10f 100644
--- a/src/dbvt/visitor.rs
+++ b/src/dbvt/visitor.rs
@@ -7,8 +7,8 @@ use std::marker::PhantomData;
 use cgmath::BaseFloat;
 
 use super::{TreeValue, Visitor};
-use crate::{Frustum, PlaneBound, Relation};
 use crate::prelude::*;
+use crate::{Frustum, PlaneBound, Relation};
 
 /// Visitor for doing continuous intersection testing on the DBVT.
 ///
diff --git a/src/frustum.rs b/src/frustum.rs
index 4031693..d6518fc 100644
--- a/src/frustum.rs
+++ b/src/frustum.rs
@@ -1,11 +1,11 @@
 //! View frustum for visibility determination
 
-use crate::Plane;
 use crate::bound::*;
-use cgmath::{Matrix, Matrix4};
-use cgmath::{Ortho, Perspective, PerspectiveFov};
+use crate::Plane;
 use cgmath::BaseFloat;
 use cgmath::Point3;
+use cgmath::{Matrix, Matrix4};
+use cgmath::{Ortho, Perspective, PerspectiveFov};
 
 /// View frustum, used for frustum culling
 #[derive(Copy, Clone, Debug, PartialEq)]
@@ -84,15 +84,16 @@ impl<S: BaseFloat> Frustum<S> {
             self.bottom,
             self.near,
             self.far,
-        ].iter()
-            .fold(Relation::In, |cur, p| {
-                use std::cmp::max;
-                let r = bound.relate_plane(*p);
-                // If any of the planes are `Out`, the bound is outside.
-                // Otherwise, if any are `Cross`, the bound is crossing.
-                // Otherwise, the bound is fully inside.
-                max(cur, r)
-            })
+        ]
+        .iter()
+        .fold(Relation::In, |cur, p| {
+            use std::cmp::max;
+            let r = bound.relate_plane(*p);
+            // If any of the planes are `Out`, the bound is outside.
+            // Otherwise, if any are `Cross`, the bound is crossing.
+            // Otherwise, the bound is fully inside.
+            max(cur, r)
+        })
     }
 }
 
diff --git a/src/lib.rs b/src/lib.rs
index aaf393d..9ffce96 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,7 +1,13 @@
 #![crate_type = "rlib"]
 #![crate_type = "dylib"]
-#![deny(missing_docs, trivial_casts, unsafe_code, unstable_features, unused_import_braces,
-unused_qualifications)]
+#![deny(
+    missing_docs,
+    trivial_casts,
+    unsafe_code,
+    unstable_features,
+    unused_import_braces,
+    unused_qualifications
+)]
 
 //! Companion library to cgmath, dealing with collision detection centric data structures and
 //! algorithms.
@@ -22,8 +28,6 @@ extern crate serde;
 #[cfg_attr(test, macro_use)]
 extern crate smallvec;
 
-
-
 // Re-exports
 
 pub use bound::*;
@@ -35,18 +39,18 @@ pub use ray::*;
 pub use traits::*;
 pub use volume::*;
 
-pub mod prelude;
+pub mod algorithm;
 pub mod dbvt;
+pub mod prelude;
 pub mod primitive;
-pub mod algorithm;
 
 // Modules
 
 mod bound;
+mod contact;
 mod frustum;
-mod traits;
+mod line;
 mod plane;
 mod ray;
-mod line;
+mod traits;
 mod volume;
-mod contact;
diff --git a/src/line.rs b/src/line.rs
index 894e385..1a16c3e 100644
--- a/src/line.rs
+++ b/src/line.rs
@@ -2,13 +2,13 @@
 
 use std::marker::PhantomData;
 
+use cgmath::prelude::*;
 use cgmath::{BaseFloat, BaseNum};
 use cgmath::{Point2, Point3};
 use cgmath::{Vector2, Vector3};
-use cgmath::prelude::*;
 
-use crate::Ray2;
 use crate::prelude::*;
+use crate::Ray2;
 
 /// A generic directed line segment from `origin` to `dest`.
 #[derive(Copy, Clone, PartialEq, Debug)]
diff --git a/src/plane.rs b/src/plane.rs
index 48e035b..9b57d81 100644
--- a/src/plane.rs
+++ b/src/plane.rs
@@ -1,12 +1,12 @@
 use std::fmt;
 
-use cgmath::{BaseFloat, Point3, Vector3, Vector4};
-use cgmath::{AbsDiffEq, RelativeEq, UlpsEq};
-use cgmath::prelude::*;
 use approx::{ulps_eq, ulps_ne};
+use cgmath::prelude::*;
+use cgmath::{AbsDiffEq, RelativeEq, UlpsEq};
+use cgmath::{BaseFloat, Point3, Vector3, Vector4};
 
-use crate::Ray3;
 use crate::prelude::*;
+use crate::Ray3;
 
 /// A 3-dimensional plane formed from the equation: `A*x + B*y + C*z - D = 0`.
 ///
@@ -108,9 +108,9 @@ impl<S: BaseFloat> Plane<S> {
 }
 
 impl<S: AbsDiffEq> AbsDiffEq for Plane<S>
-    where
-        S::Epsilon: Copy,
-        S: BaseFloat,
+where
+    S::Epsilon: Copy,
+    S: BaseFloat,
 {
     type Epsilon = S::Epsilon;
 
@@ -127,9 +127,9 @@ impl<S: AbsDiffEq> AbsDiffEq for Plane<S>
 }
 
 impl<S: RelativeEq> RelativeEq for Plane<S>
-    where
-        S::Epsilon: Copy,
-        S: BaseFloat,
+where
+    S::Epsilon: Copy,
+    S: BaseFloat,
 {
     #[inline]
     fn default_max_relative() -> S::Epsilon {
@@ -143,9 +143,9 @@ impl<S: RelativeEq> RelativeEq for Plane<S>
     }
 }
 impl<S: UlpsEq> UlpsEq for Plane<S>
-    where
-        S::Epsilon: Copy,
-        S: BaseFloat,
+where
+    S::Epsilon: Copy,
+    S: BaseFloat,
 {
     #[inline]
     fn default_max_ulps() -> u32 {
diff --git a/src/primitive/capsule.rs b/src/primitive/capsule.rs
index 03f58d1..719133b 100644
--- a/src/primitive/capsule.rs
+++ b/src/primitive/capsule.rs
@@ -1,10 +1,10 @@
-use cgmath::{BaseFloat, Point3, Vector3};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point3, Vector3};
 
-use crate::{Aabb3, Ray3};
 use crate::prelude::*;
 use crate::primitive::util::cylinder_ray_quadratic_solve;
 use crate::volume::Sphere;
+use crate::{Aabb3, Ray3};
 
 /// Capsule primitive
 /// Capsule body is aligned with the Y axis, with local origin in the center of the capsule.
@@ -54,6 +54,13 @@ where
         result.y = direction.y.signum() * self.half_height;
         transform.transform_point(result + direction.normalize_to(self.radius))
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        unimplemented!("closest_valid_normal_local is only implemented for 2D primitives for now")
+    }
 }
 
 impl<S> ComputeBound<Aabb3<S>> for Capsule<S>
@@ -201,8 +208,8 @@ where
 mod tests {
     use std;
 
-    use cgmath::{Decomposed, Quaternion, Rad, Vector3};
     use approx::assert_ulps_eq;
+    use cgmath::{Decomposed, Quaternion, Rad, Vector3};
 
     use super::*;
 
diff --git a/src/primitive/circle.rs b/src/primitive/circle.rs
index efc6653..f49fdc0 100644
--- a/src/primitive/circle.rs
+++ b/src/primitive/circle.rs
@@ -1,10 +1,10 @@
 //! Circle primitive
 
-use cgmath::{BaseFloat, Point2, Vector2};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point2, Vector2};
 
-use crate::{Aabb2, Ray2};
 use crate::prelude::*;
+use crate::{Aabb2, Ray2};
 
 /// Circle primitive
 #[derive(Debug, Clone, PartialEq)]
@@ -34,6 +34,13 @@ where
         let direction = transform.inverse_transform_vector(*direction).unwrap();
         transform.transform_point(Point2::from_vec(direction.normalize_to(self.radius)))
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        *normal
+    }
 }
 
 impl<S> ComputeBound<Aabb2<S>> for Circle<S>
@@ -91,10 +98,10 @@ where
 mod tests {
     use std;
 
-    use crate::Ray2;
-    use cgmath::{Basis2, Decomposed, Point2, Rad, Rotation2, Vector2};
     use crate::prelude::*;
+    use crate::Ray2;
     use approx::assert_ulps_eq;
+    use cgmath::{vec2, Basis2, Decomposed, Point2, Rad, Rotation2, Vector2};
 
     use super::*;
 
@@ -129,6 +136,13 @@ mod tests {
         assert_eq!(bound(-10., -10., 10., 10.), circle.compute_bound())
     }
 
+    #[test]
+    fn test_circle_closest_valid_normal() {
+        let circle = Circle::new(1.);
+        let vector = vec2(0.5f64.sqrt(), -0.5f64.sqrt());
+        assert_eq!(vector, circle.closest_valid_normal_local(&vector));
+    }
+
     #[test]
     fn test_circle_ray_discrete() {
         let circle = Circle::new(10.);
diff --git a/src/primitive/cuboid.rs b/src/primitive/cuboid.rs
index a0e448f..cd6ecdb 100644
--- a/src/primitive/cuboid.rs
+++ b/src/primitive/cuboid.rs
@@ -1,10 +1,10 @@
-use cgmath::{BaseFloat, Point3, Vector3};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point3, Vector3};
 
-use crate::{Aabb3, Ray3};
 use crate::prelude::*;
 use crate::primitive::util::get_max_point;
 use crate::volume::Sphere;
+use crate::{Aabb3, Ray3};
 
 /// Cuboid primitive.
 ///
@@ -73,6 +73,13 @@ where
     {
         get_max_point(self.corners.iter(), direction, transform)
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        unimplemented!("closest_valid_normal_local is only implemented for 2D primitives for now")
+    }
 }
 
 impl<S> ComputeBound<Aabb3<S>> for Cuboid<S>
@@ -108,7 +115,8 @@ where
         Aabb3::new(
             Point3::from_vec(-self.half_dim),
             Point3::from_vec(self.half_dim),
-        ).intersects(ray)
+        )
+        .intersects(ray)
     }
 }
 
@@ -122,7 +130,8 @@ where
         Aabb3::new(
             Point3::from_vec(-self.half_dim),
             Point3::from_vec(self.half_dim),
-        ).intersection(ray)
+        )
+        .intersection(ray)
     }
 }
 
@@ -169,6 +178,13 @@ where
     {
         self.cuboid.support_point(direction, transform)
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        unimplemented!("closest_valid_normal_local is only implemented for 2D primitives for now")
+    }
 }
 
 impl<S> ComputeBound<Aabb3<S>> for Cube<S>
@@ -212,8 +228,8 @@ where
 #[cfg(test)]
 mod tests {
 
-    use cgmath::{Decomposed, Point3, Quaternion, Rad, Vector3};
     use approx::assert_ulps_eq;
+    use cgmath::{Decomposed, Point3, Quaternion, Rad, Vector3};
 
     use super::*;
     use Ray3;
diff --git a/src/primitive/cylinder.rs b/src/primitive/cylinder.rs
index 7fa7673..f991c04 100644
--- a/src/primitive/cylinder.rs
+++ b/src/primitive/cylinder.rs
@@ -1,10 +1,10 @@
-use cgmath::{BaseFloat, Point3, Vector3};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point3, Vector3};
 
-use crate::{Aabb3, Ray3};
 use crate::prelude::*;
 use crate::primitive::util::cylinder_ray_quadratic_solve;
 use crate::volume::Sphere;
+use crate::{Aabb3, Ray3};
 
 /// Cylinder primitive
 /// Cylinder body is aligned with the Y axis, with local origin in the center of the cylinders.
@@ -71,6 +71,13 @@ where
         }
         transform.transform_point(Point3::from_vec(result))
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        unimplemented!("closest_valid_normal_local is only implemented for 2D primitives for now")
+    }
 }
 
 impl<S> ComputeBound<Aabb3<S>> for Cylinder<S>
@@ -231,8 +238,8 @@ where
 mod tests {
     use std;
 
-    use cgmath::{Decomposed, Quaternion, Rad, Vector3};
     use approx::assert_ulps_eq;
+    use cgmath::{Decomposed, Quaternion, Rad, Vector3};
 
     use super::*;
 
diff --git a/src/primitive/line.rs b/src/primitive/line.rs
index bac28e2..5df556a 100644
--- a/src/primitive/line.rs
+++ b/src/primitive/line.rs
@@ -1,8 +1,8 @@
-use cgmath::{BaseFloat, InnerSpace, Point2, Transform, Vector2};
+use cgmath::{vec2, BaseFloat, EuclideanSpace, InnerSpace, Point2, Transform, Vector2, Zero};
 
-use crate::Aabb2;
 use crate::line::Line2;
 use crate::traits::{ComputeBound, Primitive};
+use crate::Aabb2;
 
 impl<S> Primitive for Line2<S>
 where
@@ -22,6 +22,19 @@ where
             transform.transform_point(self.origin)
         }
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        let mut perp = (self.dest - self.origin).normalize();
+        perp = vec2(-perp.y, perp.x);
+        if normal.dot(perp) >= Zero::zero() {
+            perp
+        } else {
+            -perp
+        }
+    }
 }
 
 impl<S> ComputeBound<Aabb2<S>> for Line2<S>
@@ -38,8 +51,8 @@ mod tests {
 
     use super::*;
     use crate::algorithm::minkowski::GJK2;
-    use cgmath::{Basis2, Decomposed, Rad, Rotation2};
     use crate::primitive::Rectangle;
+    use cgmath::{Basis2, Decomposed, Rad, Rotation2};
 
     fn transform(x: f32, y: f32, angle: f32) -> Decomposed<Vector2<f32>, Basis2<f32>> {
         Decomposed {
@@ -49,6 +62,19 @@ mod tests {
         }
     }
 
+    #[test]
+    fn test_line_closest_valid_normal() {
+        let line = Line2::new(Point2::new(1., 1.), Point2::new(1., 2.));
+        assert_eq!(
+            vec2(-1., 0.),
+            line.closest_valid_normal_local(&vec2(-0.5, 0.75f64.sqrt()))
+        );
+        assert_eq!(
+            vec2(1., 0.),
+            line.closest_valid_normal_local(&vec2(0.5, 0.75f64.sqrt()))
+        );
+    }
+
     #[test]
     fn test_line_rectangle_intersect() {
         let line = Line2::new(Point2::new(0., -1.), Point2::new(0., 1.));
@@ -56,9 +82,8 @@ mod tests {
         let transform_1 = transform(1., 0., 0.);
         let transform_2 = transform(1.1, 0., 0.);
         let gjk = GJK2::new();
-        assert!(
-            gjk.intersect(&line, &transform_1, &rectangle, &transform_2)
-                .is_some()
-        );
+        assert!(gjk
+            .intersect(&line, &transform_1, &rectangle, &transform_2)
+            .is_some());
     }
 }
diff --git a/src/primitive/mod.rs b/src/primitive/mod.rs
index d8140c6..faea124 100644
--- a/src/primitive/mod.rs
+++ b/src/primitive/mod.rs
@@ -13,18 +13,18 @@ pub use self::quad::Quad;
 pub use self::rectangle::{Rectangle, Square};
 pub use self::sphere::Sphere;
 
-mod circle;
-mod cylinder;
 mod capsule;
+mod circle;
 mod cuboid;
+mod cylinder;
 mod line;
 mod particle;
 mod polygon;
 mod polyhedron;
 mod primitive2;
 mod primitive3;
-mod rectangle;
 mod quad;
+mod rectangle;
 mod sphere;
 
 pub(crate) mod util;
@@ -62,4 +62,11 @@ where
     {
         self.0.support_point(direction, transform)
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        self.0.closest_valid_normal_local(normal)
+    }
 }
diff --git a/src/primitive/particle.rs b/src/primitive/particle.rs
index 8b9dfc9..e668f37 100644
--- a/src/primitive/particle.rs
+++ b/src/primitive/particle.rs
@@ -3,11 +3,11 @@
 use std::marker;
 use std::ops::Range;
 
-use cgmath::{BaseFloat, Point2, Point3};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point2, Point3};
 
-use crate::Ray;
 use crate::prelude::*;
+use crate::Ray;
 
 /// Represents a particle in space.
 ///
@@ -161,9 +161,9 @@ where
 
 #[cfg(test)]
 mod tests {
-    use cgmath::{Basis2, Decomposed, Point2, Rad, Vector2};
-    use cgmath::prelude::*;
     use approx::assert_ulps_eq;
+    use cgmath::prelude::*;
+    use cgmath::{Basis2, Decomposed, Point2, Rad, Vector2};
 
     use super::*;
     use crate::primitive::Circle;
diff --git a/src/primitive/polygon.rs b/src/primitive/polygon.rs
index bf536e0..6e5f254 100644
--- a/src/primitive/polygon.rs
+++ b/src/primitive/polygon.rs
@@ -1,16 +1,19 @@
 //! Convex polygon primitive
 
-use cgmath::{BaseFloat, Point2, Vector2};
 use cgmath::prelude::*;
+use cgmath::{vec2, BaseFloat, Point2, Vector2};
+use std::cmp::Ordering;
 
-use crate::{Aabb2, Line2, Ray2};
 use crate::prelude::*;
 use crate::primitive::util::{get_bound, get_max_point};
+use crate::{Aabb2, Line2, Ray2};
 
 /// Convex polygon primitive.
 ///
 /// Can contain any number of vertices, but a high number of vertices will
 /// affect performance of course. Vertices need to be in CCW order.
+/// If using `GJKLeft2`, two vertices should not be in the same place,
+/// and all the vertices should be vertices of the resulting convex polygon.
 #[derive(Debug, Clone, PartialEq)]
 #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
 pub struct ConvexPolygon<S> {
@@ -20,6 +23,8 @@ pub struct ConvexPolygon<S> {
 
 impl<S> ConvexPolygon<S> {
     /// Create a new convex polygon from the given vertices. Vertices need to be in CCW order.
+    /// If using `GJKLeft2`, two vertices should not be in the same place.
+    /// and all the vertices should be vertices of the resulting convex polygon.
     pub fn new(vertices: Vec<Point2<S>>) -> Self {
         Self { vertices }
     }
@@ -41,6 +46,27 @@ where
             support_point(&self.vertices, direction, transform)
         }
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        let v0_iter = self.vertices.iter();
+        let v1_iter = self.vertices.iter().chain(self.vertices.iter()).skip(1);
+        let edge_normals = v0_iter
+            .zip(v1_iter)
+            .map(|(v0, v1)| v1 - v0)
+            .map(|v| vec2(v.y, -v.x).normalize());
+
+        // There better be at least 2 vertices in this convex polygon!
+        edge_normals
+            .max_by(|u, v| {
+                u.dot(*normal)
+                    .partial_cmp(&v.dot(*normal))
+                    .unwrap_or(Ordering::Equal)
+            })
+            .unwrap()
+    }
 }
 
 impl<S> ComputeBound<Aabb2<S>> for ConvexPolygon<S>
@@ -190,13 +216,14 @@ where
         vertices[vertices.len() - 1]
     } else {
         vertices[index_u]
-    }.dot(*direction)
+    }
+    .dot(*direction)
 }
 
 #[cfg(test)]
 mod tests {
-    use cgmath::{Basis2, Decomposed, Point2, Rad, Vector2};
     use approx::assert_ulps_eq;
+    use cgmath::{Basis2, Decomposed, Point2, Rad, Vector2};
 
     use super::*;
     use {Aabb2, Ray2};
@@ -264,6 +291,25 @@ mod tests {
         );
     }
 
+    #[test]
+    fn test_closest_valid_normal() {
+        let vertices = vec![
+            Point2::new(0., 0.),
+            Point2::new(-1., 1.),
+            Point2::new(-1., -1.),
+            Point2::new(1., -1.),
+        ];
+        let polygon = ConvexPolygon::new(vertices);
+        assert_ulps_eq!(
+            vec2(0.5f64.sqrt(), 0.5f64.sqrt()),
+            polygon.closest_valid_normal_local(&vec2(0.75f64.sqrt(), 0.5))
+        );
+        assert_eq!(
+            vec2(-1., 0.),
+            polygon.closest_valid_normal_local(&vec2(-1., 0.))
+        );
+    }
+
     #[test]
     fn test_ray_discrete() {
         let vertices = vec![
diff --git a/src/primitive/polyhedron.rs b/src/primitive/polyhedron.rs
index f62c70a..09377bb 100644
--- a/src/primitive/polyhedron.rs
+++ b/src/primitive/polyhedron.rs
@@ -2,13 +2,13 @@ use std::cmp::Ordering;
 use std::collections::HashMap;
 
 use bit_set::BitSet;
-use cgmath::{BaseFloat, Point3, Vector3};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point3, Vector3};
 
-use crate::{Aabb3, Plane, Ray3};
 use crate::prelude::*;
 use crate::primitive::util::barycentric_point;
 use crate::volume::Sphere;
+use crate::{Aabb3, Plane, Ray3};
 
 #[derive(Debug, Clone, PartialEq)]
 #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
@@ -141,7 +141,8 @@ where
 
     #[inline]
     fn brute_force_support_point(&self, direction: Vector3<S>) -> Point3<S> {
-        let (p, _) = self.vertices
+        let (p, _) = self
+            .vertices
             .iter()
             .map(|v| (v.position, v.position.dot(direction)))
             .fold(
@@ -287,7 +288,8 @@ where
                 vertices[a].position,
                 vertices[b].position,
                 vertices[c].position,
-            ).unwrap(),
+            )
+            .unwrap(),
             ready: false,
         };
         let face_index = faces.len();
@@ -370,16 +372,21 @@ where
         T: Transform<Point3<S>>,
     {
         let p = match self.mode {
-            PolyhedronMode::VertexOnly => self.brute_force_support_point(
-                transform.inverse_transform_vector(*direction).unwrap(),
-            ),
+            PolyhedronMode::VertexOnly => self
+                .brute_force_support_point(transform.inverse_transform_vector(*direction).unwrap()),
 
-            PolyhedronMode::HalfEdge => self.hill_climb_support_point(
-                transform.inverse_transform_vector(*direction).unwrap(),
-            ),
+            PolyhedronMode::HalfEdge => self
+                .hill_climb_support_point(transform.inverse_transform_vector(*direction).unwrap()),
         };
         transform.transform_point(p)
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        unimplemented!("closest_valid_normal_local is only implemented for 2D primitives for now")
+    }
 }
 
 impl<S> ComputeBound<Aabb3<S>> for ConvexPolyhedron<S>
@@ -427,7 +434,8 @@ where
             let v0 = f.vertices.0;
             let v1 = f.vertices.1;
             let v2 = f.vertices.2;
-            let p = (self.vertices[v0].position * u) + (self.vertices[v1].position.to_vec() * v)
+            let p = (self.vertices[v0].position * u)
+                + (self.vertices[v1].position.to_vec() * v)
                 + (self.vertices[v2].position.to_vec() * w);
             Some(p)
         })
@@ -575,13 +583,13 @@ where
 #[cfg(test)]
 mod tests {
 
-    use cgmath::{Decomposed, Point3, Quaternion, Rad, Vector3};
-    use cgmath::prelude::*;
     use approx::assert_ulps_eq;
+    use cgmath::prelude::*;
+    use cgmath::{Decomposed, Point3, Quaternion, Rad, Vector3};
 
     use super::ConvexPolyhedron;
-    use crate::{Aabb3, Ray3};
     use crate::prelude::*;
+    use crate::{Aabb3, Ray3};
 
     #[test]
     fn test_polytope_half_edge() {
diff --git a/src/primitive/primitive2.rs b/src/primitive/primitive2.rs
index 1513f38..79bae97 100644
--- a/src/primitive/primitive2.rs
+++ b/src/primitive/primitive2.rs
@@ -1,11 +1,11 @@
 //! Wrapper enum for 2D primitives
 
-use cgmath::{BaseFloat, Point2, Vector2};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point2, Vector2};
 
-use crate::{Aabb2, Line2, Ray2};
 use crate::prelude::*;
 use crate::primitive::{Circle, ConvexPolygon, Particle2, Rectangle, Square};
+use crate::{Aabb2, Line2, Ray2};
 
 /// Wrapper enum for 2D primitives, that also implements the `Primitive` trait, making it easier
 /// to use many different primitives in algorithms.
@@ -97,6 +97,23 @@ where
             Primitive2::ConvexPolygon(ref polygon) => polygon.support_point(direction, transform),
         }
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        match *self {
+            Primitive2::Particle(_) => panic!(concat!(
+                "Particles don't have valid normals. ",
+                "Please don't use GJKLeft2 where the left collider is a particle"
+            )),
+            Primitive2::Line(ref shape) => shape.closest_valid_normal_local(normal),
+            Primitive2::Circle(ref shape) => shape.closest_valid_normal_local(normal),
+            Primitive2::Rectangle(ref shape) => shape.closest_valid_normal_local(normal),
+            Primitive2::Square(ref shape) => shape.closest_valid_normal_local(normal),
+            Primitive2::ConvexPolygon(ref shape) => shape.closest_valid_normal_local(normal),
+        }
+    }
 }
 
 impl<S> DiscreteTransformed<Ray2<S>> for Primitive2<S>
diff --git a/src/primitive/primitive3.rs b/src/primitive/primitive3.rs
index dd86efc..2a52e71 100644
--- a/src/primitive/primitive3.rs
+++ b/src/primitive/primitive3.rs
@@ -1,11 +1,13 @@
 //! Wrapper enum for 3D primitives
 
-use cgmath::{BaseFloat, Point3, Vector3};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point3, Vector3};
 
-use crate::{Aabb3, Ray3};
 use crate::prelude::*;
-use crate::primitive::{Capsule, ConvexPolyhedron, Cube, Cuboid, Cylinder, Particle3, Quad, Sphere};
+use crate::primitive::{
+    Capsule, ConvexPolyhedron, Cube, Cuboid, Cylinder, Particle3, Quad, Sphere,
+};
+use crate::{Aabb3, Ray3};
 
 /// Wrapper enum for 3D primitives, that also implements the `Primitive` trait, making it easier
 /// to use many different primitives in algorithms.
@@ -167,6 +169,13 @@ where
             }
         }
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        unimplemented!("closest_valid_normal_local is only implemented for 2D primitives for now")
+    }
 }
 
 impl<S> DiscreteTransformed<Ray3<S>> for Primitive3<S>
diff --git a/src/primitive/quad.rs b/src/primitive/quad.rs
index 8ed6f89..f025d95 100644
--- a/src/primitive/quad.rs
+++ b/src/primitive/quad.rs
@@ -1,11 +1,11 @@
 //! Rectangular plane primitive
 
-use cgmath::{BaseFloat, Point3, Vector2, Vector3};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point3, Vector2, Vector3};
 
-use crate::{Aabb3, Ray3, Sphere};
 use crate::prelude::*;
 use crate::primitive::util::get_max_point;
+use crate::{Aabb3, Ray3, Sphere};
 
 /// Rectangular plane primitive. Will lie on the xy plane when not transformed.
 ///
@@ -70,6 +70,13 @@ where
     {
         get_max_point(self.corners.iter(), direction, transform)
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        unimplemented!("closest_valid_normal_local is only implemented for 2D primitives for now")
+    }
 }
 
 impl<S> ComputeBound<Aabb3<S>> for Quad<S>
@@ -125,8 +132,8 @@ mod tests {
 
     use super::*;
     use crate::algorithm::minkowski::GJK3;
-    use cgmath::{Decomposed, Quaternion};
     use crate::primitive::Cuboid;
+    use cgmath::{Decomposed, Quaternion};
 
     fn transform(x: f32, y: f32, z: f32) -> Decomposed<Vector3<f32>, Quaternion<f32>> {
         Decomposed {
@@ -143,9 +150,8 @@ mod tests {
         let transform_1 = transform(0., 0., 1.);
         let transform_2 = transform(0., 0., 1.1);
         let gjk = GJK3::new();
-        assert!(
-            gjk.intersect(&quad, &transform_1, &cuboid, &transform_2)
-                .is_some()
-        );
+        assert!(gjk
+            .intersect(&quad, &transform_1, &cuboid, &transform_2)
+            .is_some());
     }
 }
diff --git a/src/primitive/rectangle.rs b/src/primitive/rectangle.rs
index 8c1d9a0..4209cfe 100644
--- a/src/primitive/rectangle.rs
+++ b/src/primitive/rectangle.rs
@@ -1,11 +1,11 @@
 //! Rectangle primitive
 
-use cgmath::{BaseFloat, Point2, Vector2};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point2, Vector2};
 
-use crate::{Aabb2, Ray2};
 use crate::prelude::*;
 use crate::primitive::util::get_max_point;
+use crate::{Aabb2, Ray2};
 
 /// Rectangle primitive.
 ///
@@ -70,6 +70,17 @@ where
     {
         get_max_point(self.corners.iter(), direction, transform)
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        if normal.x.abs() > normal.y.abs() {
+            Vector2::new(normal.x.signum(), Zero::zero())
+        } else {
+            Vector2::new(Zero::zero(), normal.y.signum())
+        }
+    }
 }
 
 impl<S> ComputeBound<Aabb2<S>> for Rectangle<S>
@@ -149,6 +160,13 @@ where
     {
         self.rectangle.support_point(direction, transform)
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        self.rectangle.closest_valid_normal_local(normal)
+    }
 }
 
 impl<S> ComputeBound<Aabb2<S>> for Square<S>
@@ -184,8 +202,8 @@ where
 
 #[cfg(test)]
 mod tests {
-    use cgmath::{Basis2, Decomposed, Point2, Rad, Vector2};
     use approx::assert_ulps_eq;
+    use cgmath::{vec2, Basis2, Decomposed, Point2, Rad, Vector2};
 
     use super::*;
 
@@ -195,6 +213,27 @@ mod tests {
         assert_eq!(bound(-5., -5., 5., 5.), r.compute_bound())
     }
 
+    #[test]
+    fn test_rectangle_closest_valid_normal() {
+        let r = Rectangle::new(1., 1.);
+        assert_eq!(
+            vec2(1., 0.),
+            r.closest_valid_normal_local(&vec2(0.75f64.sqrt(), 0.5))
+        );
+        assert_eq!(
+            vec2(0., 1.),
+            r.closest_valid_normal_local(&vec2(-0.5, 0.75f64.sqrt()))
+        );
+        assert_eq!(
+            vec2(-1., 0.),
+            r.closest_valid_normal_local(&vec2(-0.75f64.sqrt(), -0.5))
+        );
+        assert_eq!(
+            vec2(0., -1.),
+            r.closest_valid_normal_local(&vec2(0.5, -0.75f64.sqrt()))
+        );
+    }
+
     #[test]
     fn test_rectangle_ray_discrete() {
         let rectangle = Rectangle::new(10., 10.);
diff --git a/src/primitive/sphere.rs b/src/primitive/sphere.rs
index f486ea6..1bd3edf 100644
--- a/src/primitive/sphere.rs
+++ b/src/primitive/sphere.rs
@@ -1,8 +1,8 @@
-use cgmath::{BaseFloat, Point3, Vector3};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point3, Vector3};
 
-use crate::{Aabb3, Ray3};
 use crate::prelude::*;
+use crate::{Aabb3, Ray3};
 
 /// Sphere primitive
 #[derive(Debug, Clone, PartialEq)]
@@ -32,6 +32,13 @@ where
         let direction = transform.inverse_transform_vector(*direction).unwrap();
         transform.transform_point(Point3::from_vec(direction.normalize_to(self.radius)))
     }
+
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff {
+        unimplemented!("closest_valid_normal_local is only implemented for 2D primitives for now")
+    }
 }
 
 impl<S> ComputeBound<Aabb3<S>> for Sphere<S>
@@ -101,8 +108,8 @@ where
 mod tests {
     use std;
 
-    use cgmath::{Decomposed, Point3, Quaternion, Rad, Rotation3, Vector3};
     use approx::assert_ulps_eq;
+    use cgmath::{Decomposed, Point3, Quaternion, Rad, Rotation3, Vector3};
 
     use super::*;
 
diff --git a/src/primitive/util.rs b/src/primitive/util.rs
index 9df4474..b7fb055 100644
--- a/src/primitive/util.rs
+++ b/src/primitive/util.rs
@@ -4,9 +4,9 @@
 use std::ops::Neg;
 
 use crate::{Aabb, Ray3};
-use cgmath::{BaseFloat, BaseNum, Vector2};
-use cgmath::prelude::*;
 use cgmath::num_traits::Float;
+use cgmath::prelude::*;
+use cgmath::{BaseFloat, BaseNum, Vector2};
 
 pub(crate) fn get_max_point<'a, P: 'a, T, I>(vertices: I, direction: &P::Diff, transform: &T) -> P
 where
@@ -37,6 +37,17 @@ where
     vertices.fold(A::zero(), |bound, p| bound.grow(*p))
 }
 
+/// Gets normal of counterclockwise vector by rotating it 90 degrees.
+/// The result is the same length as the vector.
+#[allow(dead_code)]
+#[inline]
+pub(crate) fn vector_normal<S>(a: &Vector2<S>) -> Vector2<S>
+where
+    S: BaseNum,
+{
+    Vector2::new(a.y, S::zero() - a.x)
+}
+
 #[allow(dead_code)]
 #[inline]
 pub(crate) fn triple_product<S>(a: &Vector2<S>, b: &Vector2<S>, c: &Vector2<S>) -> Vector2<S>
@@ -138,8 +149,8 @@ where
 mod tests {
     use std;
 
-    use cgmath::{Basis2, Decomposed, Point2, Rad, Rotation2, Vector2};
     use approx::assert_ulps_eq;
+    use cgmath::{Basis2, Decomposed, Point2, Rad, Rotation2, Vector2};
 
     use super::*;
     use crate::Aabb2;
diff --git a/src/ray.rs b/src/ray.rs
index 5677df5..9356b11 100644
--- a/src/ray.rs
+++ b/src/ray.rs
@@ -2,10 +2,10 @@
 
 use std::marker::PhantomData;
 
+use cgmath::prelude::*;
 use cgmath::{BaseFloat, BaseNum};
 use cgmath::{Point2, Point3};
 use cgmath::{Vector2, Vector3};
-use cgmath::prelude::*;
 
 use crate::traits::{Continuous, ContinuousTransformed, Discrete, DiscreteTransformed};
 
diff --git a/src/traits.rs b/src/traits.rs
index 0534d33..c55ba34 100644
--- a/src/traits.rs
+++ b/src/traits.rs
@@ -1,5 +1,5 @@
-use cgmath::BaseNum;
 use cgmath::prelude::*;
+use cgmath::BaseNum;
 
 /// An intersection test with a result.
 ///
@@ -112,6 +112,22 @@ pub trait Primitive {
     ) -> Self::Point
     where
         T: Transform<Self::Point>;
+
+    /// Get the closest valid normal on the shape to a given normal, in local space.
+    /// A normal is valid if there exists a tangent to a differentiable point on the
+    /// shape that the normal is perpendicular to, pointing outside the shape.
+    ///
+    /// ## Parameters
+    ///
+    /// - `normal`: The normal to use. This is assumed to be normalized.
+    ///
+    /// ## Returns
+    ///
+    /// Returns the closest valid normal to the given normal.
+    fn closest_valid_normal_local(
+        &self,
+        normal: &<Self::Point as EuclideanSpace>::Diff,
+    ) -> <Self::Point as EuclideanSpace>::Diff;
 }
 
 /// Discrete intersection test on transformed primitive
@@ -175,8 +191,8 @@ pub trait TranslationInterpolate<S> {
 
 mod interpolate {
     use super::{Interpolate, TranslationInterpolate};
-    use cgmath::{BaseFloat, Basis2, Basis3, Decomposed, Quaternion, Rad};
     use cgmath::prelude::*;
+    use cgmath::{BaseFloat, Basis2, Basis3, Decomposed, Quaternion, Rad};
 
     impl<S> Interpolate<S> for Quaternion<S>
     where
diff --git a/src/volume/aabb/aabb2.rs b/src/volume/aabb/aabb2.rs
index 7668dc8..45a47b4 100644
--- a/src/volume/aabb/aabb2.rs
+++ b/src/volume/aabb/aabb2.rs
@@ -3,12 +3,12 @@
 
 use std::fmt;
 
-use cgmath::{BaseFloat, BaseNum, Point2, Vector2};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, BaseNum, Point2, Vector2};
 
 use super::{max, min};
-use crate::{Line2, Ray2};
 use crate::prelude::*;
+use crate::{Line2, Ray2};
 
 /// A two-dimensional AABB, aka a rectangle.
 #[derive(Copy, Clone, PartialEq)]
@@ -103,7 +103,9 @@ impl<S: BaseNum> Contains<Aabb2<S>> for Aabb2<S> {
         let other_min = other.min();
         let other_max = other.max();
 
-        other_min.x >= self.min.x && other_min.y >= self.min.y && other_max.x <= self.max.x
+        other_min.x >= self.min.x
+            && other_min.y >= self.min.y
+            && other_max.x <= self.max.x
             && other_max.y <= self.max.y
     }
 }
diff --git a/src/volume/aabb/aabb3.rs b/src/volume/aabb/aabb3.rs
index cb6b1db..04fd3f9 100644
--- a/src/volume/aabb/aabb3.rs
+++ b/src/volume/aabb/aabb3.rs
@@ -3,12 +3,12 @@
 
 use std::fmt;
 
-use cgmath::{BaseFloat, BaseNum, Point3, Vector3};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, BaseNum, Point3, Vector3};
 
 use super::{max, min};
-use crate::{Line3, Plane, Ray3, Sphere};
 use crate::prelude::*;
+use crate::{Line3, Plane, Ray3, Sphere};
 
 /// A three-dimensional AABB, aka a rectangular prism.
 #[derive(Copy, Clone, PartialEq)]
@@ -105,8 +105,12 @@ impl<S: BaseNum> fmt::Debug for Aabb3<S> {
 impl<S: BaseNum> Contains<Point3<S>> for Aabb3<S> {
     #[inline]
     fn contains(&self, p: &Point3<S>) -> bool {
-        self.min.x <= p.x && p.x < self.max.x && self.min.y <= p.y && p.y < self.max.y
-            && self.min.z <= p.z && p.z < self.max.z
+        self.min.x <= p.x
+            && p.x < self.max.x
+            && self.min.y <= p.y
+            && p.y < self.max.y
+            && self.min.z <= p.z
+            && p.z < self.max.z
     }
 }
 
@@ -116,8 +120,11 @@ impl<S: BaseNum> Contains<Aabb3<S>> for Aabb3<S> {
         let other_min = other.min();
         let other_max = other.max();
 
-        other_min.x >= self.min.x && other_min.y >= self.min.y && other_min.z >= self.min.z
-            && other_max.x <= self.max.x && other_max.y <= self.max.y
+        other_min.x >= self.min.x
+            && other_min.y >= self.min.y
+            && other_min.z >= self.min.z
+            && other_max.x <= self.max.x
+            && other_max.y <= self.max.y
             && other_max.z <= self.max.z
     }
 }
@@ -158,7 +165,8 @@ impl<S: BaseFloat> Union<Sphere<S>> for Aabb3<S> {
             sphere.center.x - sphere.radius,
             sphere.center.y - sphere.radius,
             sphere.center.z - sphere.radius,
-        )).grow(sphere.center + Vector3::from_value(sphere.radius))
+        ))
+        .grow(sphere.center + Vector3::from_value(sphere.radius))
     }
 }
 
diff --git a/src/volume/aabb/mod.rs b/src/volume/aabb/mod.rs
index d04f08e..4194b1b 100644
--- a/src/volume/aabb/mod.rs
+++ b/src/volume/aabb/mod.rs
@@ -10,8 +10,8 @@ pub use self::aabb3::Aabb3;
 
 use std::cmp::{Ordering, PartialOrd};
 
-use cgmath::{BaseNum, Point2, Point3};
 use cgmath::prelude::*;
+use cgmath::{BaseNum, Point2, Point3};
 
 use crate::traits::Bound;
 
diff --git a/src/volume/mod.rs b/src/volume/mod.rs
index a0de5f6..ac9fa19 100644
--- a/src/volume/mod.rs
+++ b/src/volume/mod.rs
@@ -3,7 +3,7 @@ pub use self::cylinder::Cylinder;
 pub use self::obb::*;
 pub use self::sphere::*;
 
-mod cylinder;
 mod aabb;
+mod cylinder;
 mod obb;
 mod sphere;
diff --git a/src/volume/sphere.rs b/src/volume/sphere.rs
index 2e56db4..55eeb30 100644
--- a/src/volume/sphere.rs
+++ b/src/volume/sphere.rs
@@ -1,10 +1,10 @@
 //! Bounding sphere
 
-use cgmath::{BaseFloat, Point3, Vector3};
 use cgmath::prelude::*;
+use cgmath::{BaseFloat, Point3, Vector3};
 
-use crate::{Aabb3, Line3, Plane, Ray3};
 use crate::prelude::*;
+use crate::{Aabb3, Line3, Plane, Ray3};
 
 /// Bounding sphere.
 #[derive(Copy, Clone, PartialEq, Debug)]
diff --git a/tests/aabb.rs b/tests/aabb.rs
index 000b82b..3f2649d 100644
--- a/tests/aabb.rs
+++ b/tests/aabb.rs
@@ -1,9 +1,9 @@
 extern crate cgmath;
 extern crate collision;
 
+use cgmath::InnerSpace;
 use cgmath::{Point2, Point3};
 use cgmath::{Vector2, Vector3};
-use cgmath::InnerSpace;
 use collision::{Aabb, Aabb2, Aabb3};
 use collision::{Contains, Continuous, Discrete, SurfaceArea, Union};
 use collision::{Line2, Line3, Ray2, Ray3, Sphere};
@@ -224,7 +224,8 @@ fn test_corners() {
     let corners = Aabb3::new(
         Point3::new(-20isize, 30isize, 5isize),
         Point3::new(10isize, -10isize, -5isize),
-    ).to_corners();
+    )
+    .to_corners();
     assert!(corners.contains(&Point3::new(-20isize, 30isize, -5isize)));
     assert!(corners.contains(&Point3::new(10isize, 30isize, 5isize)));
     assert!(corners.contains(&Point3::new(10isize, -10isize, 5isize)));
diff --git a/tests/dbvt.rs b/tests/dbvt.rs
index 4ddcf72..32b6aab 100644
--- a/tests/dbvt.rs
+++ b/tests/dbvt.rs
@@ -2,11 +2,11 @@ extern crate cgmath;
 extern crate collision;
 extern crate rand;
 
-use cgmath::{Deg, PerspectiveFov, Point2, Point3, Vector2, Vector3};
 use cgmath::prelude::*;
-use collision::{Aabb2, Aabb3, Frustum, Projection, Ray2, Relation};
+use cgmath::{Deg, PerspectiveFov, Point2, Point3, Vector2, Vector3};
 use collision::dbvt::*;
 use collision::prelude::*;
+use collision::{Aabb2, Aabb3, Frustum, Projection, Ray2, Relation};
 use rand::Rng;
 
 #[derive(Debug, Clone)]
diff --git a/tests/frustum.rs b/tests/frustum.rs
index 7081500..a7cdf7b 100644
--- a/tests/frustum.rs
+++ b/tests/frustum.rs
@@ -11,7 +11,8 @@ fn test_contains() {
         aspect: 1f32,
         near: 1f32,
         far: 10f32,
-    }.to_frustum();
+    }
+    .to_frustum();
     assert_eq!(
         frustum.contains(&Sphere {
             center: Point3::new(0f32, 0f32, -5f32),
diff --git a/tests/plane.rs b/tests/plane.rs
index 111e85d..f120e70 100644
--- a/tests/plane.rs
+++ b/tests/plane.rs
@@ -41,7 +41,8 @@ fn test_ray_intersection() {
         Point3::new(5f64, 0f64, 5f64),
         Point3::new(5f64, 5f64, 5f64),
         Point3::new(5f64, 0f64, -1f64),
-    ).unwrap();
+    )
+    .unwrap();
     let r1: Ray3<f64> = Ray::new(
         Point3::new(0f64, 0f64, 0f64),
         Vector3::new(-1f64, 0f64, 0f64).normalize(),