diff --git a/skimage/filter/tests/test_thresholding.py b/skimage/filter/tests/test_thresholding.py
index f0a68f7d..8464cb21 100644
--- a/skimage/filter/tests/test_thresholding.py
+++ b/skimage/filter/tests/test_thresholding.py
@@ -67,6 +67,12 @@ class TestSimpleImage():
     def test_isodata_linspace(self):
         assert -63.8 < threshold_isodata(np.linspace(-127, 0, 256)) < -63.6
 
+    def test_isodata_16bit(self):
+        np.random.seed(0)
+        imfloat = np.random.rand(256, 256)
+        t = threshold_isodata(imfloat, nbins=1024)
+        assert 0.49 < t < 0.51
+
     def test_threshold_adaptive_generic(self):
         def func(arr):
             return arr.sum() / arr.shape[0]
diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py
index 6ccf8f6b..da74fe47 100644
--- a/skimage/filter/thresholding.py
+++ b/skimage/filter/thresholding.py
@@ -243,7 +243,7 @@ def threshold_isodata(image, nbins=256):
     cpmfl = np.cumsum(pmf, dtype=np.float32)
     cpmfh = np.cumsum(pmf[::-1], dtype=np.float32)[::-1]
 
-    binnums = np.arange(pmf.size, dtype=np.uint8)
+    binnums = np.arange(pmf.size, dtype=np.min_scalar_type(nbins))
     # l and h contain average value of pixels in sum of bins, calculated
     # from lower to higher and from higher to lower respectively.
     l = np.ma.divide(np.cumsum(pmf * binnums, dtype=np.float32), cpmfl)