diff --git a/doc/examples/plot_multiblock_local_binary_pattern.py b/doc/examples/plot_multiblock_local_binary_pattern.py
new file mode 100644
index 00000000..a05056b4
--- /dev/null
+++ b/doc/examples/plot_multiblock_local_binary_pattern.py
@@ -0,0 +1,87 @@
+"""
+===========================================================
+Multi-Block Local Binary Pattern for texture classification
+===========================================================
+
+In this example, we will see how to compute the multi-block
+local binary pattern at a specified image and how to visualize it.
+
+The features are calculated in a way similar to local binary
+patterns, except that block summed up pixel values
+rather than pixel values are used.
+
+`MB-LBP` is an extension of LBP that can be computed on any
+scales in a constant time using integral image. It consists of
+`9` equal-sized rectangles. They are used to compute a feature.
+Sum of pixels' intensity values in each of them are compared
+to the central rectangle and depending on comparison result,
+the feature descriptor is computed.
+
+We will start with a simple image that we will generate by our
+own to show how the `MB-LBP` works. We will create a `(9, 9)`
+rectangle with and divide it into `9` blocks. After this
+we will apply `MB-LBP` on it.
+
+
+"""
+from __future__ import print_function
+from skimage.feature import multiblock_local_binary_pattern
+import numpy as np
+from skimage.util import img_as_float
+from skimage.transform import integral_image
+
+# Create dummy matrix where first and fifth
+# rectangles have greater value than the central one
+# Therefore, the following bits should be 1.
+test_img = np.zeros((9, 9), dtype='uint8')
+test_img[3:6, 3:6] = 1
+test_img[:3, :3] = 50
+test_img[6:, 6:] = 50
+
+# MB-LBP is filled in reverse order.
+# So the first and fifth bits from the end should
+# be filled.
+correct_answer = 0b10001000
+
+# The function accepts the float images.
+# Also it has to be C-contiguous.
+test_img = img_as_float(test_img)
+int_img = integral_image(test_img)
+
+lbp_code = multiblock_local_binary_pattern(int_img, 0, 0, 3, 3)
+
+print(lbp_code == correct_answer)
+
+"""
+Now let's apply the operator to a real image and see how the visualization works.
+"""
+
+from __future__ import print_function
+from skimage.feature import (multiblock_local_binary_pattern,
+                             visualize_multiblock_lbp)
+from skimage.util import img_as_float
+from skimage.transform import integral_image
+from skimage import data
+from matplotlib import pyplot as plt
+
+test_img = data.coins()
+
+test_img = img_as_float(test_img)
+int_img = integral_image(test_img)
+
+lbp_code = multiblock_local_binary_pattern(int_img, 0, 0, 90, 90)
+
+img = visualize_multiblock_lbp(test_img, 0, 0, 90, 90,
+                               lbp_code=lbp_code)
+
+
+plt.imshow(img, interpolation='nearest')
+
+"""
+.. image:: PLOT2RST.current_figure
+
+On the above plot we see the result of computing a `MB-LBP` and visualization
+of the computed feature. The rectangles that have less intensity than the central
+rectangle are marked with cyan color. The ones that have bigger intensity values
+are marked with white color. The central rectangle is left untouched.
+"""
diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py
index dc25f55d..7cac9954 100644
--- a/skimage/feature/__init__.py
+++ b/skimage/feature/__init__.py
@@ -2,8 +2,10 @@ from ._canny import canny
 from ._daisy import daisy
 from ._hog import hog
 from .texture import (greycomatrix, greycoprops,
-                      local_binary_pattern, multiblock_local_binary_pattern,
+                      local_binary_pattern,
                       visualize_multiblock_lbp)
