diff --git a/src/geos_functions.jl b/src/geos_functions.jl
index 05a7705..9655cc1 100644
--- a/src/geos_functions.jl
+++ b/src/geos_functions.jl
@@ -1481,6 +1481,16 @@ function nearestPoints(obj1::Geometry, obj2::Geometry, context::GEOSContext = ge
     end
 end
 
+function nearestPoints(obj1::PreparedGeometry, obj2::Geometry, context::GEOSContext = get_context(obj1))
+    points = GEOSPreparedNearestPoints_r(context, obj1, obj2)
+    if points == C_NULL
+        return Point[]
+    else
+        return Point[Point(getCoordinates(points, 1, context), context),
+                     Point(getCoordinates(points, 2, context), context)]
+    end
+end
+
 # -----
 # Precision functions
 # -----
diff --git a/test/test_geos_functions.jl b/test/test_geos_functions.jl
index d94e4aa..7da3d89 100644
--- a/test/test_geos_functions.jl
+++ b/test/test_geos_functions.jl
@@ -431,21 +431,48 @@ end
     # GEOSNearestPointsTest
     geom1_ = LibGEOS.readgeom("POLYGON EMPTY")
     geom2_ = LibGEOS.readgeom("POLYGON EMPTY")
+    prepGeom1_ = LibGEOS.prepareGeom(geom1_)
     @test LibGEOS.nearestPoints(geom1_, geom2_) == Point[]
+    @test LibGEOS.nearestPoints(prepGeom1_, geom2_) == Point[]
     LibGEOS.destroyGeom(geom1_)
     LibGEOS.destroyGeom(geom2_)
+    LibGEOS.destroyPreparedGeom(prepGeom1_)
 
-    geom1_ = LibGEOS.readgeom("POLYGON((1 1,1 5,5 5,5 1,1 1))")
-    geom2_ = LibGEOS.readgeom("POLYGON((8 8, 9 9, 9 10, 8 8))")
-    coords_ = LibGEOS.nearestPoints(geom1_, geom2_)
-    @test coords_ isa Vector{Point}
-    @test length(coords_) == 2
-    @test GeoInterface.coordinates(coords_[1]) ≈ [5.0, 5.0] atol = 1e-5
-    @test GeoInterface.coordinates(coords_[2]) ≈ [8.0, 8.0] atol = 1e-5
-    LibGEOS.destroyGeom(geom1_)
-    LibGEOS.destroyGeom(geom2_)
-    LibGEOS.destroyGeom(coords_[1])
-    LibGEOS.destroyGeom(coords_[2])
+    for (wkt1_, wkt2_, wkt3_, wkt4_) in [
+        ["POLYGON((1 1,1 5,5 5,5 1,1 1))", "POLYGON((8 8, 9 9, 9 10, 8 8))", "POINT(5 5)", "POINT(8 8)"],
+        ["POLYGON((1 1,1 5,5 5,5 1,1 1))", "POINT(2 2)", "POINT(2 2)", "POINT(2 2)"],
+        ["LINESTRING(1 5,5 5,5 1,1 1)", "POINT(2 2)", "POINT(2 1)", "POINT(2 2)"],
+        ["LINESTRING(0 0,10 10)", "LINESTRING(0 10,10 0)", "POINT(5 5)", "POINT(5 5)"],
+        ["POLYGON((0 0,10 0,10 10,0 10,0 0))", "LINESTRING(8 5,12 5)", "POINT(8 5)", "POINT(8 5)"]
+        ]
+        geom1_ = LibGEOS.readgeom(wkt1_)
+        geom2_ = LibGEOS.readgeom(wkt2_)
+        geom3_ = LibGEOS.readgeom(wkt3_)
+        geom4_ = LibGEOS.readgeom(wkt4_)
+        prepGeom1_ = LibGEOS.prepareGeom(geom1_)
+
+        coords_ = LibGEOS.nearestPoints(geom1_, geom2_)
+        @test coords_ isa Vector{Point}
+        @test length(coords_) == 2
+        @test LibGEOS.equals(coords_[1], geom3_)
+        @test LibGEOS.equals(coords_[2], geom4_)
+
+        prepCoords_ = LibGEOS.nearestPoints(prepGeom1_, geom2_)
+        @test prepCoords_ isa Vector{Point}
+        @test length(prepCoords_) == 2
+        @test LibGEOS.equals(prepCoords_[1], geom3_)
+        @test LibGEOS.equals(prepCoords_[2], geom4_)
+
+        LibGEOS.destroyGeom(geom1_)
+        LibGEOS.destroyGeom(geom2_)
+        LibGEOS.destroyGeom(geom3_)
+        LibGEOS.destroyGeom(geom4_)
+        LibGEOS.destroyGeom(coords_[1])
+        LibGEOS.destroyGeom(coords_[2])
+        LibGEOS.destroyGeom(prepCoords_[1])
+        LibGEOS.destroyGeom(prepCoords_[2])
+        LibGEOS.destroyPreparedGeom(prepGeom1_)
+    end
 
     # GEOSNodeTest
     geom1_ = LibGEOS.readgeom("LINESTRING(0 0, 10 10, 10 0, 0 10)")
diff --git a/test/test_geos_operations.jl b/test/test_geos_operations.jl
index 8d408e3..ea126f8 100644
--- a/test/test_geos_operations.jl
+++ b/test/test_geos_operations.jl
@@ -141,18 +141,6 @@ end
     @test GeoInterface.coordinates(startPoint(g1)) ≈ [0, 0] atol = 1e-5
     @test GeoInterface.coordinates(endPoint(g1)) ≈ [10, 10] atol = 1e-5
 
-    # GEOSNearestPointsTest
-    g1 = readgeom("POLYGON EMPTY")
-    g2 = readgeom("POLYGON EMPTY")
-    @test length(nearestPoints(g1, g2)) == 0
-
-    g1 = readgeom("POLYGON((1 1,1 5,5 5,5 1,1 1))")
-    g2 = readgeom("POLYGON((8 8, 9 9, 9 10, 8 8))")
-    points = nearestPoints(g1, g2)
-    @test length(points) == 2
-    @test GeoInterface.coordinates(points[1])[1:2] == [5.0, 5.0]
-    @test GeoInterface.coordinates(points[2])[1:2] == [8.0, 8.0]
-
     # GEOSNodeTest
     g1 = node(readgeom("LINESTRING(0 0, 10 10, 10 0, 0 10)"))
     LibGEOS.normalize!(g1)