diff --git a/geom.go b/geom.go
index bbd1a2d58..5a2e7f4f3 100644
--- a/geom.go
+++ b/geom.go
@@ -33,6 +33,12 @@ func (g *GeoM) Reset() {
 	g.impl.Reset()
 }
 
+// Apply pre-multiplies a vector (x, y) by the matrix.
+// In other words, Apply calculates GeoM * (x, y).
+func (g *GeoM) Apply(x, y float64) (x2, y2 float64) {
+	return g.impl.Apply(x, y)
+}
+
 // Element returns a value of a matrix at (i, j).
 func (g *GeoM) Element(i, j int) float64 {
 	a, b, c, d, tx, ty := g.impl.Elements()
diff --git a/geom_test.go b/geom_test.go
index b9dd7c443..11ac27938 100644
--- a/geom_test.go
+++ b/geom_test.go
@@ -116,6 +116,67 @@ func TestGeoMConcatSelf(t *testing.T) {
 	}
 }
 
+func TestGeoMApply(t *testing.T) {
+	trans := GeoM{}
+	trans.Translate(1, 2)
+
+	scale := GeoM{}
+	scale.Scale(1.5, 2.5)
+
+	cpx := GeoM{}
+	cpx.Rotate(math.Pi)
+	cpx.Scale(1.5, 2.5)
+	cpx.Translate(-2, -3)
+
+	cases := []struct {
+		GeoM  GeoM
+		InX   float64
+		InY   float64
+		OutX  float64
+		OutY  float64
+		Delta float64
+	}{
+		{
+			GeoM:  GeoM{},
+			InX:   3.14159,
+			InY:   2.81828,
+			OutX:  3.14159,
+			OutY:  2.81828,
+			Delta: 0.00001,
+		},
+		{
+			GeoM:  trans,
+			InX:   3.14159,
+			InY:   2.81828,
+			OutX:  4.14159,
+			OutY:  4.81828,
+			Delta: 0.00001,
+		},
+		{
+			GeoM:  scale,
+			InX:   3.14159,
+			InY:   2.81828,
+			OutX:  4.71239,
+			OutY:  7.04570,
+			Delta: 0.00001,
+		},
+		{
+			GeoM:  cpx,
+			InX:   3.14159,
+			InY:   2.81828,
+			OutX:  -6.71239,
+			OutY:  -10.04570,
+			Delta: 0.00001,
+		},
+	}
+	for _, c := range cases {
+		rx, ry := c.GeoM.Apply(c.InX, c.InY)
+		if math.Abs(rx-c.OutX) > c.Delta || math.Abs(ry-c.OutY) > c.Delta {
+			t.Errorf("%v.Apply(%v, %v) = (%v, %v), want (%v, %v)", c.GeoM, c.InX, c.InY, rx, ry, c.OutX, c.OutY)
+		}
+	}
+}
+
 func BenchmarkGeoM(b *testing.B) {
 	var m GeoM
 	for i := 0; i < b.N; i++ {
diff --git a/internal/affine/geom.go b/internal/affine/geom.go
index d43d1c90f..ce1c4fa47 100644
--- a/internal/affine/geom.go
+++ b/internal/affine/geom.go
@@ -38,6 +38,13 @@ func (g *GeoM) Reset() {
 	g.inited = false
 }
 
+func (g *GeoM) Apply(x, y float64) (x2, y2 float64) {
+	if !g.inited {
+		return x, y
+	}
+	return g.a*x + g.b*y + g.tx, g.c*x + g.d*y + g.ty
+}
+
 func (g *GeoM) Elements() (a, b, c, d, tx, ty float64) {
 	if !g.inited {
 		return 1, 0, 0, 1, 0, 0