+
+from ._texture import multiblock_local_binary_pattern
 from .peak import peak_local_max
 from .corner import (corner_kitchen_rosenfeld, corner_harris,
                      corner_shi_tomasi, corner_foerstner, corner_subpix,
diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx
index d716e690..7add3fe0 100644
--- a/skimage/feature/_texture.pyx
+++ b/skimage/feature/_texture.pyx
@@ -264,3 +264,174 @@ def _local_binary_pattern(double[:, ::1] image,
                 output[r, c] = lbp
 
     return np.asarray(output)
+
+
+cdef inline Py_ssize_t _clip(Py_ssize_t x, Py_ssize_t low,
+                             Py_ssize_t high) nogil:
+    """Clips coordinate between high and low.
+
+    Parameters
+    ----------
+    x : int
+        Coordinate to be clipped.
+    low : int
+        The lower bound.
+    high : int
+        The higher bound.
+
+    Returns
+    -------
+    x : int
+        `x` clipped between `high` and `low`.
+    """
+
+    if(x > high):
+        return high
+    if(x < low):
+        return low
+    return x
+
+
+cdef inline cnp.double_t _integ(
+    cnp.double_t[:, ::1] img, Py_ssize_t r0, Py_ssize_t c0,
+        Py_ssize_t r1, Py_ssize_t c1) nogil:
+    """Integrate over the integral image in the given window
+
+    This method was created so that `multiblock_local_binary_pattern`
+    does not have to make a Python call.
+
+    Parameters
+    ----------
+    img : array
+        The integral image over which to integrate.
+    r0 : int
+        The row number of the top left corner.
+    c0 : int
+        The column number of the top left corner.
+    r1 : int
+        The row number of the bottom right corner.
+    c1 : int
+        The column number of the bottom right corner.
+
+    Returns
+    -------
+    ans : double
+        The integral over the given window.
+    """
+
+    r = _clip(r0, 0, img.shape[0] - 1)
+    c = _clip(c0, 0, img.shape[1] - 1)
+
+    r2 = _clip(r1, 0, img.shape[0] - 1)
+    c2 = _clip(c1, 0, img.shape[1] - 1)
+
+    cdef cnp.double_t ans = img[r1, c1]
+
+    if (r0 >= 1) and (c0 >= 1):
+        ans += img[r0 - 1, c0 - 1]
+
+    if (r0 >= 1):
+        ans -= img[r0 - 1, c1]
+
+    if (c0 >= 1):
+        ans -= img[r1, c0 - 1]
+
+    return ans
+
+
+def multiblock_local_binary_pattern(cnp.double_t[:, ::1] int_image,
+                                    Py_ssize_t x,
+                                    Py_ssize_t y,
+                                    Py_ssize_t width,
+                                    Py_ssize_t height):
+    """Multi-block local binary pattern.
+
+    The features are calculated in a way similar to local binary
+    patterns, except that block summed up pixel values
+    rather than pixel values are used.
+
+    MB-LBP is an extension of LBP that can be computed on any
+    scales in a constant time using integral image. It consists of
+    9 equal-sized rectangles. They are used to compute a feature.
+    Sum of pixels' intensity values in each of them are compared
+    to the central rectangle and depending on comparison result,
+    the feature descriptor is computed.
+
+    Parameters
+    ----------
+    int_image : (N, M) double array
+        Integral image.
+    x : int
+        X-coordinate of top left corner of a rectangle containing feature.
+    y : int
+        Y-coordinate of top left corner of a rectangle containing feature.
+    width : int
+        Width of one of 9 equal rectangles that will be used to compute
+        a feature.
+    height : int
+        Height of one of 9 equal rectangles that will be used to compute
+        a feature.
+
+    Returns
+    -------
+    output : int
+        8bit MB-LBP feature descriptor.
+
+    References
+    ----------
+    .. [1] Face Detection Based on Multi-Block LBP
+           Representation. Lun Zhang, Rufeng Chu, Shiming Xiang, Shengcai Liao,
+           Stan Z. Li
+           http://www.cbsr.ia.ac.cn/users/scliao/papers/Zhang-ICB07-MBLBP.pdf
+    """
+
+    # Top-left coordinates of central rectangle
+    cdef:
+        Py_ssize_t central_rect_x = x + width
+        Py_ssize_t central_rect_y = y + height
+
+    # Sum of intensity values of central rectangle
+    cdef double central_rect_val = _integ(int_image,
+                                          central_rect_y,
+                                          central_rect_x,
+                                          central_rect_y + height - 1,
+                                          central_rect_x + width - 1)
+
+    #print central_rect_x, central_rect_y
+
+    # Offsets of neighbour rectangles relative to central one.
+    # It has order starting from top left and going clockwise
+    cdef:
+        Py_ssize_t *x_offsets = [-1, 0, 1, 1, 1, 0, -1, -1]
+        Py_ssize_t *y_offsets = [-1, -1, -1, 0, 1, 1, 1, 0]
+        Py_ssize_t element_num, offset_x, offset_y
+        Py_ssize_t current_rect_x, current_rect_y
+        double current_rect_val
+        int has_greater_value
+        int lbp_code = 0
+
+    for element_num in range(8):
+
+        offset_x = x_offsets[element_num]
+        offset_y = y_offsets[element_num]
+
+
+        current_rect_x = central_rect_x + offset_x * width
+        current_rect_y = central_rect_y + offset_y * height
+
+
+        current_rect_val = _integ(int_image,
+                                  current_rect_y,
+                                  current_rect_x,
+                                  current_rect_y + height - 1,
+                                  current_rect_x + width - 1)
+
+
+        has_greater_value = current_rect_val >= central_rect_val
+
+        # If current rectangle's intensity value is bigger
+        # make corresponding bit to 1.
+        lbp_code |= has_greater_value << (7 - element_num)
+
+    return lbp_code
+
diff --git a/skimage/feature/tests/test_texture.py b/skimage/feature/tests/test_texture.py
index 8635c805..c9e2c7d5 100644
--- a/skimage/feature/tests/test_texture.py
+++ b/skimage/feature/tests/test_texture.py
@@ -6,9 +6,9 @@ from skimage.feature import (
                             multiblock_local_binary_pattern
                             )
 
-
 from skimage._shared.testing import test_parallel
 from skimage.transform import integral_image
+from skimage.util import img_as_float
 
 class TestGLCM():
 
@@ -235,7 +235,10 @@ class TestLBP():
                         [ 9, 58,  0, 57,  7, 14]])
         np.testing.assert_array_almost_equal(lbp, ref)
 
