diff --git a/geometry.go b/geometry.go
index 72d148e..4680872 100644
--- a/geometry.go
+++ b/geometry.go
@@ -401,13 +401,11 @@ func (m Matrix) Project(u Vec) Vec {
 
 // Unproject does the inverse operation to Project.
 //
-// It turns out that multiplying a vector by the inverse matrix of m can be nearly-accomplished by
-// subtracting the translate part of the matrix and multplying by the inverse of the top-left 2x2
-// matrix, and the inverse of a 2x2 matrix is simple enough to just be inlined in the computation.
-//
 // Time complexity is O(1).
 func (m Matrix) Unproject(u Vec) Vec {
-	d := (m[0] * m[3]) - (m[1] * m[2])
-	u.X, u.Y = (u.X-m[4])/d, (u.Y-m[5])/d
-	return Vec{u.X*m[3] - u.Y*m[1], u.Y*m[0] - u.X*m[2]}
+	det := m[0]*m[3] - m[2]*m[1]
+	return Vec{
+		(m[3]*(u.X-m[4]) - m[2]*(u.Y-m[5])) / det,
+		(-m[1]*(u.X-m[4]) + m[0]*(u.Y-m[5])) / det,
+	}
 }
diff --git a/geometry_test.go b/geometry_test.go
index e1c1a6f..d766b1d 100644
--- a/geometry_test.go
+++ b/geometry_test.go
@@ -2,6 +2,8 @@ package pixel_test
 
 import (
 	"fmt"
+	"github.com/stretchr/testify/assert"
+	"math"
 	"testing"
 
 	"github.com/faiface/pixel"
@@ -77,3 +79,86 @@ func TestResizeRect(t *testing.T) {
 		})
 	}
 }
+
+func TestMatrix_Unproject(t *testing.T) {
+	const delta = 1e-15
+	t.Run("for rotated matrix", func(t *testing.T) {
+		matrix := pixel.IM.
+			Rotated(pixel.ZV, math.Pi/2)
+		unprojected := matrix.Unproject(pixel.V(0, 1))
+		assert.InDelta(t, unprojected.X, 1, delta)
+		assert.InDelta(t, unprojected.Y, 0, delta)
+	})
+	t.Run("for moved matrix", func(t *testing.T) {
+		matrix := pixel.IM.
+			Moved(pixel.V(1, 2))
+		unprojected := matrix.Unproject(pixel.V(2, 5))
+		assert.InDelta(t, unprojected.X, 1, delta)
+		assert.InDelta(t, unprojected.Y, 3, delta)
+	})
+	t.Run("for scaled matrix", func(t *testing.T) {
+		matrix := pixel.IM.
+			Scaled(pixel.ZV, 2)
+		unprojected := matrix.Unproject(pixel.V(2, 4))
+		assert.InDelta(t, unprojected.X, 1, delta)
+		assert.InDelta(t, unprojected.Y, 2, delta)
+	})
+	t.Run("for scaled, rotated and moved matrix", func(t *testing.T) {
+		matrix := pixel.IM.
+			Scaled(pixel.ZV, 2).
+			Rotated(pixel.ZV, math.Pi/2).
+			Moved(pixel.V(2, 2))
+		unprojected := matrix.Unproject(pixel.V(-2, 6))
+		assert.InDelta(t, unprojected.X, 2, delta)
+		assert.InDelta(t, unprojected.Y, 2, delta)
+	})
+	t.Run("for rotated and moved matrix", func(t *testing.T) {
+		matrix := pixel.IM.
+			Rotated(pixel.ZV, math.Pi/2).
+			Moved(pixel.V(1, 1))
+		unprojected := matrix.Unproject(pixel.V(1, 2))
+		assert.InDelta(t, unprojected.X, 1, delta)
+		assert.InDelta(t, unprojected.Y, 0, delta)
+	})
+	t.Run("for projected vertices using all kinds of matrices", func(t *testing.T) {
+		namedMatrices := map[string]pixel.Matrix{
+			"IM":                        pixel.IM,
+			"Scaled":                    pixel.IM.Scaled(pixel.ZV, 0.5),
+			"Scaled x 2":                pixel.IM.Scaled(pixel.ZV, 2),
+			"Rotated":                   pixel.IM.Rotated(pixel.ZV, math.Pi/4),
+			"Moved":                     pixel.IM.Moved(pixel.V(0.5, 1)),
+			"Moved 2":                   pixel.IM.Moved(pixel.V(-1, -0.5)),
+			"Scaled and Rotated":        pixel.IM.Scaled(pixel.ZV, 0.5).Rotated(pixel.ZV, math.Pi/4),
+			"Scaled, Rotated and Moved": pixel.IM.Scaled(pixel.ZV, 0.5).Rotated(pixel.ZV, math.Pi/4).Moved(pixel.V(1, 2)),
+			"Rotated and Moved":         pixel.IM.Rotated(pixel.ZV, math.Pi/4).Moved(pixel.V(1, 2)),
+		}
+		vertices := [...]pixel.Vec{
+			pixel.V(0, 0),
+			pixel.V(5, 0),
+			pixel.V(5, 10),
+			pixel.V(0, 10),
+			pixel.V(-5, 10),
+			pixel.V(-5, 0),
+			pixel.V(-5, -10),
+			pixel.V(0, -10),
+			pixel.V(5, -10),
+		}
+		for matrixName, matrix := range namedMatrices {
+			for _, vertex := range vertices {
+				testCase := fmt.Sprintf("for matrix %s and vertex %v", matrixName, vertex)
+				t.Run(testCase, func(t *testing.T) {
+					projected := matrix.Project(vertex)
+					unprojected := matrix.Unproject(projected)
+					assert.InDelta(t, vertex.X, unprojected.X, delta)
+					assert.InDelta(t, vertex.Y, unprojected.Y, delta)
+				})
+			}
+		}
+	})
+	t.Run("for singular matrix", func(t *testing.T) {
+		matrix := pixel.Matrix{0, 0, 0, 0, 0, 0}
+		unprojected := matrix.Unproject(pixel.ZV)
+		assert.True(t, math.IsNaN(unprojected.X))
+		assert.True(t, math.IsNaN(unprojected.Y))
+	})
+}