diff --git a/R/geom_pointdensity.R b/R/geom_pointdensity.R
index de2d7ac..4a0ae65 100644
--- a/R/geom_pointdensity.R
+++ b/R/geom_pointdensity.R
@@ -138,8 +138,8 @@ StatPointdensity <- ggproto("StatPointdensity", Stat,
                                 ix <- findInterval(data$x, dens$x)
                                 iy <- findInterval(data$y, dens$y)
                                 ii <- cbind(ix, iy)
-                                data$density <- dens$z[ii]
-
+                                data$density[finites] <- dens$z[ii]
+                                data$density[!finites] <- min(dens$z)
                               } else {
 
                                 if (is.character(method)) {
diff --git a/src/count_neighbors.c b/src/count_neighbors.c
index c32bd3f..abe90b3 100644
--- a/src/count_neighbors.c
+++ b/src/count_neighbors.c
@@ -22,11 +22,18 @@ SEXP count_neighbors_( SEXP x, SEXP y, SEXP r2, SEXP xy ) {
       int s = 0;
       double xi = xp[i];
       double yi = yp[i];
-      for( int j = 0; j < l; j++ ) {
-        double dx = xi - xp[j];
-        double dy = yi - yp[j];
-        if( yxp*dx*dx + xyp*dy*dy <= r2p )
-          s++;
+      for( int j = 0; j < l; j++ ) {\
+        double xj = xp[j];
+        double yj = yp[j];
+        double dx = xi - xj;
+        double dy = yi - yj;
+        if((xi == xj && yi == yj) ||
+           (xi == xj && xyp * dy * dy <= r2p) ||
+           (yi == yj && xyp * dx * dx <= r2p) ||
+           (yxp*dx*dx + xyp*dy*dy <= r2p) ){
+           s++;
+        }
+
       }
       resp[i] = s;
     }