diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx
index 67f425f7..8d49b377 100644
--- a/skimage/feature/_texture.pyx
+++ b/skimage/feature/_texture.pyx
@@ -129,6 +129,9 @@ def _local_binary_pattern(double[:, ::1] image,
     cdef Py_ssize_t rot_index, n_ones
     cdef cnp.int8_t first_zero, first_one
 
+    # To compute the variance features
+    cdef double sum_, var_, texture_i
+
     for r in range(image.shape[0]):
         for c in range(image.shape[1]):
             for i in range(P):
@@ -144,8 +147,24 @@ def _local_binary_pattern(double[:, ::1] image,
 
             lbp = 0
 
-            # if method == 'uniform' or method == 'var':
-            if method == 'U' or method == 'N' or method == 'V':
+            # if method == 'var':
+            if method == 'V':
+                # Compute the variance without passing from numpy.
+                # Following the LBP paper, we're taking a biased estimate
+                # of the variance (ddof=0)
+                sum_ = 0.0
+                var_ = 0.0
+                for i in range(P):
+                    texture_i = texture[i]
+                    sum_ += texture_i
+                    var_ += texture_i * texture_i
+                var_ = (var_ - (sum_ * sum_) / P) / P
+                if var_ != 0:
+                    lbp = var_
+                else:
+                    lbp = np.nan
+            # if method == 'uniform':
+            elif method == 'U' or method == 'N':
                 # determine number of 0 - 1 changes
                 changes = 0
                 for i in range(P - 1):
@@ -186,7 +205,7 @@ def _local_binary_pattern(double[:, ::1] image,
 
                     if changes <= 2:
                         # We have a uniform pattern
-                        n_ones = 0  # determies the number of ones
+                        n_ones = 0  # determines the number of ones
                         first_one = -1  # position was the first one
                         first_zero = -1  # position of the first zero
                         for i in range(P):
@@ -215,13 +234,6 @@ def _local_binary_pattern(double[:, ::1] image,
                             lbp += signed_texture[i]
                     else:
                         lbp = P + 1
-
-                    if method == 'V':
-                        var = np.var(texture)
-                        if var != 0:
-                            lbp /= var
-                        else:
-                            lbp = np.nan
             else:
                 # method == 'default'
                 for i in range(P):
diff --git a/skimage/feature/tests/test_texture.py b/skimage/feature/tests/test_texture.py
index e4fb6acb..56e99b95 100644
--- a/skimage/feature/tests/test_texture.py
+++ b/skimage/feature/tests/test_texture.py
@@ -183,21 +183,27 @@ class TestLBP():
         np.testing.assert_array_equal(lbp, ref)
 
     def test_var(self):
-        lbp = local_binary_pattern(self.image, 8, 1, 'var')
-        ref = np.array([[0.        , 0.00072786, 0.        , 0.00115377,
-                         0.00032355, 0.00224467],
-                        [0.00051758, 0.        , 0.0026383 , 0.00163246,
-                         0.00027414, 0.00041124],
-                        [0.00192834, 0.00130368, 0.00042095, 0.00171894,
-                         0.        , 0.00063726],
-                        [0.00023048, 0.00019464 , 0.00082291, 0.00225386,
-                         0.00076696, 0.        ],
-                        [0.00097253, 0.00013236, 0.0009134 , 0.0014467 ,
-                         0.        , 0.00082472],
-                        [0.00024701, 0.0012277 , 0.        , 0.00109869,
-                         0.00015445, 0.00035881]])
-        np.testing.assert_array_almost_equal(lbp, ref)
+        # Test idea: mean of variance is estimate of overall variance.
 
+        # Fix random seed for test stability.
+        np.random.seed(13141516)
+
+        # Create random image with known variance.
+        image = np.random.random((500, 500))
+        target_std = 0.3
+        image = image / image.std() * target_std
+
+        # Use P=4 to avoid interpolation effects
+        P, R = 4, 1
+        lbp = local_binary_pattern(image, P, R, 'var')
+
+        # Take central part to avoid border effect.
+        lbp = lbp[5:-5,5:-5]
+
+        # The LBP variance is biased (ddof=0), correct for that.
+        expected = target_std**2 * (P-1)/P
+
+        np.testing.assert_almost_equal(lbp.mean(), expected, 4)
 
     def test_nri_uniform(self):
         lbp = local_binary_pattern(self.image, 8, 1, 'nri_uniform')