diff --git a/.travis.yml b/.travis.yml
index 3a127d3d..69534244 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -19,6 +19,9 @@ install:
     - sudo easy_install$PYSUF pip
     - sudo pip-$PYVER install cython
     - sudo apt-get install libfreeimage3
+    - if [[ $PYVER == '2.7' ]]; then sudo apt-get install $PYTHON-matplotlib; fi
+    - if [[ $PYVER == '3.2' ]]; then sudo pip-$PYVER install git+git://github.com/matplotlib/matplotlib.git@v1.2.x; fi
+    - sudo pip-$PYVER install flake8 --use-mirrors
     - $PYTHON setup.py build
     - sudo $PYTHON setup.py install
 script:
@@ -26,4 +29,8 @@ script:
     - mkdir for_test
     - cd for_test
     - nosetests-$PYVER --exe -v --cover-package=skimage skimage
-
+    # Change back to repository root directory and run all doc examples
+    - cd ..
+    - "echo 'backend : Agg' > matplotlibrc"
+    - for f in doc/examples/*.py; do $PYTHON "$f"; if [ $? -ne 0 ]; then exit $?; fi done
+    - flake8 --exit-zero skimage doc/examples viewer_examples
diff --git a/doc/examples/applications/plot_geometric.py b/doc/examples/applications/plot_geometric.py
index 337ecad7..141ef51f 100644
--- a/doc/examples/applications/plot_geometric.py
+++ b/doc/examples/applications/plot_geometric.py
@@ -7,6 +7,8 @@ In this example, we will see how to use geometric transformations in the context
 of image processing.
 """
 
+from __future__ import print_function
+
 import math
 import numpy as np
 import matplotlib.pyplot as plt
@@ -31,7 +33,7 @@ First we create a transformation using explicit parameters:
 
 tform = tf.SimilarityTransform(scale=1, rotation=math.pi / 2,
                                translation=(0, 1))
-print tform._matrix
+print(tform._matrix)
 
 """
 Alternatively you can define a transformation by the transformation matrix
@@ -49,8 +51,8 @@ systems:
 """
 
 coord = [1, 0]
-print tform2(coord)
-print tform2.inverse(tform(coord))
+print(tform2(coord))
+print(tform2.inverse(tform(coord)))
 
 """
 Image warping
diff --git a/doc/examples/plot_circular_hough_transform.py b/doc/examples/plot_circular_hough_transform.py
index 3344e0d1..fe4bea26 100755
--- a/doc/examples/plot_circular_hough_transform.py
+++ b/doc/examples/plot_circular_hough_transform.py
@@ -36,10 +36,11 @@ from skimage import data, filter, color
 from skimage.transform import hough_circle
 from skimage.feature import peak_local_max
 from skimage.draw import circle_perimeter
+from skimage.util import img_as_ubyte
 
 
 # Load picture and detect edges
-image = data.coins()[0:95, 70:370]
+image = img_as_ubyte(data.coins()[0:95, 70:370])
 edges = filter.canny(image, sigma=3, low_threshold=10, high_threshold=50)
 
 fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(6, 6))
diff --git a/doc/examples/plot_entropy.py b/doc/examples/plot_entropy.py
index ce7bcb2d..db230d8b 100644
--- a/doc/examples/plot_entropy.py
+++ b/doc/examples/plot_entropy.py
@@ -10,10 +10,11 @@ import matplotlib.pyplot as plt
 from skimage import data
 from skimage.filter.rank import entropy
 from skimage.morphology import disk
+from skimage.util import img_as_ubyte
 
 
 # defining a 8- and a 16-bit test images
-a8 = data.camera()
+a8 = img_as_ubyte(data.camera())
 a16 = a8.astype(np.uint16) * 4
 
 ent8 = entropy(a8, disk(5)) # pixel value contain 10x the local entropy
diff --git a/doc/examples/plot_gabor.py b/doc/examples/plot_gabor.py
index c335016d..db1ff7d6 100644
--- a/doc/examples/plot_gabor.py
+++ b/doc/examples/plot_gabor.py
@@ -12,6 +12,8 @@ kernels. The mean and variance of the filtered images are then used as features
 for classification, which is based on the least squared error for simplicity.
 
 """
+from __future__ import print_function
+
 import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
@@ -69,19 +71,19 @@ ref_feats[0, :, :] = compute_feats(brick, kernels)
 ref_feats[1, :, :] = compute_feats(grass, kernels)
 ref_feats[2, :, :] = compute_feats(wall, kernels)
 
-print 'Rotated images matched against references using Gabor filter banks:'
+print('Rotated images matched against references using Gabor filter banks:')
 
-print 'original: brick, rotated: 30deg, match result:',
+print('original: brick, rotated: 30deg, match result: ', end='')
 feats = compute_feats(nd.rotate(brick, angle=190, reshape=False), kernels)
-print image_names[match(feats, ref_feats)]
+print(image_names[match(feats, ref_feats)])
 
-print 'original: brick, rotated: 70deg, match result:',
+print('original: brick, rotated: 70deg, match result: ', end='')
 feats = compute_feats(nd.rotate(brick, angle=70, reshape=False), kernels)
-print image_names[match(feats, ref_feats)]
+print(image_names[match(feats, ref_feats)])
 
-print 'original: grass, rotated: 145deg, match result:',
+print('original: grass, rotated: 145deg, match result: ', end='')
 feats = compute_feats(nd.rotate(grass, angle=145, reshape=False), kernels)
-print image_names[match(feats, ref_feats)]
+print(image_names[match(feats, ref_feats)])
 
 
 def power(image, kernel):
diff --git a/doc/examples/plot_ihc_color_separation.py b/doc/examples/plot_ihc_color_separation.py
index b0b4a16c..f89f28fe 100644
--- a/doc/examples/plot_ihc_color_separation.py
+++ b/doc/examples/plot_ihc_color_separation.py
@@ -22,6 +22,7 @@ import matplotlib.pyplot as plt
 from skimage import data
 from skimage.color import rgb2hed
 
+
 ihc_rgb = data.immunohistochemistry()
 ihc_hed = rgb2hed(ihc_rgb)
 
diff --git a/doc/examples/plot_local_binary_pattern.py b/doc/examples/plot_local_binary_pattern.py
index d32a642f..8c9cabae 100644
--- a/doc/examples/plot_local_binary_pattern.py
+++ b/doc/examples/plot_local_binary_pattern.py
@@ -9,9 +9,12 @@ textures. For simplicity the histogram distributions are then tested against
 each other using the Kullback-Leibler-Divergence.
 
 """
+from __future__ import print_function
+
 import numpy as np
 import matplotlib
 import matplotlib.pyplot as plt
+
 from skimage.transform import rotate
 from skimage.feature import local_binary_pattern
 from skimage import data
@@ -57,13 +60,13 @@ refs = {
 }
 
 # classify rotated textures
-print 'Rotated images matched against references using LBP:'
-print 'original: brick, rotated: 30deg, match result:',
-print match(refs, rotate(brick, angle=30, resize=False))
-print 'original: brick, rotated: 70deg, match result:',
-print match(refs, rotate(brick, angle=70, resize=False))
-print 'original: grass, rotated: 145deg, match result:',
-print match(refs, rotate(grass, angle=145, resize=False))
+print('Rotated images matched against references using LBP:')
+print('original: brick, rotated: 30deg, match result: ', end='')
+print(match(refs, rotate(brick, angle=30, resize=False)))
+print('original: brick, rotated: 70deg, match result: ', end='')
+print(match(refs, rotate(brick, angle=70, resize=False)))
+print('original: grass, rotated: 145deg, match result: ', end='')
+print(match(refs, rotate(grass, angle=145, resize=False)))
 
 # plot histograms of LBP of textures
 fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows=2, ncols=3,
diff --git a/doc/examples/plot_local_equalize.py b/doc/examples/plot_local_equalize.py
index efa04c33..1fd8325f 100644
--- a/doc/examples/plot_local_equalize.py
+++ b/doc/examples/plot_local_equalize.py
@@ -19,15 +19,14 @@ References
 .. [2] http://en.wikipedia.org/wiki/Adaptive_histogram_equalization
 
 """
+import numpy as np
+import matplotlib.pyplot as plt
 
 from skimage import data
 from skimage.util.dtype import dtype_range
+from skimage.util import img_as_ubyte
 from skimage import exposure
 from skimage.morphology import disk
-
-import matplotlib.pyplot as plt
-
-import numpy as np
 from skimage.filter import rank
 
 
@@ -58,7 +57,7 @@ def plot_img_and_hist(img, axes, bins=256):
 
 
 # Load an example image
-img = data.moon()
+img = img_as_ubyte(data.moon())
 
 # Contrast stretching
 p2 = np.percentile(img, 2)
diff --git a/doc/examples/plot_local_otsu.py b/doc/examples/plot_local_otsu.py
index aee0ef6a..3a321780 100644
--- a/doc/examples/plot_local_otsu.py
+++ b/doc/examples/plot_local_otsu.py
@@ -18,12 +18,12 @@ The example compares the local threshold with the global threshold.
 import matplotlib.pyplot as plt
 
 from skimage import data
-from skimage.morphology.selem import disk
-import skimage.filter.rank as rank
-from skimage.filter import threshold_otsu
+from skimage.morphology import disk
+from skimage.filter import threshold_otsu, rank
+from skimage.util import img_as_ubyte
 
 
-p8 = data.page()
+p8 = img_as_ubyte(data.page())
 
 radius = 10
 selem = disk(radius)
diff --git a/doc/examples/plot_marked_watershed.py b/doc/examples/plot_marked_watershed.py
index 31e7a245..d7f2c354 100644
--- a/doc/examples/plot_marked_watershed.py
+++ b/doc/examples/plot_marked_watershed.py
@@ -16,13 +16,14 @@ See Wikipedia_ for more details on the algorithm.
 
 from scipy import ndimage
 import matplotlib.pyplot as plt
+
 from skimage.morphology import watershed, disk
 from skimage import data
-
-# original data
 from skimage.filter import rank
+from skimage.util import img_as_ubyte
 
-image = data.camera()
+
+image = img_as_ubyte(data.camera())
 
 # denoise image
 denoised = rank.median(image, disk(2))
diff --git a/doc/examples/plot_polygon.py b/doc/examples/plot_polygon.py
index 078bb6cf..821c4558 100644
--- a/doc/examples/plot_polygon.py
+++ b/doc/examples/plot_polygon.py
@@ -7,8 +7,11 @@ This example shows how to approximate (Douglas-Peucker algorithm) and subdivide
 (B-Splines) polygonal chains.
 
 """
+from __future__ import print_function
+
 import numpy as np
 import matplotlib.pyplot as plt
+
 from skimage.draw import ellipse
 from skimage.measure import find_contours, approximate_polygon, \
     subdivide_polygon
@@ -45,7 +48,7 @@ for _ in range(5):
 # approximate subdivided polygon with Douglas-Peucker algorithm
 appr_hand = approximate_polygon(new_hand, tolerance=0.02)
 
-print "Number of coordinates:", len(hand), len(new_hand), len(appr_hand)
+print("Number of coordinates:", len(hand), len(new_hand), len(appr_hand))
 
 fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(9, 4))
 
@@ -70,7 +73,7 @@ for contour in find_contours(img, 0):
     ax2.plot(coords[:, 1], coords[:, 0], '-r', linewidth=2)
     coords2 = approximate_polygon(contour, tolerance=39.5)
     ax2.plot(coords2[:, 1], coords2[:, 0], '-g', linewidth=2)
-    print "Number of coordinates:", len(contour), len(coords), len(coords2)
+    print("Number of coordinates:", len(contour), len(coords), len(coords2))
 
 ax2.axis((0, 800, 0, 800))
 
diff --git a/doc/examples/plot_segmentations.py b/doc/examples/plot_segmentations.py
index 418a0903..55bb6082 100644
--- a/doc/examples/plot_segmentations.py
+++ b/doc/examples/plot_segmentations.py
@@ -58,6 +58,7 @@ of Quickshift, while ``n_segments`` chooses the number of centers for kmeans.
     Pascal Fua, and Sabine Suesstrunk, SLIC Superpixels Compared to
     State-of-the-art Superpixel Methods, TPAMI, May 2012.
 """
+from __future__ import print_function
 
 import matplotlib.pyplot as plt
 import numpy as np
diff --git a/skimage/color/tests/test_colorconv.py b/skimage/color/tests/test_colorconv.py
index 3365bb1c..7f1d7c73 100644
--- a/skimage/color/tests/test_colorconv.py
+++ b/skimage/color/tests/test_colorconv.py
@@ -126,7 +126,7 @@ class TestColorconv(TestCase):
 
     # RGB<->HED roundtrip with ubyte image
     def test_hed_rgb_roundtrip(self):
-        img_rgb = self.img_rgb
+        img_rgb = img_as_ubyte(self.img_rgb)
         assert_equal(img_as_ubyte(hed2rgb(rgb2hed(img_rgb))), img_rgb)
 
     # RGB<->HED roundtrip with float image
diff --git a/skimage/data/__init__.py b/skimage/data/__init__.py
index 07d2688d..864e3638 100644
--- a/skimage/data/__init__.py
+++ b/skimage/data/__init__.py
@@ -169,7 +169,7 @@ def immunohistochemistry():
     No known copyright restrictions.
 
     """
-    return load("ihc.jpg")
+    return load("ihc.png")
 
 
 def chelsea():
@@ -183,4 +183,4 @@ def chelsea():
     No copyright restrictions.  CC0 by the photographer (Stefan van der Walt).
 
     """
-    return load("chelsea.jpg")
+    return load("chelsea.png")
diff --git a/skimage/data/chelsea.jpg b/skimage/data/chelsea.jpg
deleted file mode 100644
index 7692c1f5..00000000
Binary files a/skimage/data/chelsea.jpg and /dev/null differ
diff --git a/skimage/data/chelsea.png b/skimage/data/chelsea.png
new file mode 100644
index 00000000..fae212d2
Binary files /dev/null and b/skimage/data/chelsea.png differ
diff --git a/skimage/data/ihc.jpg b/skimage/data/ihc.jpg
deleted file mode 100644
index e64d2a0d..00000000
Binary files a/skimage/data/ihc.jpg and /dev/null differ
diff --git a/skimage/data/ihc.png b/skimage/data/ihc.png
new file mode 100644
index 00000000..59df35fd
Binary files /dev/null and b/skimage/data/ihc.png differ