-    def test_multiblock_lbp(self):
+
+class TestMBLBP():
+
+    def test_single_mblbp(self):
 
         # Create dummy matrix where first and fifth
         # rectangles have greater value than the central one
@@ -245,11 +248,19 @@ class TestLBP():
         test_img[:3, :3] = 255
         test_img[6:, 6:] = 255
 
+        # MB-LBP is filled in reverse order.
+        # So the first and fifth bits from the end should
+        # be filled.
+        correct_answer = 0b10001000
+
+        # The function accepts the float images.
+        # Also it has to be C-contiguous.
+        test_img = img_as_float(test_img)
         int_img = integral_image(test_img)
 
         lbp_code = multiblock_local_binary_pattern(int_img, 0, 0, 3, 3)
 
-        np.testing.assert_equal(lbp_code, 17)
+        np.testing.assert_equal(lbp_code, correct_answer)
 
 
 if __name__ == '__main__':
diff --git a/skimage/feature/texture.py b/skimage/feature/texture.py
index 3d82bf44..7d827666 100644
--- a/skimage/feature/texture.py
+++ b/skimage/feature/texture.py
@@ -4,10 +4,9 @@ Methods to characterize image textures.
 
 import numpy as np
 from .._shared.utils import assert_nD
+from ..util import img_as_float
 from ._texture import _glcm_loop, _local_binary_pattern
 
-from ..transform import integrate
-
 
 def greycomatrix(image, distances, angles, levels=256, symmetric=False,
                  normed=False):
@@ -294,8 +293,9 @@ def local_binary_pattern(image, P, R, method='default'):
     output = _local_binary_pattern(image, P, R, methods[method.lower()])
     return output
 
-def multiblock_local_binary_pattern(int_image, x, y, width, height):
-    """Multi-block local binary pattern.
+
+def visualize_multiblock_lbp(img, x, y, width, height, lbp_code=0):
+    """Multi-block local binary pattern visualization.
 
     MB-LBP is an extension of LBP that can be computed on many
     scales in a constant time using integral image. It consists of
@@ -304,10 +304,15 @@ def multiblock_local_binary_pattern(int_image, x, y, width, height):
     depending on comparison result, the feature descriptor is
     computed.
 
+    The blocks visualized in the following manner: the center block
+    is left untouched. The blocks that have higher are covered with
+    transparent white rectangles. The blocks that have less intensity
+    are covered with cyan rectangles.
+
     Parameters
     ----------
-    int_image : (N, M) array
-        Integral image.
+    img :
+        Image on which to visualize the pattern.
     x : int
         X-coordinate of top left corner of a rectangle containing feature.
     y : int
@@ -318,11 +323,14 @@ def multiblock_local_binary_pattern(int_image, x, y, width, height):
     height : int
         Height of one of 9 equal rectangles that will be used to compute
         a feature.
