diff --git a/src/sample_point_cloud.cpp b/src/sample_point_cloud.cpp
index d7a7398..9a99561 100644
--- a/src/sample_point_cloud.cpp
+++ b/src/sample_point_cloud.cpp
@@ -35,7 +35,15 @@ IGL_INLINE void blue_noise_downsample(const Eigen::MatrixBase<DerivedV> & X,
   
     // Make a uniform random sampling with 30*expected_number_of_points.
     const int nx = X.rows();
-  
+
+    // An empty input has nothing to sample. Return early before computing the
+    // grid bounds, since colwise().minCoeff()/maxCoeff() are undefined on a
+    // zero-row matrix and would otherwise crash.
+    if (nx == 0) {
+        XI.resize(0);
+        return;
+    }
+
     // Rescale so that s = 1
     Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> Xs = ((X.rowwise()-X.colwise().minCoeff())/s).template cast<int>();
     const int w = Xs.maxCoeff() + 1;
diff --git a/tests/test_examples.py b/tests/test_examples.py
index d73ba5c..7a559bc 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -81,6 +81,14 @@ def test_mesh_sampling(self):
         if bc1.shape == bc3.shape:
             self.assertFalse(np.all(bc1 == bc3))
 
+    def test_poisson_disk_downsample_empty(self):
+        import point_cloud_utils as pcu
+        import numpy as np
+
+        empty = np.zeros((0, 3), dtype=np.float64)
+        s_idx = pcu.downsample_point_cloud_poisson_disk(empty, 0.1, random_seed=1234567)
+        self.assertEqual(s_idx.shape[0], 0)
+
     def test_downsample_point_cloud_on_voxel_grid(self):
         import point_cloud_utils as pcu
         import numpy as np