diff --git a/cogl/cogl-journal.c b/cogl/cogl-journal.c
index 7370dd6b6..ea152f31e 100644
--- a/cogl/cogl-journal.c
+++ b/cogl/cogl-journal.c
@@ -1699,7 +1699,7 @@ try_checking_point_hits_entry_after_clipping (CoglJournalEntry *entry,
       software_clip_entry (entry, vertices, &clip_bounds);
       entry_to_screen_polygon (entry, vertices, poly);
 
-      *hit = _cogl_util_point_in_poly (x, y, poly, sizeof (float) * 4, 4);
+      *hit = _cogl_util_point_in_screen_poly (x, y, poly, sizeof (float) * 4, 4);
       return TRUE;
     }
 
@@ -1753,7 +1753,7 @@ _cogl_journal_try_read_pixel (CoglJournal *journal,
 
       entry_to_screen_polygon (entry, vertices, poly);
 
-      if (!_cogl_util_point_in_poly (x, y, poly, sizeof (float) * 4, 4))
+      if (!_cogl_util_point_in_screen_poly (x, y, poly, sizeof (float) * 4, 4))
         continue;
 
       /* FIXME: the journal should have a back pointer to the
diff --git a/cogl/cogl-point-in-poly-private.h b/cogl/cogl-point-in-poly-private.h
index 0cebdf8b5..6b534b6d5 100644
--- a/cogl/cogl-point-in-poly-private.h
+++ b/cogl/cogl-point-in-poly-private.h
@@ -28,11 +28,11 @@
 G_BEGIN_DECLS
 
 int
-_cogl_util_point_in_poly (float point_x,
-                          float point_y,
-                          void *vertices,
-                          size_t stride,
-                          int n_vertices);
+_cogl_util_point_in_screen_poly (float point_x,
+                                 float point_y,
+                                 void *vertices,
+                                 size_t stride,
+                                 int n_vertices);
 
 G_END_DECLS
 
diff --git a/cogl/cogl-point-in-poly.c b/cogl/cogl-point-in-poly.c
index fc38d1d67..87e090eca 100644
--- a/cogl/cogl-point-in-poly.c
+++ b/cogl/cogl-point-in-poly.c
@@ -38,14 +38,40 @@
 #include "config.h"
 #endif
 
+#include "cogl-util.h"
+
 #include <glib.h>
 
+/* We've made a notable change to the original algorithm referenced
+ * above to make sure we have reliable results for screen aligned
+ * rectangles even though there may be some numerical in-precision in
+ * how the vertices of the polygon were calculated.
+ *
+ * We've avoided introducing an epsilon factor to the comparisons
+ * since we feel there's a risk of changing some semantics in ways that
+ * might not be desirable. One of those is that if you transform two
+ * polygons which share an edge and test a point close to that edge
+ * then this algorithm will currently give a positive result for only
+ * one polygon.
+ *
+ * Another concern is the way this algorithm resolves the corner case
+ * where the horizontal ray being cast to count edge crossings may
+ * cross directly through a vertex. The solution is based on the "idea
+ * of Simulation of Simplicity" and "pretends to shift the ray
+ * infinitesimally down so that it either clearly intersects, or
+ * clearly doesn't touch". I'm not familiar with the idea myself so I
+ * expect a misplaced epsilon is likely to break that aspect of the
+ * algorithm.
+ *
+ * The simple solution we've gone for is to pixel align the polygon
+ * vertices which should eradicate most noise due to in-precision.
+ */
 int
-_cogl_util_point_in_poly (float point_x,
-                          float point_y,
-                          void *vertices,
-                          size_t stride,
-                          int n_vertices)
+_cogl_util_point_in_screen_poly (float point_x,
+                                 float point_y,
+                                 void *vertices,
+                                 size_t stride,
+                                 int n_vertices)
 {
   int i, j, c = 0;
 
@@ -58,6 +84,11 @@ _cogl_util_point_in_poly (float point_x,
       float vert_yj = *(float *)((guint8 *)vertices + j * stride +
                                  sizeof (float));
 
+      vert_xi = COGL_UTIL_NEARBYINT (vert_xi);
+      vert_xj = COGL_UTIL_NEARBYINT (vert_xj);
+      vert_yi = COGL_UTIL_NEARBYINT (vert_yi);
+      vert_yj = COGL_UTIL_NEARBYINT (vert_yj);
+
       if (((vert_yi > point_y) != (vert_yj > point_y)) &&
            (point_x < (vert_xj - vert_xi) * (point_y - vert_yi) /
             (vert_yj - vert_yi) + vert_xi) )