+    lbp_code : int
+        The descriptor of feature to visualize. If not provided,
+        the descriptor with 0 value will be used.
 
     Returns
     -------
-    output : int
-        8bit MB-LBP feature descriptor.
+    output :
+        Float image with visualization.
 
     References
     ----------
@@ -332,55 +340,21 @@ def multiblock_local_binary_pattern(int_image, x, y, width, height):
            http://www.cbsr.ia.ac.cn/users/scliao/papers/Zhang-ICB07-MBLBP.pdf
     """
 
-    # Top-left coordinates of central rectangle
-    central_rect_x = x + width
-    central_rect_y = y + height
+    # Default colors for regions.
+    # White is for the blocks that are brighter.
+    # Cyan is for the blocks that has less intensity.
+    color_greater_block = np.asarray([1, 1, 1], dtype='float64')
+    color_less_block = np.asarray([0, 0.69, 0.96], dtype='float64')
 
-    # Sum of intensity values of central rectangle
-    central_rect_val = integrate(int_image,
-                                 central_rect_y,
-                                 central_rect_x,
-                                 central_rect_y + height - 1,
-                                 central_rect_x + width - 1)
+    # Copy array to avoid the changes to the original one
+    output = np.copy(img)
 
-    # Offsets of neighbour rectangles relative to central one.
-    # It has order starting from top left and going clockwise
-    neighbour_rect_offsets = ((-1, -1), (0, -1), (1, -1),
-                              (1, 0), (1, 1), (0, 1),
-                              (-1, 1), (-1, 0))
+    # As the visualization uses RGB color we need 3 bands.
+    if len(img.shape) < 3:
+        output = np.dstack((img,) * 3)
 
-    lbp_code = 0
-
-    for element_num, offset in enumerate(neighbour_rect_offsets):
-
-        offset_x, offset_y = offset
-
-        current_rect_x = central_rect_x + offset_x * width
-        current_rect_y = central_rect_y + offset_y * height
-        current_rect_val = integrate(int_image,
-                                     current_rect_y,
-                                     current_rect_x,
-                                     current_rect_y + height - 1,
-                                     current_rect_x + width - 1)
-
-        has_greater_value = current_rect_val >= central_rect_val
-
-        # If current rectangle's intensity value is bigger
-        # make corresponding bit to 1.
-        lbp_code |= has_greater_value << element_num
-
-        print lbp_code
-
-    return lbp_code
-
-def visualize_multiblock_lbp(img, x, y, width, height, lbp_code=0):
-
-    import matplotlib.patches as patches
-    import matplotlib.pyplot as plt
-
-    plt.imshow(img)
-    img_desc = plt.gca()
-    plt.set_cmap('gray')
+    # Colors are specified in floats
+    output = img_as_float(output)
 
     # Offsets of neighbour rectangles relative to central one.
     # It has order starting from top left and going clockwise
@@ -396,27 +370,19 @@ def visualize_multiblock_lbp(img, x, y, width, height, lbp_code=0):
 
         offset_x, offset_y = offset
 
-        current_rect_x = central_rect_x + offset_x * width
-        current_rect_y = central_rect_y + offset_y * height
+        curr_x = central_rect_x + offset_x * width
+        curr_y = central_rect_y + offset_y * height
 
-        has_greater_value = lbp_code & (1 << element_num)
-
-        # Hatch the rectangles that has less
-        # intensity than the central rectangle.
-        hatch = '\\'
+        has_greater_value = lbp_code & (1 << (7-element_num))
 
+        # Mix-in the visualization colors
         if has_greater_value:
-            hatch = ''
+            output[curr_y:curr_y+height, curr_x:curr_x+width] = \
+                0.5 * output[curr_y:curr_y+height, curr_x:curr_x+width] \
+                + 0.5 * color_greater_block
+        else:
+            output[curr_y:curr_y+height, curr_x:curr_x+width] = \
+                0.5 * output[curr_y:curr_y+height, curr_x:curr_x+width] \
+                + 0.5 * color_less_block
 
-        img_desc.add_patch(
-            patches.Rectangle(
-                (current_rect_x, current_rect_y),
-                width,
-                height,
-                fill=False,
-                hatch=hatch,
-                color='w'
-            )
-        )
-
-    plt.show()
+    return output