diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index bc58df84..c94755d8 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -215,3 +215,6 @@ - Jim Fienup, Alexander Iacchetta In-depth review of sub-pixel shift registration + +- Damian Eads + Structuring elements in morphology module. diff --git a/RELEASE.txt b/RELEASE.txt index 915134ad..b5b34cb9 100644 --- a/RELEASE.txt +++ b/RELEASE.txt @@ -22,9 +22,6 @@ How to make a new release of ``skimage`` - Edit ``doc/source/_static/docversions.js`` and commit - Build a clean version of the docs. Run ``python setup.py install`` in the root dir, then ``rm -rf build; make html`` in the docs. - - Run ``make html`` again to copy the newly generated ``random.js`` into - place. Double check ``random.js``, otherwise the skimage.org front - page gets broken! - Build using ``make gh-pages``. - Update the symlink to ``stable``. - Push upstream: ``git push origin gh-pages`` in ``doc/gh-pages``. diff --git a/TODO.txt b/TODO.txt index 2b890b89..02bd9ead 100644 --- a/TODO.txt +++ b/TODO.txt @@ -8,6 +8,10 @@ Version 0.14 * Remove deprecated ``skimage.restoration.nl_means_denoising``. * Remove deprecated ``skimage.filters.gaussian_filter``. * Remove deprecated ``skimage.filters.gabor_filter``. +* Remove deprecated ``skimage.measure.LineModel`` and + add an alias LineModel = LineModelND. While the deprecated LineModel has for + parameters `(dist, theta)`, LineModelND has the more general parameters + `(origin, direction)`. Version 0.13 diff --git a/doc/Makefile b/doc/Makefile index 589e2937..63523365 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -13,7 +13,7 @@ PAPEROPT_a4 = -D latex_paper_size=a4 PAPEROPT_letter = -D latex_paper_size=letter ALLSPHINXOPTS = -d $(SPHINXCACHE) $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source DEST = build -.PHONY: all help clean html dirhtml pickle json htmlhelp qthelp latex changes linkcheck doctest gitwash gh-pages release_notes random_gallery +.PHONY: all help clean html dirhtml pickle json htmlhelp qthelp latex changes linkcheck doctest gitwash gh-pages release_notes all: html @@ -41,18 +41,16 @@ api: $(PYTHON) tools/build_modref_templates.py @echo "Build API docs...done." -random_gallery: - @cd source && $(PYTHON) random_gallery.py - release_notes: @echo "Copying release notes" @tail -n +4 `ls release/*.txt | sort -k 2 -t . -n | tail -n 1` > release/_release_notes_for_docs.txt -html: api release_notes random_gallery +html: api release_notes $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(DEST)/html cp -r source/plots $(DEST)/html @echo @echo "Build finished. The HTML pages are in build/html." + $(PYTHON) source/random_gallery.py dirhtml: $(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(DEST)/dirhtml diff --git a/doc/README.md b/doc/README.md index 372d84e5..d9fec89a 100644 --- a/doc/README.md +++ b/doc/README.md @@ -25,9 +25,10 @@ sudo tlmgr install ucs dvipng ## Fixing Warnings ## - "citation not found: R###" - $ cd doc/build; grep -rin R### . There is probably an underscore after a reference - in the first line of a docstring (e.g. [1]_) + in the first line of a docstring (e.g. [1]_). + Use this method to find the source file: + $ cd doc/build; grep -rin R#### - "Duplicate citation R###, other instance in..."" There is probably a [2] without a [1] in one of diff --git a/doc/examples/applications/plot_coins_segmentation.py b/doc/examples/applications/plot_coins_segmentation.py index fe165760..404c1b98 100644 --- a/doc/examples/applications/plot_coins_segmentation.py +++ b/doc/examples/applications/plot_coins_segmentation.py @@ -35,13 +35,15 @@ background with the coins: """ -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3)) +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3), sharex=True, sharey=True) ax1.imshow(coins > 100, cmap=plt.cm.gray, interpolation='nearest') ax1.set_title('coins > 100') ax1.axis('off') +ax1.set_adjustable('box-forced') ax2.imshow(coins > 150, cmap=plt.cm.gray, interpolation='nearest') ax2.set_title('coins > 150') ax2.axis('off') +ax2.set_adjustable('box-forced') margins = dict(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1) fig.subplots_adjust(**margins) @@ -162,12 +164,14 @@ segmentation = ndi.binary_fill_holes(segmentation - 1) labeled_coins, _ = ndi.label(segmentation) image_label_overlay = label2rgb(labeled_coins, image=coins) -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3)) +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3), sharex=True, sharey=True) ax1.imshow(coins, cmap=plt.cm.gray, interpolation='nearest') ax1.contour(segmentation, [0.5], linewidths=1.2, colors='y') ax1.axis('off') +ax1.set_adjustable('box-forced') ax2.imshow(image_label_overlay, interpolation='nearest') ax2.axis('off') +ax2.set_adjustable('box-forced') fig.subplots_adjust(**margins) diff --git a/doc/examples/applications/plot_morphology.py b/doc/examples/applications/plot_morphology.py index 68ecc6d6..33246793 100644 --- a/doc/examples/applications/plot_morphology.py +++ b/doc/examples/applications/plot_morphology.py @@ -42,13 +42,15 @@ Let's also define a convenience function for plotting comparisons: def plot_comparison(original, filtered, filter_name): - fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4)) + fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4), sharex=True, sharey=True) ax1.imshow(original, cmap=plt.cm.gray) ax1.set_title('original') ax1.axis('off') + ax1.set_adjustable('box-forced') ax2.imshow(filtered, cmap=plt.cm.gray) ax2.set_title(filter_name) ax2.axis('off') + ax2.set_adjustable('box-forced') """ Erosion diff --git a/doc/examples/applications/plot_rank_filters.py b/doc/examples/applications/plot_rank_filters.py index e2dbd71b..b1e0cc79 100644 --- a/doc/examples/applications/plot_rank_filters.py +++ b/doc/examples/applications/plot_rank_filters.py @@ -70,24 +70,31 @@ noisy_image = img_as_ubyte(data.camera()) noisy_image[noise > 0.99] = 255 noisy_image[noise < 0.01] = 0 -fig, ax = plt.subplots(2, 2, figsize=(10, 7)) +fig, ax = plt.subplots(2, 2, figsize=(10, 7), sharex=True, sharey=True) ax1, ax2, ax3, ax4 = ax.ravel() ax1.imshow(noisy_image, vmin=0, vmax=255, cmap=plt.cm.gray) ax1.set_title('Noisy image') ax1.axis('off') +ax1.set_adjustable('box-forced') ax2.imshow(median(noisy_image, disk(1)), vmin=0, vmax=255, cmap=plt.cm.gray) ax2.set_title('Median $r=1$') ax2.axis('off') +ax2.set_adjustable('box-forced') + ax3.imshow(median(noisy_image, disk(5)), vmin=0, vmax=255, cmap=plt.cm.gray) ax3.set_title('Median $r=5$') ax3.axis('off') +ax3.set_adjustable('box-forced') + ax4.imshow(median(noisy_image, disk(20)), vmin=0, vmax=255, cmap=plt.cm.gray) ax4.set_title('Median $r=20$') ax4.axis('off') +ax4.set_adjustable('box-forced') + """ @@ -109,17 +116,19 @@ image. from skimage.filters.rank import mean -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 7]) +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 7], sharex=True, sharey=True) loc_mean = mean(noisy_image, disk(10)) ax1.imshow(noisy_image, vmin=0, vmax=255, cmap=plt.cm.gray) ax1.set_title('Original') ax1.axis('off') +ax1.set_adjustable('box-forced') ax2.imshow(loc_mean, vmin=0, vmax=255, cmap=plt.cm.gray) ax2.set_title('Local mean $r=10$') ax2.axis('off') +ax2.set_adjustable('box-forced') """ @@ -143,22 +152,26 @@ noisy_image = img_as_ubyte(data.camera()) bilat = mean_bilateral(noisy_image.astype(np.uint16), disk(20), s0=10, s1=10) -fig, ax = plt.subplots(2, 2, figsize=(10, 7)) +fig, ax = plt.subplots(2, 2, figsize=(10, 7), sharex='row', sharey='row') ax1, ax2, ax3, ax4 = ax.ravel() ax1.imshow(noisy_image, cmap=plt.cm.gray) ax1.set_title('Original') ax1.axis('off') +ax1.set_adjustable('box-forced') ax2.imshow(bilat, cmap=plt.cm.gray) ax2.set_title('Bilateral mean') ax2.axis('off') +ax2.set_adjustable('box-forced') ax3.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray) ax3.axis('off') +ax3.set_adjustable('box-forced') ax4.imshow(bilat[200:350, 350:450], cmap=plt.cm.gray) ax4.axis('off') +ax4.set_adjustable('box-forced') """ @@ -236,15 +249,17 @@ noisy_image = img_as_ubyte(data.camera()) auto = autolevel(noisy_image.astype(np.uint16), disk(20)) -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 7]) +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 7], sharex=True, sharey=True) ax1.imshow(noisy_image, cmap=plt.cm.gray) ax1.set_title('Original') ax1.axis('off') +ax1.set_adjustable('box-forced') ax2.imshow(auto, cmap=plt.cm.gray) ax2.set_title('Local autolevel') ax2.axis('off') +ax2.set_adjustable('box-forced') """ @@ -271,22 +286,29 @@ loc_perc_autolevel1 = autolevel_percentile(image, selem=selem, p0=.01, p1=.99) loc_perc_autolevel2 = autolevel_percentile(image, selem=selem, p0=.05, p1=.95) loc_perc_autolevel3 = autolevel_percentile(image, selem=selem, p0=.1, p1=.9) -fig, axes = plt.subplots(nrows=3, figsize=(7, 8)) +fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(7, 8), sharex=True, sharey=True) ax0, ax1, ax2 = axes plt.gray() -ax0.imshow(np.hstack((image, loc_autolevel)), cmap=plt.cm.gray) -ax0.set_title('Original / auto-level') +title_list = ['Original', + 'auto_level', + 'auto-level 0%', + 'auto-level 1%', + 'auto-level 5%', + 'auto-level 10%'] +image_list = [image, + loc_autolevel, + loc_perc_autolevel0, + loc_perc_autolevel1, + loc_perc_autolevel2, + loc_perc_autolevel3] +axes_list = axes.ravel().tolist() -ax1.imshow( - np.hstack((loc_perc_autolevel0, loc_perc_autolevel1)), vmin=0, vmax=255) -ax1.set_title('Percentile auto-level 0%,1%') -ax2.imshow( - np.hstack((loc_perc_autolevel2, loc_perc_autolevel3)), vmin=0, vmax=255) -ax2.set_title('Percentile auto-level 5% and 10%') - -for ax in axes: - ax.axis('off') +for i in range(0,len(image_list)): + axes_list[i].imshow(image_list[i], cmap=plt.cm.gray, vmin=0, vmax=255) + axes_list[i].set_title(title_list[i]) + axes_list[i].axis('off') + axes_list[i].set_adjustable('box-forced') """ @@ -304,22 +326,26 @@ noisy_image = img_as_ubyte(data.camera()) enh = enhance_contrast(noisy_image, disk(5)) -fig, ax = plt.subplots(2, 2, figsize=[10, 7]) +fig, ax = plt.subplots(2, 2, figsize=[10, 7], sharex='row', sharey='row') ax1, ax2, ax3, ax4 = ax.ravel() ax1.imshow(noisy_image, cmap=plt.cm.gray) ax1.set_title('Original') ax1.axis('off') +ax1.set_adjustable('box-forced') ax2.imshow(enh, cmap=plt.cm.gray) ax2.set_title('Local morphological contrast enhancement') ax2.axis('off') +ax2.set_adjustable('box-forced') ax3.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray) ax3.axis('off') +ax3.set_adjustable('box-forced') ax4.imshow(enh[200:350, 350:450], cmap=plt.cm.gray) ax4.axis('off') +ax4.set_adjustable('box-forced') """ @@ -336,22 +362,22 @@ noisy_image = img_as_ubyte(data.camera()) penh = enhance_contrast_percentile(noisy_image, disk(5), p0=.1, p1=.9) -fig, ax = plt.subplots(2, 2, figsize=[10, 7]) +fig, ax = plt.subplots(2, 2, figsize=[10, 7], sharex='row', sharey='row') ax1, ax2, ax3, ax4 = ax.ravel() ax1.imshow(noisy_image, cmap=plt.cm.gray) ax1.set_title('Original') -ax1.axis('off') ax2.imshow(penh, cmap=plt.cm.gray) ax2.set_title('Local percentile morphological\n contrast enhancement') -ax2.axis('off') ax3.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray) -ax3.axis('off') ax4.imshow(penh[200:350, 350:450], cmap=plt.cm.gray) -ax4.axis('off') + +for ax in ax.ravel(): + ax.axis('off') + ax.set_adjustable('box-forced') """ @@ -393,24 +419,24 @@ loc_otsu = p8 >= t_loc_otsu t_glob_otsu = threshold_otsu(p8) glob_otsu = p8 >= t_glob_otsu -fig, ax = plt.subplots(2, 2) +fig, ax = plt.subplots(2, 2, sharex=True, sharey=True) ax1, ax2, ax3, ax4 = ax.ravel() fig.colorbar(ax1.imshow(p8, cmap=plt.cm.gray), ax=ax1) ax1.set_title('Original') -ax1.axis('off') fig.colorbar(ax2.imshow(t_loc_otsu, cmap=plt.cm.gray), ax=ax2) ax2.set_title('Local Otsu ($r=%d$)' % radius) -ax2.axis('off') ax3.imshow(p8 >= t_loc_otsu, cmap=plt.cm.gray) ax3.set_title('Original >= local Otsu' % t_glob_otsu) -ax3.axis('off') ax4.imshow(glob_otsu, cmap=plt.cm.gray) ax4.set_title('Global Otsu ($t=%d$)' % t_glob_otsu) -ax4.axis('off') + +for ax in ax.ravel(): + ax.axis('off') + ax.set_adjustable('box-forced') """ @@ -429,15 +455,17 @@ m = (np.tile(x, (n, 1)) * np.linspace(0.1, 1, n) * 128 + 128).astype(np.uint8) radius = 10 t = rank.otsu(m, disk(radius)) -fig, (ax1, ax2) = plt.subplots(1, 2) +fig, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True) ax1.imshow(m) ax1.set_title('Original') ax1.axis('off') +ax1.set_adjustable('box-forced') ax2.imshow(m >= t, interpolation='nearest') ax2.set_title('Local Otsu ($r=%d$)' % radius) ax2.axis('off') +ax2.set_adjustable('box-forced') """ @@ -468,25 +496,24 @@ opening = minimum(maximum(noisy_image, disk(5)), disk(5)) grad = gradient(noisy_image, disk(5)) # display results -fig, ax = plt.subplots(2, 2, figsize=[10, 7]) +fig, ax = plt.subplots(2, 2, figsize=[10, 7], sharex=True, sharey=True) ax1, ax2, ax3, ax4 = ax.ravel() ax1.imshow(noisy_image, cmap=plt.cm.gray) ax1.set_title('Original') -ax1.axis('off') ax2.imshow(closing, cmap=plt.cm.gray) ax2.set_title('Gray-level closing') -ax2.axis('off') ax3.imshow(opening, cmap=plt.cm.gray) ax3.set_title('Gray-level opening') -ax3.axis('off') ax4.imshow(grad, cmap=plt.cm.gray) ax4.set_title('Morphological gradient') -ax4.axis('off') +for ax in ax.ravel(): + ax.axis('off') + ax.set_adjustable('box-forced') """ .. image:: PLOT2RST.current_figure @@ -518,15 +545,17 @@ import matplotlib.pyplot as plt image = data.camera() -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4)) +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), sharex=True, sharey=True) fig.colorbar(ax1.imshow(image, cmap=plt.cm.gray), ax=ax1) ax1.set_title('Image') ax1.axis('off') +ax1.set_adjustable('box-forced') fig.colorbar(ax2.imshow(entropy(image, disk(5)), cmap=plt.cm.jet), ax=ax2) ax2.set_title('Entropy') ax2.axis('off') +ax2.set_adjustable('box-forced') """ @@ -680,10 +709,15 @@ Comparison of outcome of the three methods: """ -fig, ax = plt.subplots() -ax.imshow(np.hstack((rc, rndi))) -ax.set_title('filters.rank.median vs. scipy.ndimage.percentile') -ax.axis('off') +fig, (ax0, ax1) = plt.subplots(ncols=2, sharex=True, sharey=True) +ax0.set_title('filters.rank.median') +ax0.imshow(rc) +ax0.axis('off') +ax0.set_adjustable('box-forced') +ax1.set_title('scipy.ndimage.percentile') +ax1.imshow(rndi) +ax1.axis('off') +ax1.set_adjustable('box-forced') """ .. image:: PLOT2RST.current_figure diff --git a/doc/examples/plot_adapt_rgb.py b/doc/examples/plot_adapt_rgb.py index d267cc4a..24d8f078 100644 --- a/doc/examples/plot_adapt_rgb.py +++ b/doc/examples/plot_adapt_rgb.py @@ -48,8 +48,8 @@ import matplotlib.pyplot as plt image = data.astronaut() fig = plt.figure(figsize=(14, 7)) -ax_each = fig.add_subplot(121) -ax_hsv = fig.add_subplot(122) +ax_each = fig.add_subplot(121, adjustable='box-forced') +ax_hsv = fig.add_subplot(122, sharex=ax_each, sharey=ax_each, adjustable='box-forced') # We use 1 - sobel_each(image) # but this will not work if image is not normalized @@ -107,7 +107,7 @@ def sobel_gray(image): return filters.sobel(image) fig = plt.figure(figsize=(7, 7)) -ax = fig.add_subplot(111) +ax = fig.add_subplot(111, sharex=ax_each, sharey=ax_each, adjustable='box-forced') # We use 1 - sobel_gray(image) # but this will not work if image is not normalized diff --git a/doc/examples/plot_blob.py b/doc/examples/plot_blob.py index 33cc44c2..bb5dd089 100644 --- a/doc/examples/plot_blob.py +++ b/doc/examples/plot_blob.py @@ -61,8 +61,12 @@ titles = ['Laplacian of Gaussian', 'Difference of Gaussian', 'Determinant of Hessian'] sequence = zip(blobs_list, colors, titles) + +fig,axes = plt.subplots(1, 3, sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +axes = axes.ravel() for blobs, color, title in sequence: - fig, ax = plt.subplots(1, 1) + ax = axes[0] + axes = axes[1:] ax.set_title(title) ax.imshow(image, interpolation='nearest') for blob in blobs: diff --git a/doc/examples/plot_canny.py b/doc/examples/plot_canny.py index 7d770787..06d2a05d 100644 --- a/doc/examples/plot_canny.py +++ b/doc/examples/plot_canny.py @@ -35,7 +35,7 @@ edges1 = feature.canny(im) edges2 = feature.canny(im, sigma=3) # display results -fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(8, 3)) +fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(8, 3), sharex=True, sharey=True) ax1.imshow(im, cmap=plt.cm.jet) ax1.axis('off') diff --git a/doc/examples/plot_circular_elliptical_hough_transform.py b/doc/examples/plot_circular_elliptical_hough_transform.py index 684ad30b..5d5a8b2e 100755 --- a/doc/examples/plot_circular_elliptical_hough_transform.py +++ b/doc/examples/plot_circular_elliptical_hough_transform.py @@ -138,7 +138,7 @@ image_rgb[cy, cx] = (0, 0, 255) edges = color.gray2rgb(edges) edges[cy, cx] = (250, 0, 0) -fig2, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(8, 4)) +fig2, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(8, 4), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax1.set_title('Original picture') ax1.imshow(image_rgb) diff --git a/doc/examples/plot_denoise.py b/doc/examples/plot_denoise.py index 1313bc04..f39591d2 100644 --- a/doc/examples/plot_denoise.py +++ b/doc/examples/plot_denoise.py @@ -38,7 +38,7 @@ astro = astro[220:300, 220:320] noisy = astro + 0.6 * astro.std() * np.random.random(astro.shape) noisy = np.clip(noisy, 0, 1) -fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5)) +fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) plt.gray() diff --git a/doc/examples/plot_edge_filter.py b/doc/examples/plot_edge_filter.py index 96b654a8..7ad1a6eb 100644 --- a/doc/examples/plot_edge_filter.py +++ b/doc/examples/plot_edge_filter.py @@ -19,7 +19,7 @@ image = camera() edge_roberts = roberts(image) edge_sobel = sobel(image) -fig, (ax0, ax1) = plt.subplots(ncols=2) +fig, (ax0, ax1) = plt.subplots(ncols=2, sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax0.imshow(edge_roberts, cmap=plt.cm.gray) ax0.set_title('Roberts Edge Detection') @@ -66,7 +66,7 @@ diff_scharr_prewitt = edge_scharr - edge_prewitt diff_scharr_sobel = edge_scharr - edge_sobel max_diff = np.max(np.maximum(diff_scharr_prewitt, diff_scharr_sobel)) -fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(nrows=2, ncols=2) +fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(nrows=2, ncols=2, sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax0.imshow(img, cmap=plt.cm.gray) ax0.set_title('Original image') diff --git a/doc/examples/plot_entropy.py b/doc/examples/plot_entropy.py index f33f26ff..9bef036e 100644 --- a/doc/examples/plot_entropy.py +++ b/doc/examples/plot_entropy.py @@ -3,30 +3,71 @@ Entropy ======= -Image entropy is a quantity which is used to describe the amount of information -coded in an image. +In information theory, information entropy is the log-base-2 of the number of +possible outcomes for a message. + +For an image, local entropy is related to the complexity contained in a given +neighborhood, typically defined by a structuring element. The entropy filter can +detect subtle variations in the local gray level distribution. + +In the first example, the image is composed of two surfaces with two slightly +different distributions. The image has a uniform random distribution in the +range [-14, +14] in the middle of the image and a uniform random distribution in +the range [-15, 15] at the image borders, both centered at a gray value of 128. +To detect the central square, we compute the local entropy measure using a +circular structuring element of a radius big enough to capture the local gray +level distribution. The second example shows how to detect texture in the camera +image using a smaller structuring element. """ import matplotlib.pyplot as plt +import numpy as np from skimage import data +from skimage.util import img_as_ubyte from skimage.filters.rank import entropy from skimage.morphology import disk -from skimage.util import img_as_ubyte +# First example: object detection. + +noise_mask = 28 * np.ones((128, 128), dtype=np.uint8) +noise_mask[32:-32, 32:-32] = 30 + +noise = (noise_mask * np.random.random(noise_mask.shape) - 0.5 * + noise_mask).astype(np.uint8) +img = noise + 128 + +entr_img = entropy(img, disk(10)) + +fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(8, 3)) + +ax0.imshow(noise_mask, cmap=plt.cm.gray) +ax0.set_xlabel("Noise mask") +ax1.imshow(img, cmap=plt.cm.gray) +ax1.set_xlabel("Noisy image") +ax2.imshow(entr_img) +ax2.set_xlabel("Local entropy") + +fig.tight_layout() + +# Second example: texture detection. image = img_as_ubyte(data.camera()) -fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 4)) +fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 4), sharex=True, + sharey=True, + subplot_kw={"adjustable": "box-forced"}) img0 = ax0.imshow(image, cmap=plt.cm.gray) -ax0.set_title('Image') -ax0.axis('off') +ax0.set_title("Image") +ax0.axis("off") fig.colorbar(img0, ax=ax0) img1 = ax1.imshow(entropy(image, disk(5)), cmap=plt.cm.jet) -ax1.set_title('Entropy') -ax1.axis('off') +ax1.set_title("Entropy") +ax1.axis("off") fig.colorbar(img1, ax=ax1) +fig.tight_layout() + plt.show() diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py index 2d04e22a..2289e270 100644 --- a/doc/examples/plot_equalize.py +++ b/doc/examples/plot_equalize.py @@ -40,6 +40,7 @@ def plot_img_and_hist(img, axes, bins=256): # Display image ax_img.imshow(img, cmap=plt.cm.gray) ax_img.set_axis_off() + ax_img.set_adjustable('box-forced') # Display histogram ax_hist.hist(img.ravel(), bins=bins, histtype='step', color='black') @@ -70,7 +71,13 @@ img_eq = exposure.equalize_hist(img) img_adapteq = exposure.equalize_adapthist(img, clip_limit=0.03) # Display results -fig, axes = plt.subplots(nrows=2, ncols=4, figsize=(8, 5)) +fig = plt.figure(figsize=(8, 5)) +axes = np.zeros((2,4), dtype=np.object) +axes[0,0] = fig.add_subplot(2, 4, 1) +for i in range(1,4): + axes[0,i] = fig.add_subplot(2, 4, 1+i, sharex=axes[0,0], sharey=axes[0,0]) +for i in range(0,4): + axes[1,i] = fig.add_subplot(2, 4, 5+i) ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0]) ax_img.set_title('Low contrast image') diff --git a/doc/examples/plot_hog.py b/doc/examples/plot_hog.py index 15b279bd..ab71bb36 100644 --- a/doc/examples/plot_hog.py +++ b/doc/examples/plot_hog.py @@ -90,11 +90,12 @@ image = color.rgb2gray(data.astronaut()) fd, hog_image = hog(image, orientations=8, pixels_per_cell=(16, 16), cells_per_block=(1, 1), visualise=True) -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), sharex=True, sharey=True) ax1.axis('off') ax1.imshow(image, cmap=plt.cm.gray) ax1.set_title('Input image') +ax1.set_adjustable('box-forced') # Rescale histogram for better display hog_image_rescaled = exposure.rescale_intensity(hog_image, in_range=(0, 0.02)) @@ -102,4 +103,5 @@ hog_image_rescaled = exposure.rescale_intensity(hog_image, in_range=(0, 0.02)) ax2.axis('off') ax2.imshow(hog_image_rescaled, cmap=plt.cm.gray) ax2.set_title('Histogram of Oriented Gradients') +ax1.set_adjustable('box-forced') plt.show() diff --git a/doc/examples/plot_holes_and_peaks.py b/doc/examples/plot_holes_and_peaks.py index 90af1cbe..23c119e8 100644 --- a/doc/examples/plot_holes_and_peaks.py +++ b/doc/examples/plot_holes_and_peaks.py @@ -21,16 +21,13 @@ image = data.moon() # Rescale image intensity so that we can see dim features. image = rescale_intensity(image, in_range=(50, 200)) - -# convenience function for plotting images -def imshow(image, title, **kwargs): - fig, ax = plt.subplots(figsize=(5, 4)) - ax.imshow(image, **kwargs) - ax.axis('off') - ax.set_title(title) +fig,ax = plt.subplots(2, 2, figsize=(5, 4), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +ax = ax.ravel() -imshow(image, 'Original image') +ax[0].imshow(image) +ax[0].set_title('Original image') +ax[0].axis('off') """ .. image:: PLOT2RST.current_figure @@ -52,8 +49,9 @@ mask = image filled = reconstruction(seed, mask, method='erosion') -imshow(filled, 'after filling holes', vmin=image.min(), vmax=image.max()) - +ax[1].imshow(filled) +ax[1].set_title('after filling holes') +ax[1].axis('off') """ .. image:: PLOT2RST.current_figure @@ -63,8 +61,9 @@ isolate the dark regions by subtracting the reconstructed image from the original image. """ -imshow(image - filled, 'holes') -# plt.title('holes') +ax[2].imshow(image-filled) +ax[2].set_title('holes') +ax[2].axis('off') """ .. image:: PLOT2RST.current_figure @@ -79,7 +78,10 @@ intensity instead of the maximum. The remainder of the process is the same. seed = np.copy(image) seed[1:-1, 1:-1] = image.min() rec = reconstruction(seed, mask, method='dilation') -imshow(image - rec, 'peaks') + +ax[3].imshow(image-rec) +ax[3].set_title('peaks') +ax[3].axis('off') plt.show() """ diff --git a/doc/examples/plot_ihc_color_separation.py b/doc/examples/plot_ihc_color_separation.py index 3e297386..6191a800 100644 --- a/doc/examples/plot_ihc_color_separation.py +++ b/doc/examples/plot_ihc_color_separation.py @@ -26,7 +26,7 @@ from skimage.color import rgb2hed ihc_rgb = data.immunohistochemistry() ihc_hed = rgb2hed(ihc_rgb) -fig, axes = plt.subplots(2, 2, figsize=(7, 6)) +fig, axes = plt.subplots(2, 2, figsize=(7, 6), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax0, ax1, ax2, ax3 = axes.ravel() ax0.imshow(ihc_rgb) @@ -61,7 +61,9 @@ h = rescale_intensity(ihc_hed[:, :, 0], out_range=(0, 1)) d = rescale_intensity(ihc_hed[:, :, 2], out_range=(0, 1)) zdh = np.dstack((np.zeros_like(h), d, h)) -fig, ax = plt.subplots() +#fig, ax = plt.subplots() +fig = plt.figure() +ax = plt.subplot(1, 1, 1, sharex=ax0, sharey=ax0, adjustable='box-forced') ax.imshow(zdh) ax.set_title("Stain separated image (rescaled)") ax.axis('off') diff --git a/doc/examples/plot_join_segmentations.py b/doc/examples/plot_join_segmentations.py index 5fded86e..0ad7e57c 100644 --- a/doc/examples/plot_join_segmentations.py +++ b/doc/examples/plot_join_segmentations.py @@ -40,7 +40,7 @@ seg2 = slic(coins, n_segments=117, max_iter=160, sigma=1, compactness=0.75, segj = join_segmentations(seg1, seg2) # show the segmentations -fig, axes = plt.subplots(ncols=4, figsize=(9, 2.5)) +fig, axes = plt.subplots(ncols=4, figsize=(9, 2.5), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) axes[0].imshow(coins, cmap=plt.cm.gray, interpolation='nearest') axes[0].set_title('Image') diff --git a/doc/examples/plot_line_hough_transform.py b/doc/examples/plot_line_hough_transform.py index 7ee3f1f7..15464768 100644 --- a/doc/examples/plot_line_hough_transform.py +++ b/doc/examples/plot_line_hough_transform.py @@ -77,30 +77,30 @@ image[idx, idx] = 255 h, theta, d = hough_line(image) -fig, ax = plt.subplots(1, 3, figsize=(8, 4)) +fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8,4)) -ax[0].imshow(image, cmap=plt.cm.gray) -ax[0].set_title('Input image') -ax[0].axis('image') +ax1.imshow(image, cmap=plt.cm.gray) +ax1.set_title('Input image') +ax1.set_axis_off() -ax[1].imshow(np.log(1 + h), +ax2.imshow(np.log(1 + h), extent=[np.rad2deg(theta[-1]), np.rad2deg(theta[0]), d[-1], d[0]], cmap=plt.cm.gray, aspect=1/1.5) -ax[1].set_title('Hough transform') -ax[1].set_xlabel('Angles (degrees)') -ax[1].set_ylabel('Distance (pixels)') -ax[1].axis('image') +ax2.set_title('Hough transform') +ax2.set_xlabel('Angles (degrees)') +ax2.set_ylabel('Distance (pixels)') +ax2.axis('image') -ax[2].imshow(image, cmap=plt.cm.gray) +ax3.imshow(image, cmap=plt.cm.gray) rows, cols = image.shape for _, angle, dist in zip(*hough_line_peaks(h, theta, d)): y0 = (dist - 0 * np.cos(angle)) / np.sin(angle) y1 = (dist - cols * np.cos(angle)) / np.sin(angle) - ax[2].plot((0, cols), (y0, y1), '-r') -ax[2].axis((0, cols, rows, 0)) -ax[2].set_title('Detected lines') -ax[2].axis('image') + ax3.plot((0, cols), (y0, y1), '-r') +ax3.axis((0, cols, rows, 0)) +ax3.set_title('Detected lines') +ax3.set_axis_off() # Line finding, using the Probabilistic Hough Transform @@ -109,22 +109,25 @@ edges = canny(image, 2, 1, 25) lines = probabilistic_hough_line(edges, threshold=10, line_length=5, line_gap=3) -fig2, ax = plt.subplots(1, 3, figsize=(8, 3)) +fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8,4), sharex=True, sharey=True) -ax[0].imshow(image, cmap=plt.cm.gray) -ax[0].set_title('Input image') -ax[0].axis('image') +ax1.imshow(image, cmap=plt.cm.gray) +ax1.set_title('Input image') +ax1.set_axis_off() +ax1.set_adjustable('box-forced') -ax[1].imshow(edges, cmap=plt.cm.gray) -ax[1].set_title('Canny edges') -ax[1].axis('image') +ax2.imshow(edges, cmap=plt.cm.gray) +ax2.set_title('Canny edges') +ax2.set_axis_off() +ax2.set_adjustable('box-forced') -ax[2].imshow(edges * 0) +ax3.imshow(edges * 0) for line in lines: p0, p1 = line - ax[2].plot((p0[0], p1[0]), (p0[1], p1[1])) + ax3.plot((p0[0], p1[0]), (p0[1], p1[1])) -ax[2].set_title('Probabilistic Hough') -ax[2].axis('image') +ax3.set_title('Probabilistic Hough') +ax3.set_axis_off() +ax3.set_adjustable('box-forced') plt.show() diff --git a/doc/examples/plot_local_equalize.py b/doc/examples/plot_local_equalize.py index cce8fa42..194bfe6a 100644 --- a/doc/examples/plot_local_equalize.py +++ b/doc/examples/plot_local_equalize.py @@ -72,7 +72,14 @@ img_eq = rank.equalize(img, selem=selem) # Display results -fig, axes = plt.subplots(2, 3, figsize=(8, 5)) +fig = plt.figure(figsize=(8, 5)) +axes = np.zeros((2, 3), dtype=np.object) +axes[0,0] = plt.subplot(2, 3, 1, adjustable='box-forced') +axes[0,1] = plt.subplot(2, 3, 2, sharex=axes[0,0], sharey=axes[0,0], adjustable='box-forced') +axes[0,2] = plt.subplot(2, 3, 3, sharex=axes[0,0], sharey=axes[0,0], adjustable='box-forced') +axes[1,0] = plt.subplot(2, 3, 4) +axes[1,1] = plt.subplot(2, 3, 5) +axes[1,2] = plt.subplot(2, 3, 6) ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0]) ax_img.set_title('Low contrast image') diff --git a/doc/examples/plot_local_otsu.py b/doc/examples/plot_local_otsu.py index 210d5c03..d853c25a 100644 --- a/doc/examples/plot_local_otsu.py +++ b/doc/examples/plot_local_otsu.py @@ -37,7 +37,7 @@ threshold_global_otsu = threshold_otsu(img) global_otsu = img >= threshold_global_otsu -fig, ax = plt.subplots(2, 2, figsize=(8, 5)) +fig, ax = plt.subplots(2, 2, figsize=(8, 5), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax1, ax2, ax3, ax4 = ax.ravel() fig.colorbar(ax1.imshow(img, cmap=plt.cm.gray), diff --git a/doc/examples/plot_log_gamma.py b/doc/examples/plot_log_gamma.py index 04a1edbb..70d5881c 100644 --- a/doc/examples/plot_log_gamma.py +++ b/doc/examples/plot_log_gamma.py @@ -54,7 +54,14 @@ gamma_corrected = exposure.adjust_gamma(img, 2) logarithmic_corrected = exposure.adjust_log(img, 1) # Display results -fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(8, 5)) +fig = plt.figure(figsize=(8, 5)) +axes = np.zeros((2,3), dtype=np.object) +axes[0, 0] = plt.subplot(2, 3, 1, adjustable='box-forced') +axes[0, 1] = plt.subplot(2, 3, 2, sharex=axes[0, 0], sharey=axes[0, 0], adjustable='box-forced') +axes[0, 2] = plt.subplot(2, 3, 3, sharex=axes[0, 0], sharey=axes[0, 0], adjustable='box-forced') +axes[1, 0] = plt.subplot(2, 3, 4) +axes[1, 1] = plt.subplot(2, 3, 5) +axes[1, 2] = plt.subplot(2, 3, 6) ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0]) ax_img.set_title('Low contrast image') diff --git a/doc/examples/plot_marked_watershed.py b/doc/examples/plot_marked_watershed.py index 8180f982..3365740c 100644 --- a/doc/examples/plot_marked_watershed.py +++ b/doc/examples/plot_marked_watershed.py @@ -45,7 +45,8 @@ gradient = rank.gradient(denoised, disk(2)) labels = watershed(gradient, markers) # display results -fig, axes = plt.subplots(ncols=4, figsize=(8, 2.7)) +fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(8, 8), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +axes = axes.ravel() ax0, ax1, ax2, ax3 = axes ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest') @@ -61,6 +62,5 @@ ax3.set_title("Segmented") for ax in axes: ax.axis('off') -fig.subplots_adjust(hspace=0.01, wspace=0.01, top=0.9, bottom=0, - left=0, right=1) +fig.tight_layout() plt.show() diff --git a/doc/examples/plot_medial_transform.py b/doc/examples/plot_medial_transform.py index 436e439e..8b9775cd 100644 --- a/doc/examples/plot_medial_transform.py +++ b/doc/examples/plot_medial_transform.py @@ -54,7 +54,7 @@ skel, distance = medial_axis(data, return_distance=True) # Distance to the background for pixels of the skeleton dist_on_skel = distance * skel -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax1.imshow(data, cmap=plt.cm.gray, interpolation='nearest') ax1.axis('off') ax2.imshow(dist_on_skel, cmap=plt.cm.spectral, interpolation='nearest') diff --git a/doc/examples/plot_nonlocal_means.py b/doc/examples/plot_nonlocal_means.py index 207f0ca4..6082ee9c 100644 --- a/doc/examples/plot_nonlocal_means.py +++ b/doc/examples/plot_nonlocal_means.py @@ -26,7 +26,7 @@ noisy = np.clip(noisy, 0, 1) denoise = denoise_nl_means(noisy, 7, 9, 0.08) -fig, ax = plt.subplots(ncols=2, figsize=(8, 4)) +fig, ax = plt.subplots(ncols=2, figsize=(8, 4), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax[0].imshow(noisy) ax[0].axis('off') diff --git a/doc/examples/plot_otsu.py b/doc/examples/plot_otsu.py index 14a7e7cb..c279a65b 100644 --- a/doc/examples/plot_otsu.py +++ b/doc/examples/plot_otsu.py @@ -28,7 +28,12 @@ image = camera() thresh = threshold_otsu(image) binary = image > thresh -fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8, 2.5)) +#fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8, 2.5)) +fig = plt.figure(figsize=(8, 2.5)) +ax1 = plt.subplot(1, 3, 1, adjustable='box-forced') +ax2 = plt.subplot(1, 3, 2) +ax3 = plt.subplot(1, 3, 3, sharex=ax1, sharey=ax1, adjustable='box-forced') + ax1.imshow(image, cmap=plt.cm.gray) ax1.set_title('Original') ax1.axis('off') diff --git a/doc/examples/plot_peak_local_max.py b/doc/examples/plot_peak_local_max.py index aba2aab0..8a690c5e 100644 --- a/doc/examples/plot_peak_local_max.py +++ b/doc/examples/plot_peak_local_max.py @@ -25,7 +25,7 @@ image_max = ndi.maximum_filter(im, size=20, mode='constant') coordinates = peak_local_max(im, min_distance=20) # display results -fig, ax = plt.subplots(1, 3, figsize=(8, 3)) +fig, ax = plt.subplots(1, 3, figsize=(8, 3), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax1, ax2, ax3 = ax.ravel() ax1.imshow(im, cmap=plt.cm.gray) ax1.axis('off') diff --git a/doc/examples/plot_phase_unwrap.py b/doc/examples/plot_phase_unwrap.py index edb061aa..4de93dee 100644 --- a/doc/examples/plot_phase_unwrap.py +++ b/doc/examples/plot_phase_unwrap.py @@ -26,7 +26,7 @@ image_wrapped = np.angle(np.exp(1j * image)) # Perform phase unwrapping image_unwrapped = unwrap_phase(image_wrapped) -fig, ax = plt.subplots(2, 2) +fig, ax = plt.subplots(2, 2, sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax1, ax2, ax3, ax4 = ax.ravel() fig.colorbar(ax1.imshow(image, cmap='gray', vmin=0, vmax=4 * np.pi), ax=ax1) diff --git a/doc/examples/plot_radon_transform.py b/doc/examples/plot_radon_transform.py index 91775ba0..cc327d72 100644 --- a/doc/examples/plot_radon_transform.py +++ b/doc/examples/plot_radon_transform.py @@ -101,7 +101,7 @@ error = reconstruction_fbp - image print('FBP rms reconstruction error: %.3g' % np.sqrt(np.mean(error**2))) imkwargs = dict(vmin=-0.2, vmax=0.2) -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5)) +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax1.set_title("Reconstruction\nFiltered back projection") ax1.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r) ax2.set_title("Reconstruction error\nFiltered back projection") @@ -152,7 +152,7 @@ error = reconstruction_sart - image print('SART (1 iteration) rms reconstruction error: %.3g' % np.sqrt(np.mean(error**2))) -fig, ax = plt.subplots(2, 2, figsize=(8, 8.5)) +fig, ax = plt.subplots(2, 2, figsize=(8, 8.5), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax1, ax2, ax3, ax4 = ax.ravel() ax1.set_title("Reconstruction\nSART") ax1.imshow(reconstruction_sart, cmap=plt.cm.Greys_r) diff --git a/doc/examples/plot_random_walker_segmentation.py b/doc/examples/plot_random_walker_segmentation.py index 446e938f..81e0bd91 100644 --- a/doc/examples/plot_random_walker_segmentation.py +++ b/doc/examples/plot_random_walker_segmentation.py @@ -38,15 +38,18 @@ markers[data > 1.3] = 2 labels = random_walker(data, markers, beta=10, mode='bf') # Plot results -fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8, 3.2)) +fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8, 3.2), sharex=True, sharey=True) ax1.imshow(data, cmap='gray', interpolation='nearest') ax1.axis('off') +ax1.set_adjustable('box-forced') ax1.set_title('Noisy data') ax2.imshow(markers, cmap='hot', interpolation='nearest') ax2.axis('off') +ax2.set_adjustable('box-forced') ax2.set_title('Markers') ax3.imshow(labels, cmap='gray', interpolation='nearest') ax3.axis('off') +ax3.set_adjustable('box-forced') ax3.set_title('Segmentation') fig.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, diff --git a/doc/examples/plot_rank_mean.py b/doc/examples/plot_rank_mean.py index e8f2394c..51f69788 100644 --- a/doc/examples/plot_rank_mean.py +++ b/doc/examples/plot_rank_mean.py @@ -34,19 +34,16 @@ bilateral_result = rank.mean_bilateral(image, selem=selem, s0=500, s1=500) normal_result = rank.mean(image, selem=selem) -fig, axes = plt.subplots(nrows=3, figsize=(8, 10)) -ax0, ax1, ax2 = axes +fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(8, 10), sharex=True, sharey=True) +ax = axes.ravel() -ax0.imshow(np.hstack((image, percentile_result))) -ax0.set_title('Percentile mean') -ax0.axis('off') +titles = ['Original', 'Percentile mean', 'Bilateral mean', 'Local mean'] +imgs = [image, percentile_result, bilateral_result, normal_result] +for n in range(0, len(imgs)): + ax[n].imshow(imgs[n]) + ax[n].set_title(titles[n]) + ax[n].set_adjustable('box-forced') + ax[n].axis('off') -ax1.imshow(np.hstack((image, bilateral_result))) -ax1.set_title('Bilateral mean') -ax1.axis('off') - -ax2.imshow(np.hstack((image, normal_result))) -ax2.set_title('Local mean') -ax2.axis('off') plt.show() diff --git a/doc/examples/plot_ransac3D.py b/doc/examples/plot_ransac3D.py new file mode 100644 index 00000000..644cf154 --- /dev/null +++ b/doc/examples/plot_ransac3D.py @@ -0,0 +1,40 @@ +""" +============================================ +Robust 3D line model estimation using RANSAC +============================================ + +In this example we see how to robustly fit a 3D line model to faulty data using +the RANSAC algorithm. + +""" +import numpy as np +from matplotlib import pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from skimage.measure import LineModelND, ransac + +np.random.seed(seed=1) + +# generate coordinates of line +point = np.array([0, 0, 0], dtype='float') +direction = np.array([1, 1, 1], dtype='float') / np.sqrt(3) +xyz = point + 10 * np.arange(-100, 100)[..., np.newaxis] * direction + +# add gaussian noise to coordinates +noise = np.random.normal(size=xyz.shape) +xyz += 0.5 * noise +xyz[::2] += 20 * noise[::2] +xyz[::4] += 100 * noise[::4] + +# robustly fit line only using inlier data with RANSAC algorithm +model_robust, inliers = ransac(xyz, LineModelND, min_samples=2, + residual_threshold=1, max_trials=1000) +outliers = inliers == False + +fig = plt.figure() +ax = fig.add_subplot(111, projection='3d') +ax.scatter(xyz[inliers][:, 0], xyz[inliers][:, 1], xyz[inliers][:, 2], c='b', + marker='o', label='Inlier data') +ax.scatter(xyz[outliers][:, 0], xyz[outliers][:, 1], xyz[outliers][:, 2], c='r', + marker='o', label='Outlier data') +ax.legend(loc='lower left') +plt.show() diff --git a/doc/examples/plot_regional_maxima.py b/doc/examples/plot_regional_maxima.py index dd676014..a2f9c570 100644 --- a/doc/examples/plot_regional_maxima.py +++ b/doc/examples/plot_regional_maxima.py @@ -36,19 +36,22 @@ Subtracting the dilated image leaves an image with just the coins and a flat, black background, as shown below. """ -fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 2.5)) +fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8, 2.5), sharex=True, sharey=True) ax1.imshow(image) ax1.set_title('original image') ax1.axis('off') +ax1.set_adjustable('box-forced') ax2.imshow(dilated, vmin=image.min(), vmax=image.max()) ax2.set_title('dilated') ax2.axis('off') +ax2.set_adjustable('box-forced') ax3.imshow(image - dilated) ax3.set_title('image - dilated') ax3.axis('off') +ax3.set_adjustable('box-forced') fig.tight_layout() @@ -76,7 +79,7 @@ mask, seed, and dilated images along a slice of the image (indicated by red line). """ -fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 2.5)) +fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8, 2.5)) yslice = 197 diff --git a/doc/examples/plot_register_translation.py b/doc/examples/plot_register_translation.py index 05cf78bf..c558e107 100644 --- a/doc/examples/plot_register_translation.py +++ b/doc/examples/plot_register_translation.py @@ -34,7 +34,10 @@ print(shift) # pixel precision first shift, error, diffphase = register_translation(image, offset_image) -fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 3)) +fig = plt.figure(figsize=(8, 3)) +ax1 = plt.subplot(1, 3, 1, adjustable='box-forced') +ax2 = plt.subplot(1, 3, 2, sharex=ax1, sharey=ax1, adjustable='box-forced') +ax3 = plt.subplot(1, 3, 3) ax1.imshow(image) ax1.set_axis_off() @@ -60,7 +63,10 @@ print(shift) # subpixel precision shift, error, diffphase = register_translation(image, offset_image, 100) -fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 3)) +fig = plt.figure(figsize=(8, 3)) +ax1 = plt.subplot(1, 3, 1, adjustable='box-forced') +ax2 = plt.subplot(1, 3, 2, sharex=ax1, sharey=ax1, adjustable='box-forced') +ax3 = plt.subplot(1, 3, 3) ax1.imshow(image) ax1.set_axis_off() diff --git a/doc/examples/plot_restoration.py b/doc/examples/plot_restoration.py index 52cbec27..221a53b7 100644 --- a/doc/examples/plot_restoration.py +++ b/doc/examples/plot_restoration.py @@ -42,7 +42,7 @@ astro += 0.1 * astro.std() * np.random.standard_normal(astro.shape) deconvolved, _ = restoration.unsupervised_wiener(astro, psf) -fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 5)) +fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 5), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) plt.gray() diff --git a/doc/examples/plot_seam_carving.py b/doc/examples/plot_seam_carving.py index f1dcb948..129b0da0 100644 --- a/doc/examples/plot_seam_carving.py +++ b/doc/examples/plot_seam_carving.py @@ -39,29 +39,28 @@ plt.figure() plt.title('Resized Image') plt.imshow(resized) - """ .. image:: PLOT2RST.current_figure """ out = transform.seam_carve(img, eimg, 'vertical', 200) plt.figure() -plt.title('Resized using Seam-Carving') +plt.title('Resized using Seam Carving') plt.imshow(out) """ .. image:: PLOT2RST.current_figure -As you can see, resizing as distorted the rocket and the objects around, -whereas seam carving has reszied by removing the empty spaces in between. +Resizing distorts the rocket and surrounding objects, whereas seam carving +removes empty spaces and preserves object proportions. Object Removal -------------- -Seam Carving can also be used to remove atrifacts from images. To do that, we -have to ensure that pixels to be removes get less importance. In the following -code I approximately mark the rocket with a mask, and then decrease the -importance of those pixels +Seam carving can also be used to remove artifacts from images. +This requires weighting the artifact with low values. Recall lower weights are +preferentially removed in seam carving. The following code masks the rocket's +region with low weights, indicating it should be removed. """ @@ -78,6 +77,7 @@ plt.figure() plt.title('Object Marked') plt.imshow(masked_img) + """ .. image:: PLOT2RST.current_figure """ @@ -90,6 +90,7 @@ out = transform.seam_carve(img, eimg, 'vertical', 90) resized = transform.resize(img, out.shape) plt.imshow(out) plt.show() + """ .. image:: PLOT2RST.current_figure """ diff --git a/doc/examples/plot_segmentations.py b/doc/examples/plot_segmentations.py index 7f0c5752..af601f53 100644 --- a/doc/examples/plot_segmentations.py +++ b/doc/examples/plot_segmentations.py @@ -79,7 +79,7 @@ print("Felzenszwalb's number of segments: %d" % len(np.unique(segments_fz))) print("Slic number of segments: %d" % len(np.unique(segments_slic))) print("Quickshift number of segments: %d" % len(np.unique(segments_quick))) -fig, ax = plt.subplots(1, 3) +fig, ax = plt.subplots(1, 3, sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) fig.set_size_inches(8, 3, forward=True) fig.subplots_adjust(0.05, 0.05, 0.95, 0.95, 0.05, 0.05) diff --git a/doc/examples/plot_skeleton.py b/doc/examples/plot_skeleton.py index a4a3a60c..2bba3567 100644 --- a/doc/examples/plot_skeleton.py +++ b/doc/examples/plot_skeleton.py @@ -47,7 +47,7 @@ image[circle2] = 0 skeleton = skeletonize(image) # display results -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5)) +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax1.imshow(image, cmap=plt.cm.gray) ax1.axis('off') diff --git a/doc/examples/plot_ssim.py b/doc/examples/plot_ssim.py index 1592136b..112a6971 100644 --- a/doc/examples/plot_ssim.py +++ b/doc/examples/plot_ssim.py @@ -45,7 +45,7 @@ def mse(x, y): img_noise = img + noise img_const = img + abs(noise) -fig, (ax0, ax1, ax2) = plt.subplots(nrows=1, ncols=3, figsize=(8, 4)) +fig, (ax0, ax1, ax2) = plt.subplots(nrows=1, ncols=3, figsize=(8, 4), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) mse_none = mse(img, img) ssim_none = ssim(img, img, dynamic_range=img.max() - img.min()) diff --git a/doc/examples/plot_swirl.py b/doc/examples/plot_swirl.py index 8f79db94..b93efe9f 100644 --- a/doc/examples/plot_swirl.py +++ b/doc/examples/plot_swirl.py @@ -74,7 +74,7 @@ from skimage.transform import swirl image = data.checkerboard() swirled = swirl(image, rotation=0, strength=10, radius=120, order=2) -fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 3)) +fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 3), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax0.imshow(image, cmap=plt.cm.gray, interpolation='none') ax0.axis('off') diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py index 08581625..c7b4f01a 100644 --- a/doc/examples/plot_template.py +++ b/doc/examples/plot_template.py @@ -33,7 +33,10 @@ result = match_template(image, coin) ij = np.unravel_index(np.argmax(result), result.shape) x, y = ij[::-1] -fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 3)) +fig = plt.figure(figsize=(8, 3)) +ax1 = plt.subplot(1, 3, 1) +ax2 = plt.subplot(1, 3, 2, adjustable='box-forced') +ax3 = plt.subplot(1, 3, 3, sharex=ax2, sharey=ax2, adjustable='box-forced') ax1.imshow(coin) ax1.set_axis_off() diff --git a/doc/examples/plot_tinting_grayscale_images.py b/doc/examples/plot_tinting_grayscale_images.py index 008ed065..2ba54396 100644 --- a/doc/examples/plot_tinting_grayscale_images.py +++ b/doc/examples/plot_tinting_grayscale_images.py @@ -28,9 +28,11 @@ image = color.gray2rgb(grayscale_image) red_multiplier = [1, 0, 0] yellow_multiplier = [1, 1, 0] -fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4)) +fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4), sharex=True, sharey=True) ax1.imshow(red_multiplier * image) ax2.imshow(yellow_multiplier * image) +ax1.set_adjustable('box-forced') +ax2.set_adjustable('box-forced') """ .. image:: PLOT2RST.current_figure @@ -104,13 +106,14 @@ and a non-zero saturation: hue_rotations = np.linspace(0, 1, 6) -fig, axes = plt.subplots(nrows=2, ncols=3) +fig, axes = plt.subplots(nrows=2, ncols=3, sharex=True, sharey=True) for ax, hue in zip(axes.flat, hue_rotations): # Turn down the saturation to give it that vintage look. tinted_image = colorize(image, hue, saturation=0.3) ax.imshow(tinted_image, vmin=0, vmax=1) ax.set_axis_off() + ax.set_adjustable('box-forced') fig.tight_layout() """ @@ -142,9 +145,11 @@ textured_regions = noisy > 4 masked_image = image.copy() masked_image[textured_regions, :] *= red_multiplier -fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4)) +fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(8, 4), sharex=True, sharey=True) ax1.imshow(sliced_image) ax2.imshow(masked_image) +ax1.set_adjustable('box-forced') +ax2.set_adjustable('box-forced') plt.show() diff --git a/doc/examples/plot_view_as_blocks.py b/doc/examples/plot_view_as_blocks.py index eeb510ed..1a85bc41 100644 --- a/doc/examples/plot_view_as_blocks.py +++ b/doc/examples/plot_view_as_blocks.py @@ -44,21 +44,26 @@ max_view = np.max(flatten_view, axis=2) median_view = np.median(flatten_view, axis=2) # -- display resampled images -fig, axes = plt.subplots(2, 2, figsize=(8, 8)) +fig, axes = plt.subplots(2, 2, figsize=(8, 8), sharex=True, sharey=True) ax0, ax1, ax2, ax3 = axes.ravel() ax0.set_title("Original rescaled with\n spline interpolation (order=3)") l_resized = ndi.zoom(l, 2, order=3) -ax0.imshow(l_resized, cmap=cm.Greys_r) +#ax0.imshow(l_resized, cmap=cm.Greys_r) +ax0.imshow(l_resized, extent=(0, 128, 128, 0), interpolation='nearest', cmap=cm.Greys_r) +ax0.set_axis_off() ax1.set_title("Block view with\n local mean pooling") -ax1.imshow(mean_view, cmap=cm.Greys_r) +ax1.imshow(mean_view, interpolation='nearest', cmap=cm.Greys_r) +ax1.set_axis_off() ax2.set_title("Block view with\n local max pooling") -ax2.imshow(max_view, cmap=cm.Greys_r) +ax2.imshow(max_view, interpolation='nearest', cmap=cm.Greys_r) +ax2.set_axis_off() ax3.set_title("Block view with\n local median pooling") -ax3.imshow(median_view, cmap=cm.Greys_r) +ax3.imshow(median_view, interpolation='nearest', cmap=cm.Greys_r) +ax3.set_axis_off() fig.subplots_adjust(hspace=0.4, wspace=0.4) plt.show() diff --git a/doc/examples/plot_watershed.py b/doc/examples/plot_watershed.py index 5d65238f..664c4c77 100644 --- a/doc/examples/plot_watershed.py +++ b/doc/examples/plot_watershed.py @@ -48,7 +48,7 @@ local_maxi = peak_local_max(distance, indices=False, footprint=np.ones((3, 3)), markers = ndi.label(local_maxi)[0] labels = watershed(-distance, markers, mask=image) -fig, axes = plt.subplots(ncols=3, figsize=(8, 2.7)) +fig, axes = plt.subplots(ncols=3, figsize=(8, 2.7), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) ax0, ax1, ax2 = axes ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest') diff --git a/setup.py b/setup.py index 6c55c6d6..53b5aa71 100644 --- a/setup.py +++ b/setup.py @@ -87,14 +87,15 @@ if __name__ == "__main__": except ImportError: if len(sys.argv) >= 2 and ('--help' in sys.argv[1:] or sys.argv[1] in ('--help-commands', - 'egg_info', '--version', + '--version', 'clean')): # For these actions, NumPy is not required. # # They are required to succeed without Numpy for example when # pip is used to install scikit-image when Numpy is not yet # present in the system. - pass + from setuptools import setup + extra = {} else: print('To install scikit-image from source, you will need numpy.\n' + 'Install numpy with pip:\n' + diff --git a/skimage/color/colorconv.py b/skimage/color/colorconv.py index 86353e0b..137c14c0 100644 --- a/skimage/color/colorconv.py +++ b/skimage/color/colorconv.py @@ -768,8 +768,6 @@ def gray2rgb(image, alpha=None): else: return np.concatenate(3 * (image,), axis=-1) - return image - else: raise ValueError("Input image expected to be RGB, RGBA or gray.") diff --git a/skimage/color/tests/test_colorconv.py b/skimage/color/tests/test_colorconv.py index 71e2d849..4c991416 100644 --- a/skimage/color/tests/test_colorconv.py +++ b/skimage/color/tests/test_colorconv.py @@ -439,6 +439,8 @@ def test_gray2rgb(): assert_equal(y.shape, (3, 1, 3)) assert_equal(y.dtype, x.dtype) + assert_equal(y[..., 0], x) + assert_equal(y[0, 0, :], [0, 0, 0]) x = np.array([[0, 128, 255]], dtype=np.uint8) z = gray2rgb(x) diff --git a/skimage/data/rank_filter_tests.npz b/skimage/data/rank_filter_tests.npz index 3142fcc1..d5169bc5 100644 Binary files a/skimage/data/rank_filter_tests.npz and b/skimage/data/rank_filter_tests.npz differ diff --git a/skimage/exposure/_adapthist.py b/skimage/exposure/_adapthist.py index 757dafff..62ce2ecd 100644 --- a/skimage/exposure/_adapthist.py +++ b/skimage/exposure/_adapthist.py @@ -244,6 +244,8 @@ def clip_histogram(hist, clip_limit): n_excess -= mid.size * clip_limit - mid.sum() hist[mid_mask] = clip_limit + prev_n_excess = n_excess + while n_excess > 0: # Redistribute remaining excess index = 0 while n_excess > 0 and index < hist.size: @@ -256,6 +258,10 @@ def clip_histogram(hist, clip_limit): hist[under_mask] += 1 n_excess -= under_mask.sum() index += 1 + # bail if we have not distributed any excess + if prev_n_excess == n_excess: + break + prev_n_excess = n_excess return hist diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py index 1b9c778a..4b88b67f 100644 --- a/skimage/exposure/tests/test_exposure.py +++ b/skimage/exposure/tests/test_exposure.py @@ -213,7 +213,7 @@ def test_adapthist_grayscale(): img = np.dstack((img, img, img)) with expected_warnings(['precision loss|non-contiguous input', 'deprecated']): - adapted_old = exposure.equalize_adapthist(img, 10, 9, clip_limit=0.01, + adapted_old = exposure.equalize_adapthist(img, 10, 9, clip_limit=0.001, nbins=128) adapted = exposure.equalize_adapthist(img, kernel_size=(57, 51), clip_limit=0.01, nbins=128) assert img.shape == adapted.shape diff --git a/skimage/feature/texture.py b/skimage/feature/texture.py index c26eb59d..4be6c749 100644 --- a/skimage/feature/texture.py +++ b/skimage/feature/texture.py @@ -298,7 +298,7 @@ def local_binary_pattern(image, P, R, method='default'): def multiblock_lbp(int_image, r, c, width, height): - """Multi-block local binary pattern (MB-LBP) [1]_. + """Multi-block local binary pattern (MB-LBP). The features are calculated similarly to local binary patterns (LBPs), (See :py:meth:`local_binary_pattern`) except that summed blocks are diff --git a/skimage/filters/__init__.py b/skimage/filters/__init__.py index 9dd0b4df..708f5d53 100644 --- a/skimage/filters/__init__.py +++ b/skimage/filters/__init__.py @@ -5,7 +5,7 @@ from .edges import (sobel, hsobel, vsobel, sobel_h, sobel_v, prewitt, hprewitt, vprewitt, prewitt_h, prewitt_v, roberts, roberts_positive_diagonal, roberts_negative_diagonal, roberts_pos_diag, - roberts_neg_diag) + roberts_neg_diag, laplace) from ._rank_order import rank_order from ._gabor import gabor_kernel, gabor from .thresholding import (threshold_adaptive, threshold_otsu, threshold_yen, @@ -62,6 +62,7 @@ __all__ = ['inverse', 'roberts_negative_diagonal', 'roberts_pos_diag', 'roberts_neg_diag', + 'laplace', 'denoise_tv_chambolle', 'denoise_bilateral', 'denoise_tv_bregman', diff --git a/skimage/filters/_gabor.py b/skimage/filters/_gabor.py index 211571f9..8bb7b61b 100644 --- a/skimage/filters/_gabor.py +++ b/skimage/filters/_gabor.py @@ -56,7 +56,7 @@ def gabor_kernel(frequency, theta=0, bandwidth=1, sigma_x=None, sigma_y=None, Examples -------- - >>> from skimage.filter import gabor_kernel + >>> from skimage.filters import gabor_kernel >>> from skimage import io >>> from matplotlib import pyplot as plt # doctest: +SKIP @@ -148,7 +148,7 @@ def gabor(image, frequency, theta=0, bandwidth=1, sigma_x=None, Examples -------- - >>> from skimage.filter import gabor + >>> from skimage.filters import gabor >>> from skimage import data, io >>> from matplotlib import pyplot as plt # doctest: +SKIP diff --git a/skimage/filters/edges.py b/skimage/filters/edges.py index 1eeed4de..4cf083a3 100644 --- a/skimage/filters/edges.py +++ b/skimage/filters/edges.py @@ -14,6 +14,7 @@ from .. import img_as_float from .._shared.utils import assert_nD, deprecated from scipy.ndimage import convolve, binary_erosion, generate_binary_structure +from ..restoration.uft import laplacian EROSION_SELEM = generate_binary_structure(2, 2) @@ -175,6 +176,7 @@ def hsobel(image, mask=None): Parameters ---------- + image : 2-D array Image to process. mask : 2-D array, optional @@ -762,3 +764,35 @@ def roberts_negative_diagonal(image, mask=None): """ return np.abs(roberts_neg_diag(image, mask)) + +def laplace(image, ksize=3, mask=None): + """Find the edges of an image using the Laplace operator. + + Parameters + ---------- + image : ndarray + Image to process. + ksize : int, optional + Define the size of the discrete Laplacian operator such that it + will have a size of (ksize,) * image.ndim. + mask : ndarray, optional + An optional mask to limit the application to a certain area. + Note that pixels surrounding masked regions are also masked to + prevent masked regions from affecting the result. + + Returns + ------- + output : ndarray + The Laplace edge map. + + Notes + ----- + The Laplacian operator is generated using the function + skimage.restoration.uft.laplacian(). + + """ + image = img_as_float(image) + # Create the discrete Laplacian operator - We keep only the real part of the filter + _, laplace_op = laplacian(image.ndim, (ksize, ) * image.ndim) + result = convolve(image, laplace_op) + return _mask_filter_result(result, mask) diff --git a/skimage/filters/rank/__init__.py b/skimage/filters/rank/__init__.py index 37b5f9cd..948adee6 100644 --- a/skimage/filters/rank/__init__.py +++ b/skimage/filters/rank/__init__.py @@ -1,7 +1,7 @@ from .generic import (autolevel, bottomhat, equalize, gradient, maximum, mean, - subtract_mean, median, minimum, modal, enhance_contrast, - pop, threshold, tophat, noise_filter, entropy, otsu, - sum, windowed_histogram) + geometric_mean, subtract_mean, median, minimum, modal, + enhance_contrast, pop, threshold, tophat, noise_filter, + entropy, otsu, sum, windowed_histogram) from ._percentile import (autolevel_percentile, gradient_percentile, mean_percentile, subtract_mean_percentile, enhance_contrast_percentile, percentile, @@ -17,6 +17,7 @@ __all__ = ['autolevel', 'gradient_percentile', 'maximum', 'mean', + 'geometric_mean', 'mean_percentile', 'mean_bilateral', 'subtract_mean', diff --git a/skimage/filters/rank/generic.py b/skimage/filters/rank/generic.py index fd718b98..5c6a352b 100644 --- a/skimage/filters/rank/generic.py +++ b/skimage/filters/rank/generic.py @@ -25,8 +25,9 @@ from . import generic_cy __all__ = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', 'mean', - 'subtract_mean', 'median', 'minimum', 'modal', 'enhance_contrast', - 'pop', 'threshold', 'tophat', 'noise_filter', 'entropy', 'otsu'] + 'geometric_mean', 'subtract_mean', 'median', 'minimum', 'modal', + 'enhance_contrast', 'pop', 'threshold', 'tophat', 'noise_filter', + 'entropy', 'otsu'] def _handle_input(image, selem, out, mask, out_dtype=None, pixel_size=1): @@ -342,6 +343,48 @@ def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False): return _apply_scalar_per_pixel(generic_cy._mean, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) +def geometric_mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Return local geometric mean of an image. + + Parameters + ---------- + image : 2-D array (uint8, uint16) + Input image. + selem : 2-D array + The neighborhood expressed as a 2-D array of 1's and 0's. + out : 2-D array (same dtype as input) + If None, a new array is allocated. + mask : ndarray + Mask array that defines (>0) area of the image included in the local + neighborhood. If None, the complete image is used (default). + shift_x, shift_y : int + Offset added to the structuring element center point. Shift is bounded + to the structuring element sizes (center must be inside the given + structuring element). + + Returns + ------- + out : 2-D array (same dtype as input image) + Output image. + + Examples + -------- + >>> from skimage import data + >>> from skimage.morphology import disk + >>> from skimage.filters.rank import mean + >>> img = data.camera() + >>> avg = geometric_mean(img, disk(5)) + + References + ---------- + .. [1] Gonzalez, R. C. and Wood, R. E. "Digital Image Processing (3rd Edition)." + Prentice-Hall Inc, 2006. + + """ + + return _apply_scalar_per_pixel(generic_cy._geometric_mean, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + def subtract_mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False): diff --git a/skimage/filters/rank/generic_cy.pyx b/skimage/filters/rank/generic_cy.pyx index 1356f369..11e6ba8c 100644 --- a/skimage/filters/rank/generic_cy.pyx +++ b/skimage/filters/rank/generic_cy.pyx @@ -4,10 +4,11 @@ #cython: wraparound=False cimport numpy as cnp -from libc.math cimport log +from libc.math cimport log, exp from .core_cy cimport dtype_t, dtype_t_out, _core +from ..._shared.interpolation cimport round cdef inline void _kernel_autolevel(dtype_t_out* out, Py_ssize_t odepth, Py_ssize_t* histo, @@ -133,6 +134,25 @@ cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth, out[0] = 0 +cdef inline void _kernel_geometric_mean(dtype_t_out* out, Py_ssize_t odepth, + Py_ssize_t* histo, + double pop, dtype_t g, + Py_ssize_t max_bin, Py_ssize_t mid_bin, + double p0, double p1, + Py_ssize_t s0, Py_ssize_t s1) nogil: + + cdef Py_ssize_t i + cdef double mean = 0. + + if pop: + for i in range(max_bin): + if histo[i]: + mean += (histo[i] * log(i+1)) + out[0] = round(exp(mean / pop)-1) + else: + out[0] = 0 + + cdef inline void _kernel_subtract_mean(dtype_t_out* out, Py_ssize_t odepth, Py_ssize_t* histo, double pop, dtype_t g, @@ -467,6 +487,16 @@ def _mean(dtype_t[:, ::1] image, shift_x, shift_y, 0, 0, 0, 0, max_bin) +def _geometric_mean(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, :, ::1] out, + signed char shift_x, signed char shift_y, Py_ssize_t max_bin): + + _core(_kernel_geometric_mean[dtype_t_out, dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + def _subtract_mean(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, diff --git a/skimage/filters/rank/tests/test_rank.py b/skimage/filters/rank/tests/test_rank.py index 58e338f9..3fcc79ac 100644 --- a/skimage/filters/rank/tests/test_rank.py +++ b/skimage/filters/rank/tests/test_rank.py @@ -39,6 +39,8 @@ def check_all(): rank.maximum(image, selem)) assert_equal(refs["mean"], rank.mean(image, selem)) + assert_equal(refs["geometric_mean"], + rank.geometric_mean(image, selem)), assert_equal(refs["mean_percentile"], rank.mean_percentile(image, selem)) assert_equal(refs["mean_bilateral"], @@ -102,6 +104,13 @@ def test_random_sizes(): rank.mean(image=image8, selem=elem, mask=mask, out=out8, shift_x=+1, shift_y=+1) assert_equal(image8.shape, out8.shape) + + rank.geometric_mean(image=image8, selem=elem, mask=mask, out=out8, + shift_x=0, shift_y=0) + assert_equal(image8.shape, out8.shape) + rank.geometric_mean(image=image8, selem=elem, mask=mask, out=out8, + shift_x=+1, shift_y=+1) + assert_equal(image8.shape, out8.shape) image16 = np.ones((m, n), dtype=np.uint16) out16 = np.empty_like(image8, dtype=np.uint16) @@ -112,6 +121,13 @@ def test_random_sizes(): shift_x=+1, shift_y=+1) assert_equal(image16.shape, out16.shape) + rank.geometric_mean(image=image16, selem=elem, mask=mask, out=out16, + shift_x=0, shift_y=0) + assert_equal(image16.shape, out16.shape) + rank.geometric_mean(image=image16, selem=elem, mask=mask, out=out16, + shift_x=+1, shift_y=+1) + assert_equal(image16.shape, out16.shape) + rank.mean_percentile(image=image16, mask=mask, out=out16, selem=elem, shift_x=0, shift_y=0, p0=.1, p1=.9) assert_equal(image16.shape, out16.shape) @@ -292,8 +308,8 @@ def test_compare_8bit_unsigned_vs_signed(): assert_equal(image_u, img_as_ubyte(image_s)) methods = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', - 'mean', 'subtract_mean', 'median', 'minimum', 'modal', - 'enhance_contrast', 'pop', 'threshold', 'tophat'] + 'mean', 'geometric_mean', 'subtract_mean', 'median', 'minimum', + 'modal', 'enhance_contrast', 'pop', 'threshold', 'tophat'] for method in methods: func = getattr(rank, method) @@ -338,6 +354,9 @@ def test_trivial_selem8(): rank.mean(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) + rank.geometric_mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_equal(image, out) rank.minimum(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) @@ -361,6 +380,9 @@ def test_trivial_selem16(): rank.mean(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) + rank.geometric_mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_equal(image, out) rank.minimum(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) @@ -407,6 +429,9 @@ def test_smallest_selem16(): rank.mean(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) + rank.geometric_mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_equal(image, out) rank.minimum(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) @@ -432,6 +457,9 @@ def test_empty_selem(): rank.mean(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(res, out) + rank.geometric_mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_equal(res, out) rank.minimum(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(res, out) @@ -514,6 +542,9 @@ def test_selem_dtypes(): rank.mean(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) + rank.geometric_mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_equal(image, out) rank.mean_percentile(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) diff --git a/skimage/filters/tests/test_edges.py b/skimage/filters/tests/test_edges.py index d2c0f469..0fbaa4ac 100644 --- a/skimage/filters/tests/test_edges.py +++ b/skimage/filters/tests/test_edges.py @@ -335,6 +335,34 @@ def test_vprewitt_horizontal(): assert_allclose(result, 0) +def test_laplace_zeros(): + """Laplace on a square image.""" + # Create a synthetic 2D image + image = np.zeros((9, 9)) + image[3:-3, 3:-3] = 1 + result = filters.laplace(image) + res_chk = np.array([[ 0., 0., 0., 0., 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0., 0., 0., 0., 0.], + [ 0., 0., 0., -1., -1., -1., 0., 0., 0.], + [ 0., 0., -1., 2., 1., 2., -1., 0., 0.], + [ 0., 0., -1., 1., 0., 1., -1., 0., 0.], + [ 0., 0., -1., 2., 1., 2., -1., 0., 0.], + [ 0., 0., 0., -1., -1., -1., 0., 0., 0.], + [ 0., 0., 0., 0., 0., 0., 0., 0., 0.], + [ 0., 0., 0., 0., 0., 0., 0., 0., 0.]]) + assert_allclose(result, res_chk) + + +def test_laplace_mask(): + """Laplace on a masked array should be zero.""" + # Create a synthetic 2D image + image = np.zeros((9, 9)) + image[3:-3, 3:-3] = 1 + # Define the mask + result = filters.laplace(image, ksize=3, mask=np.zeros((9, 9), bool)) + assert (np.all(result == 0)) + + def test_horizontal_mask_line(): """Horizontal edge filters mask pixels surrounding input mask.""" vgrad, _ = np.mgrid[:1:11j, :1:11j] # vertical gradient with spacing 0.1 @@ -351,7 +379,6 @@ def test_horizontal_mask_line(): result = grad_func(vgrad, mask) yield assert_close, result, expected - def test_vertical_mask_line(): """Vertical edge filters mask pixels surrounding input mask.""" _, hgrad = np.mgrid[:1:11j, :1:11j] # horizontal gradient with spacing 0.1 diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py index 9731d6da..9e8df953 100755 --- a/skimage/measure/__init__.py +++ b/skimage/measure/__init__.py @@ -7,7 +7,7 @@ from ._polygon import approximate_polygon, subdivide_polygon from ._pnpoly import points_in_poly, grid_points_in_poly from ._moments import moments, moments_central, moments_normalized, moments_hu from .profile import profile_line -from .fit import LineModel, CircleModel, EllipseModel, ransac +from .fit import LineModel, LineModelND, CircleModel, EllipseModel, ransac from .block import block_reduce from ._label import label @@ -19,6 +19,7 @@ __all__ = ['find_contours', 'approximate_polygon', 'subdivide_polygon', 'LineModel', + 'LineModelND', 'CircleModel', 'EllipseModel', 'ransac', diff --git a/skimage/measure/_label.py b/skimage/measure/_label.py index e3ac6034..9dcefed4 100644 --- a/skimage/measure/_label.py +++ b/skimage/measure/_label.py @@ -1,7 +1,7 @@ from ._ccomp import label as _label def label(input, neighbors=None, background=None, return_num=False, - connectivity=None): + connectivity=None): return _label(input, neighbors, background, return_num, connectivity) label.__doc__ = _label.__doc__ diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py index a334fbbd..96fb1102 100644 --- a/skimage/measure/_marching_cubes.py +++ b/skimage/measure/_marching_cubes.py @@ -228,10 +228,14 @@ def correct_mesh_orientation(volume, verts, faces, spacing=(1., 1., 1.), import scipy.ndimage as ndi # Calculate gradient of `volume`, then interpolate to vertices in `verts` - grad_x, grad_y, grad_z = np.gradient(volume, *spacing) + grad_x, grad_y, grad_z = np.gradient(volume) # Fancy indexing to define two vector arrays from triangle vertices actual_verts = verts[faces] + actual_verts[:, 0] /= spacing[0] + actual_verts[:, 1] /= spacing[1] + actual_verts[:, 2] /= spacing[2] + a = actual_verts[:, 0, :] - actual_verts[:, 1, :] b = actual_verts[:, 0, :] - actual_verts[:, 2, :] diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index 1655ec61..f38fda22 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -157,7 +157,7 @@ class _RegionProperties(object): @property def euler_number(self): euler_array = self.filled_image != self.image - _, num = label(euler_array, neighbors=8, return_num=True) + _, num = label(euler_array, neighbors=8, return_num=True, background=-1) return -num + 1 @property @@ -473,7 +473,7 @@ def regionprops(label_image, intensity_image=None, cache=True): Examples -------- >>> from skimage import data, util - >>> from skimage.morphology import label + >>> from skimage.measure import label >>> img = util.img_as_ubyte(data.coins()) > 110 >>> label_img = label(img, connectivity=img.ndim) >>> props = regionprops(label_img) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index f8e65117..b227d59f 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -2,6 +2,7 @@ import math import warnings import numpy as np from scipy import optimize +from .._shared.utils import skimage_deprecation def _check_data_dim(data, dim): @@ -9,6 +10,16 @@ def _check_data_dim(data, dim): raise ValueError('Input data must have shape (N, %d).' % dim) +def _check_data_atleast_2D(data): + if data.ndim < 2 or data.shape[1] < 2: + raise ValueError('Input data must be at least 2D.') + + +def _norm_along_axis(x, axis): + """NumPy < 1.8 does not support the `axis` argument for `np.linalg.norm`.""" + return np.sqrt(np.einsum('ij,ij->i', x, x)) + + class BaseModel(object): def __init__(self): @@ -39,6 +50,8 @@ class LineModel(BaseModel): A minimum number of 2 points is required to solve for the parameters. + **Deprecated class**. Use ``LineModelND`` instead. + Attributes ---------- params : tuple @@ -46,6 +59,11 @@ class LineModel(BaseModel): """ + def __init__(self): + self.params = None + warnings.warn(skimage_deprecation('`LineModel` is deprecated, ' + 'use `LineModelND` instead.')) + def estimate(self, data): """Estimate line model from data using total least squares. @@ -156,6 +174,158 @@ class LineModel(BaseModel): return (dist - x * math.cos(theta)) / math.sin(theta) +class LineModelND(BaseModel): + """Total least squares estimator for N-dimensional lines. + + Lines are defined by a point (origin) and a unit vector (direction) + according to the following vector equation:: + + X = origin + lambda * direction + + Attributes + ---------- + params : tuple + Line model parameters in the following order `origin`, `direction`. + + """ + + def estimate(self, data): + """Estimate line model from data. + + Parameters + ---------- + data : (N, dim) array + N points in a space of dimensionality dim >= 2. + + Returns + ------- + success : bool + True, if model estimation succeeds. + """ + + _check_data_atleast_2D(data) + + X0 = data.mean(axis=0) + + if data.shape[0] == 2: # well determined + u = data[1] - data[0] + norm = np.linalg.norm(u) + if norm > 0: + u /= norm + elif data.shape[0] > 2: # over-determined + data = data - X0 + # first principal component + # Note: without full_matrices=False Python dies with joblib + # parallel_for. + _, _, u = np.linalg.svd(data, full_matrices=False) + u = u[0] + else: # under-determined + raise ValueError('At least 2 input points needed.') + + self.params = (X0, u) + + return True + + def residuals(self, data): + """Determine residuals of data to model. + + For each point the shortest distance to the line is returned. + It is obtained by projecting the data onto the line. + + Parameters + ---------- + data : (N, dim) array + N points in a space of dimension dim. + + Returns + ------- + residuals : (N, ) array + Residual for each data point. + """ + + X0, u = self.params + return _norm_along_axis((data - X0) - + np.dot(data - X0, u)[..., np.newaxis] * u, + axis=1) + + def predict(self, x, axis=0, params=None): + """Predict intersection of the estimated line model with a hyperplane + orthogonal to a given axis. + + Parameters + ---------- + x : array + coordinates along an axis. + axis : int + axis orthogonal to the hyperplane intersecting the line. + params : (2, ) array, optional + Optional custom parameter set in the form (`origin`, `direction`). + + Returns + ------- + y : array + Predicted coordinates. + + If the line is parallel to the given axis, a ValueError is raised. + """ + + if params is None: + params = self.params + + X0, u = params + + if u[axis] == 0: + # line parallel to axis + raise ValueError('Line parallel to axis %s' % axis) + + l = (x - X0[axis]) / u[axis] + return X0 + l[..., np.newaxis] * u + + def predict_x(self, y, params=None): + """Predict x-coordinates for 2D lines using the estimated model. + + Alias for:: + + predict(y, axis=1)[:, 0] + + Parameters + ---------- + y : array + y-coordinates. + params : (2, ) array, optional + Optional custom parameter set in the form (`origin`, `direction`). + + Returns + ------- + x : array + Predicted x-coordinates. + + """ + return self.predict(y, axis=1, params=params)[:, 0] + + def predict_y(self, x, params=None): + """Predict y-coordinates for 2D lines using the estimated model. + + Alias for:: + + predict(x, axis=0)[:, 1] + + Parameters + ---------- + x : array + x-coordinates. + params : (2, ) array, optional + Optional custom parameter set in the form (`origin`, `direction`). + + Returns + ------- + y : array + Predicted y-coordinates. + + """ + return self.predict(x, axis=0, params=params)[:, 1] + + class CircleModel(BaseModel): """Total least squares estimator for 2D circles. diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py index 7f98c971..794d173e 100644 --- a/skimage/measure/tests/test_fit.py +++ b/skimage/measure/tests/test_fit.py @@ -1,18 +1,18 @@ import numpy as np from numpy.testing import assert_equal, assert_raises, assert_almost_equal -from skimage.measure import LineModel, CircleModel, EllipseModel, ransac +from skimage.measure import LineModelND, CircleModel, EllipseModel, ransac from skimage.transform import AffineTransform from skimage.measure.fit import _dynamic_max_trials from skimage._shared._warnings import expected_warnings def test_line_model_invalid_input(): - assert_raises(ValueError, LineModel().estimate, np.empty((5, 3))) + assert_raises(ValueError, LineModelND().estimate, np.empty((1, 3))) def test_line_model_predict(): - model = LineModel() - model.params = (10, 1) + model = LineModelND() + model.params = ((0, 0), (1, 1)) x = np.arange(-10, 10) y = model.predict_y(x) assert_almost_equal(x, model.predict_x(y)) @@ -20,38 +20,92 @@ def test_line_model_predict(): def test_line_model_estimate(): # generate original data without noise - model0 = LineModel() - model0.params = (10, 1) + model0 = LineModelND() + model0.params = ((0, 0), (1, 1)) x0 = np.arange(-100, 100) y0 = model0.predict_y(x0) - data0 = np.column_stack([x0, y0]) + + data = np.column_stack([x0, y0]) + + # estimate parameters of noisy data + model_est = LineModelND() + model_est.estimate(data) + + # test whether estimated parameters almost equal original parameters + x = np.random.rand(100, 2) + assert_almost_equal(model0.predict(x), model_est.predict(x), 1) + + +def test_line_model_residuals(): + model = LineModelND() + model.params = (np.array([0, 0]), np.array([0, 1])) + assert_equal(model.residuals(np.array([[0, 0]])), 0) + assert_equal(model.residuals(np.array([[0, 10]])), 0) + assert_equal(model.residuals(np.array([[10, 0]])), 10) + model.params = (np.array([-2, 0]), np.array([1, 1]) / np.sqrt(2)) + assert_equal(model.residuals(np.array([[0, 0]])), np.sqrt(2)) + assert_almost_equal(model.residuals(np.array([[-4, 0]])), np.sqrt(2)) + + +def test_line_model_under_determined(): + data = np.empty((1, 2)) + assert_raises(ValueError, LineModelND().estimate, data) + + +def test_line_modelND_invalid_input(): + assert_raises(ValueError, LineModelND().estimate, np.empty((5, 1))) + + +def test_line_modelND_predict(): + model = LineModelND() + model.params = (np.array([0, 0]), np.array([0.2, 0.98])) + x = np.arange(-10, 10) + y = model.predict_y(x) + assert_almost_equal(x, model.predict_x(y)) + + +def test_line_modelND_estimate(): + # generate original data without noise + model0 = LineModelND() + model0.params = (np.array([0,0,0], dtype='float'), + np.array([1,1,1], dtype='float')/np.sqrt(3)) + # we scale the unit vector with a factor 10 when generating points on the + # line in order to compensate for the scale of the random noise + data0 = (model0.params[0] + + 10 * np.arange(-100,100)[...,np.newaxis] * model0.params[1]) # add gaussian noise to data np.random.seed(1234) data = data0 + np.random.normal(size=data0.shape) # estimate parameters of noisy data - model_est = LineModel() + model_est = LineModelND() model_est.estimate(data) - # test whether estimated parameters almost equal original parameters - assert_almost_equal(model0.params, model_est.params, 1) + # test whether estimated parameters are correct + # we use the following geometric property: two aligned vectors have + # a cross-product equal to zero + # test if direction vectors are aligned + assert_almost_equal(np.linalg.norm(np.cross(model0.params[1], + model_est.params[1])), 0, 1) + # test if origins are aligned with the direction + a = model_est.params[0] - model0.params[0] + if np.linalg.norm(a) > 0: + a /= np.linalg.norm(a) + assert_almost_equal(np.linalg.norm(np.cross(model0.params[1], a)), 0, 1) -def test_line_model_residuals(): - model = LineModel() - model.params = (0, 0) - assert_equal(abs(model.residuals(np.array([[0, 0]]))), 0) - assert_equal(abs(model.residuals(np.array([[0, 10]]))), 0) - assert_equal(abs(model.residuals(np.array([[10, 0]]))), 10) - model.params = (5, np.pi / 4) - assert_equal(abs(model.residuals(np.array([[0, 0]]))), 5) - assert_almost_equal(abs(model.residuals(np.array([[np.sqrt(50), 0]]))), 0) +def test_line_modelND_residuals(): + model = LineModelND() + model.params = (np.array([0, 0, 0]), np.array([0, 0, 1])) + assert_equal(abs(model.residuals(np.array([[0, 0, 0]]))), 0) + assert_equal(abs(model.residuals(np.array([[0, 0, 1]]))), 0) + assert_equal(abs(model.residuals(np.array([[10, 0, 0]]))), 10) -def test_line_model_under_determined(): - data = np.empty((1, 2)) - assert_raises(ValueError, LineModel().estimate, data) +def test_line_modelND_under_determined(): + data = np.empty((1, 3)) + assert_raises(ValueError, LineModelND().estimate, data) def test_circle_model_invalid_input(): @@ -189,7 +243,7 @@ def test_ransac_is_data_valid(): np.random.seed(1) is_data_valid = lambda data: data.shape[0] > 2 - model, inliers = ransac(np.empty((10, 2)), LineModel, 2, np.inf, + model, inliers = ransac(np.empty((10, 2)), LineModelND, 2, np.inf, is_data_valid=is_data_valid) assert_equal(model, None) assert_equal(inliers, None) @@ -200,7 +254,7 @@ def test_ransac_is_model_valid(): def is_model_valid(model, data): return False - model, inliers = ransac(np.empty((10, 2)), LineModel, 2, np.inf, + model, inliers = ransac(np.empty((10, 2)), LineModelND, 2, np.inf, is_model_valid=is_model_valid) assert_equal(model, None) assert_equal(inliers, None) @@ -248,8 +302,8 @@ def test_ransac_invalid_input(): def test_deprecated_params_attribute(): - model = LineModel() - model.params = (10, 1) + model = LineModelND() + model.params = ((0, 0), (1, 1)) x = np.arange(-10, 10) y = model.predict_y(x) with expected_warnings(['`_params`']): diff --git a/skimage/measure/tests/test_regionprops.py b/skimage/measure/tests/test_regionprops.py index 514b84d3..6865818e 100644 --- a/skimage/measure/tests/test_regionprops.py +++ b/skimage/measure/tests/test_regionprops.py @@ -128,14 +128,12 @@ def test_equiv_diameter(): def test_euler_number(): - with expected_warnings(['`background`|CObject type']): - en = regionprops(SAMPLE)[0].euler_number + en = regionprops(SAMPLE)[0].euler_number assert en == 0 SAMPLE_mod = SAMPLE.copy() SAMPLE_mod[7, -3] = 0 - with expected_warnings(['`background`|CObject type']): - en = regionprops(SAMPLE_mod)[0].euler_number + en = regionprops(SAMPLE_mod)[0].euler_number assert en == -1 @@ -374,9 +372,8 @@ def test_equals(): r2 = regions[0] r3 = regions[1] - with expected_warnings(['`background`|CObject type']): - assert_equal(r1 == r2, True, "Same regionprops are not equal") - assert_equal(r1 != r3, True, "Different regionprops are equal") + assert_equal(r1 == r2, True, "Same regionprops are not equal") + assert_equal(r1 != r3, True, "Different regionprops are equal") if __name__ == "__main__": diff --git a/skimage/morphology/__init__.py b/skimage/morphology/__init__.py index 9bf311c4..159a4ac2 100644 --- a/skimage/morphology/__init__.py +++ b/skimage/morphology/__init__.py @@ -8,7 +8,7 @@ from .watershed import watershed from ._skeletonize import skeletonize, medial_axis from .convex_hull import convex_hull_image, convex_hull_object from .greyreconstruct import reconstruction -from .misc import remove_small_objects +from .misc import remove_small_objects, remove_small_holes from ..measure._label import label from .._shared.utils import deprecated as _deprecated @@ -40,4 +40,5 @@ __all__ = ['binary_erosion', 'convex_hull_image', 'convex_hull_object', 'reconstruction', - 'remove_small_objects'] + 'remove_small_objects', + 'remove_small_holes'] diff --git a/skimage/morphology/convex_hull.py b/skimage/morphology/convex_hull.py index 5926e1e0..bb9ca734 100644 --- a/skimage/morphology/convex_hull.py +++ b/skimage/morphology/convex_hull.py @@ -1,11 +1,16 @@ __all__ = ['convex_hull_image', 'convex_hull_object'] import numpy as np -from ..measure import grid_points_in_poly +from ..measure._pnpoly import grid_points_in_poly from ._convex_hull import possible_hull from ..measure._label import label from ..util import unique_rows +try: + from scipy.spatial import Delaunay +except ImportError: + Delaunay = None + def convex_hull_image(image): """Compute the convex hull image of a binary image. @@ -15,12 +20,12 @@ def convex_hull_image(image): Parameters ---------- - image : ndarray + image : (M, N) array Binary input image. This array is cast to bool before processing. Returns ------- - hull : ndarray of bool + hull : (M, N) array of bool Binary image with pixels in convex hull set to True. References @@ -29,12 +34,13 @@ def convex_hull_image(image): """ - image = image.astype(bool) + if Delaunay is None: + raise ImportError("Could not import scipy.spatial.Delaunay, " + "only available in scipy >= 0.9.") # Here we do an optimisation by choosing only pixels that are # the starting or ending pixel of a row or column. This vastly - # limits the number of coordinates to examine for the virtual - # hull. + # limits the number of coordinates to examine for the virtual hull. coords = possible_hull(image.astype(np.uint8)) N = len(coords) @@ -48,12 +54,6 @@ def convex_hull_image(image): # scipy.spatial.Delaunay, so we remove them. coords = unique_rows(coords_corners) - try: - from scipy.spatial import Delaunay - except ImportError: - raise ImportError('Could not import scipy.spatial, only available in ' - 'scipy >= 0.9.') - # Subtract offset offset = coords.mean(axis=0) coords -= offset diff --git a/skimage/morphology/misc.py b/skimage/morphology/misc.py index e743a398..1a024f38 100644 --- a/skimage/morphology/misc.py +++ b/skimage/morphology/misc.py @@ -37,7 +37,12 @@ def default_selem(func): return func(image, selem=selem, *args, **kwargs) return func_out - + +def _check_dtype_supported(ar): + # Should use `issubdtype` for bool below, but there's a bug in numpy 1.7 + if not (ar.dtype == bool or np.issubdtype(ar.dtype, np.integer)): + raise TypeError("Only bool or integer image types are supported. " + "Got %s." % ar.dtype) def remove_small_objects(ar, min_size=64, connectivity=1, in_place=False): """Remove connected components smaller than the specified size. @@ -88,10 +93,8 @@ def remove_small_objects(ar, min_size=64, connectivity=1, in_place=False): >>> d is a True """ - # Should use `issubdtype` for bool below, but there's a bug in numpy 1.7 - if not (ar.dtype == bool or np.issubdtype(ar.dtype, np.integer)): - raise TypeError("Only bool or integer image types are supported. " - "Got %s." % ar.dtype) + # Raising type error if not int or bool + _check_dtype_supported(ar) if in_place: out = ar @@ -124,3 +127,89 @@ def remove_small_objects(ar, min_size=64, connectivity=1, in_place=False): out[too_small_mask] = 0 return out + +def remove_small_holes(ar, min_size=64, connectivity=1, in_place=False): + """Remove continguous holes smaller than the specified size. + + Parameters + ---------- + ar : ndarray (arbitrary shape, int or bool type) + The array containing the connected components of interest. + min_size : int, optional (default: 64) + The hole component size. + connectivity : int, {1, 2, ..., ar.ndim}, optional (default: 1) + The connectivity defining the neighborhood of a pixel. + in_place : bool, optional (default: False) + If `True`, remove the connected components in the input array itself. + Otherwise, make a copy. + + Raises + ------ + TypeError + If the input array is of an invalid type, such as float or string. + ValueError + If the input array contains negative values. + + Returns + ------- + out : ndarray, same shape and type as input `ar` + The input array with small holes within connected components removed. + + Examples + -------- + >>> from skimage import morphology + >>> a = np.array([[1, 1, 1, 1, 1, 0], + ... [1, 1, 1, 0, 1, 0], + ... [1, 0, 0, 1, 1, 0], + ... [1, 1, 1, 1, 1, 0]], bool) + >>> b = morphology.remove_small_holes(a, 2) + >>> b + array([[ True, True, True, True, True, False], + [ True, True, True, True, True, False], + [ True, False, False, True, True, False], + [ True, True, True, True, True, False]], dtype=bool) + >>> c = morphology.remove_small_holes(a, 2, connectivity=2) + >>> c + array([[ True, True, True, True, True, False], + [ True, True, True, False, True, False], + [ True, False, False, True, True, False], + [ True, True, True, True, True, False]], dtype=bool) + >>> d = morphology.remove_small_holes(a, 2, in_place=True) + >>> d is a + True + + Notes + ----- + + If the array type is int, it is assumed that it contains already-labeled + objects. The labels are not kept in the output image (this function always + outputs a bool image). It is suggested that labeling is completed after + using this function. + """ + _check_dtype_supported(ar) + + #Creates warning if image is an integer image + if ar.dtype != bool: + warnings.warn("Any labeled images will be returned as a boolean array. " + "Did you mean to use a boolean array?", UserWarning) + + if in_place: + out = ar + else: + out = ar.copy() + + #Creating the inverse of ar + if in_place: + out = np.logical_not(out,out) + else: + out = np.logical_not(out) + + #removing small objects from the inverse of ar + out = remove_small_objects(out, min_size, connectivity, in_place) + + if in_place: + out = np.logical_not(out,out) + else: + out = np.logical_not(out) + + return out diff --git a/skimage/morphology/selem.py b/skimage/morphology/selem.py index 53972cdd..8e516103 100644 --- a/skimage/morphology/selem.py +++ b/skimage/morphology/selem.py @@ -1,12 +1,8 @@ -""" -:author: Damian Eads, 2009 -:license: modified BSD -""" - import numpy as np from scipy import ndimage as ndi from .. import draw + def square(width, dtype=np.uint8): """Generates a flat, square-shaped structuring element. @@ -37,7 +33,7 @@ def rectangle(width, height, dtype=np.uint8): """Generates a flat, rectangular-shaped structuring element. Every pixel in the rectangle generated for a given width and given height - belongs to the neighboorhood. + belongs to the neighborhood. Parameters ---------- @@ -65,7 +61,7 @@ def diamond(radius, dtype=np.uint8): """Generates a flat, diamond-shaped structuring element. A pixel is part of the neighborhood (i.e. labeled 1) if - the city block/manhattan distance between it and the center of + the city block/Manhattan distance between it and the center of the neighborhood is no greater than radius. Parameters @@ -193,7 +189,7 @@ def octahedron(radius, dtype=np.uint8): This is the 3D equivalent of a diamond. A pixel is part of the neighborhood (i.e. labeled 1) if - the city block/manhattan distance between it and the center of + the city block/Manhattan distance between it and the center of the neighborhood is no greater than radius. Parameters diff --git a/skimage/morphology/tests/test_convex_hull.py b/skimage/morphology/tests/test_convex_hull.py index 67850cf2..c0cc954e 100644 --- a/skimage/morphology/tests/test_convex_hull.py +++ b/skimage/morphology/tests/test_convex_hull.py @@ -139,5 +139,6 @@ def test_object(): assert_raises(ValueError, convex_hull_object, image, 7) + if __name__ == "__main__": np.testing.run_module_suite() diff --git a/skimage/morphology/tests/test_misc.py b/skimage/morphology/tests/test_misc.py index 8789b22a..b7222cb5 100644 --- a/skimage/morphology/tests/test_misc.py +++ b/skimage/morphology/tests/test_misc.py @@ -1,7 +1,8 @@ import numpy as np from numpy.testing import (assert_array_equal, assert_equal, assert_raises, assert_warns) -from skimage.morphology import remove_small_objects +from skimage.morphology import remove_small_objects, remove_small_holes +from ..._shared._warnings import expected_warnings test_image = np.array([[0, 0, 0, 1, 0], [1, 1, 1, 0, 0], @@ -60,7 +61,8 @@ def test_single_label_warning(): image = np.array([[0, 0, 0, 1, 0], [1, 1, 1, 0, 0], [1, 1, 1, 0, 0]], int) - assert_warns(UserWarning, remove_small_objects, image, min_size=6) + with expected_warnings(['use a boolean array?']): + remove_small_objects(image, min_size=6) def test_float_input(): @@ -72,6 +74,102 @@ def test_negative_input(): negative_int = np.random.randint(-4, -1, size=(5, 5)) assert_raises(ValueError, remove_small_objects, negative_int) +test_holes_image = np.array([[0,0,0,0,0,0,1,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,0,0,1,1,0,0,0,0], + [0,1,1,1,0,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,0,0,0,0,0,0,1,1,1], + [0,0,0,0,0,0,0,1,0,1], + [0,0,0,0,0,0,0,1,1,1]], bool) +def test_one_connectivity_holes(): + expected = np.array([[0,0,0,0,0,0,1,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,0,0,0,0,0,0,1,1,1], + [0,0,0,0,0,0,0,1,1,1], + [0,0,0,0,0,0,0,1,1,1]], bool) + observed = remove_small_holes(test_holes_image, min_size=3) + assert_array_equal(observed, expected) + + +def test_two_connectivity_holes(): + expected = np.array([[0,0,0,0,0,0,1,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,0,0,1,1,0,0,0,0], + [0,1,1,1,0,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,0,0,0,0,0,0,1,1,1], + [0,0,0,0,0,0,0,1,1,1], + [0,0,0,0,0,0,0,1,1,1]], bool) + observed = remove_small_holes(test_holes_image, min_size=3, connectivity=2) + assert_array_equal(observed, expected) + +def test_in_place_holes(): + observed = remove_small_holes(test_holes_image, min_size=3, in_place=True) + assert_equal(observed is test_holes_image, True, + "remove_small_holes in_place argument failed.") + +def test_labeled_image_holes(): + labeled_holes_image = np.array([[0,0,0,0,0,0,1,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,0,0,1,1,0,0,0,0], + [0,1,1,1,0,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,0,0,0,0,0,0,2,2,2], + [0,0,0,0,0,0,0,2,0,2], + [0,0,0,0,0,0,0,2,2,2]], dtype=int) + expected = np.array([[0,0,0,0,0,0,1,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,0,0,0,0,0,0,1,1,1], + [0,0,0,0,0,0,0,1,1,1], + [0,0,0,0,0,0,0,1,1,1]], dtype=bool) + observed = remove_small_holes(labeled_holes_image, min_size=3) + assert_array_equal(observed, expected) + + +def test_uint_image_holes(): + labeled_holes_image = np.array([[0,0,0,0,0,0,1,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,0,0,1,1,0,0,0,0], + [0,1,1,1,0,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,0,0,0,0,0,0,2,2,2], + [0,0,0,0,0,0,0,2,0,2], + [0,0,0,0,0,0,0,2,2,2]], dtype=np.uint8) + expected = np.array([[0,0,0,0,0,0,1,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,0,0,0,0,0,0,1,1,1], + [0,0,0,0,0,0,0,1,1,1], + [0,0,0,0,0,0,0,1,1,1]], dtype=bool) + observed = remove_small_holes(labeled_holes_image, min_size=3) + assert_array_equal(observed, expected) + + +def test_label_warning_holes(): + labeled_holes_image = np.array([[0,0,0,0,0,0,1,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,1,0,0,1,1,0,0,0,0], + [0,1,1,1,0,1,0,0,0,0], + [0,1,1,1,1,1,0,0,0,0], + [0,0,0,0,0,0,0,2,2,2], + [0,0,0,0,0,0,0,2,0,2], + [0,0,0,0,0,0,0,2,2,2]], dtype=int) + with expected_warnings(['use a boolean array?']): + remove_small_holes(labeled_holes_image, min_size=3) + +def test_float_input_holes(): + float_test = np.random.rand(5, 5) + assert_raises(TypeError, remove_small_holes, float_test) + if __name__ == "__main__": np.testing.run_module_suite() diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index 1a7e1245..f36091bd 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -41,7 +41,7 @@ def _find_boundaries_subpixel(label_img): for index in np.ndindex(label_img_expanded.shape): if edges[index]: values = np.unique(windows[index].ravel()) - if len(values) > 2: # single value and max_label + if len(values) > 2: # single value and max_label boundaries[index] = True return boundaries @@ -51,9 +51,9 @@ def find_boundaries(label_img, connectivity=1, mode='thick', background=0): Parameters ---------- - label_img : array of int - An array in which different regions are labeled with different - integers. + label_img : array of int or bool + An array in which different regions are labeled with either different + integers or boolean values. connectivity: int in {1, ..., `label_img.ndim`}, optional A pixel is considered a boundary pixel if any of its neighbors has a different label. `connectivity` controls which pixels are @@ -144,7 +144,20 @@ def find_boundaries(label_img, connectivity=1, mode='thick', background=0): [0, 0, 0, 1, 0, 1, 0], [0, 0, 0, 1, 1, 1, 0], [0, 0, 0, 0, 0, 0, 0]], dtype=uint8) + >>> bool_image = np.array([[False, False, False, False, False], + ... [False, False, False, False, False], + ... [False, False, True, True, True], + ... [False, False, True, True, True], + ... [False, False, True, True, True]], dtype=np.bool) + >>> find_boundaries(bool_image) + array([[False, False, False, False, False], + [False, False, True, True, True], + [False, True, True, True, True], + [False, True, True, False, False], + [False, True, True, False, False]], dtype=bool) """ + if label_img.dtype == 'bool': + label_img = label_img.astype(np.uint8) ndim = label_img.ndim selem = ndi.generate_binary_structure(ndim, connectivity) if mode != 'subpixel': diff --git a/skimage/segmentation/tests/test_boundaries.py b/skimage/segmentation/tests/test_boundaries.py index e6e401ac..b8f7a4ae 100644 --- a/skimage/segmentation/tests/test_boundaries.py +++ b/skimage/segmentation/tests/test_boundaries.py @@ -25,6 +25,19 @@ def test_find_boundaries(): assert_array_equal(result, ref) +def test_find_boundaries_bool(): + image = np.zeros((5, 5), dtype=np.bool) + image[2:5, 2:5] = True + + ref = np.array([[False, False, False, False, False], + [False, False, True, True, True], + [False, True, True, True, True], + [False, True, True, False, False], + [False, True, True, False, False]], dtype=np.bool) + result = find_boundaries(image) + assert_array_equal(result, ref) + + def test_mark_boundaries(): image = np.zeros((10, 10)) label_image = np.zeros((10, 10), dtype=np.uint8) @@ -61,6 +74,27 @@ def test_mark_boundaries(): assert_array_equal(result, ref) +def test_mark_boundaries_bool(): + image = np.zeros((10, 10), dtype=np.bool) + label_image = np.zeros((10, 10), dtype=np.uint8) + label_image[2:7, 2:7] = 1 + + ref = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 0, 0], + [0, 1, 1, 0, 0, 0, 1, 1, 0, 0], + [0, 1, 1, 0, 0, 0, 1, 1, 0, 0], + [0, 1, 1, 0, 0, 0, 1, 1, 0, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) + + marked = mark_boundaries(image, label_image, color=white, mode='thick') + result = np.mean(marked, axis=-1) + assert_array_equal(result, ref) + + def test_mark_boundaries_subpixel(): labels = np.array([[0, 0, 0, 0], [0, 0, 5, 0], diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index 9b40b3a9..3df5863d 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -1,6 +1,6 @@ -from ._hough_transform import (hough_ellipse, hough_line, - probabilistic_hough_line) -from .hough_transform import hough_circle, hough_line_peaks +from .hough_transform import (hough_line, hough_line_peaks, + probabilistic_hough_line, hough_circle, + hough_ellipse) from .radon_transform import radon, iradon, iradon_sart from .finite_radon_transform import frt2, ifrt2 from .integral import integral_image, integrate diff --git a/skimage/transform/_hough_transform.pyx b/skimage/transform/_hough_transform.pyx index 6ce69d8e..a189d7b9 100644 --- a/skimage/transform/_hough_transform.pyx +++ b/skimage/transform/_hough_transform.pyx @@ -12,9 +12,6 @@ from libc.stdlib cimport rand from ..draw import circle_perimeter -cdef double PI_2 = 1.5707963267948966 -cdef double NEG_PI_2 = -PI_2 - from .._shared.interpolation cimport round @@ -278,16 +275,10 @@ def hough_line(cnp.ndarray img, .. plot:: hough_tf.py """ - if img.ndim != 2: - raise ValueError('The input image must be 2D.') - # Compute the array of angles and their sine and cosine cdef cnp.ndarray[ndim=1, dtype=cnp.double_t] ctheta cdef cnp.ndarray[ndim=1, dtype=cnp.double_t] stheta - if theta is None: - theta = np.linspace(NEG_PI_2, PI_2, 180) - ctheta = np.cos(theta) stheta = np.sin(theta) @@ -354,12 +345,6 @@ def probabilistic_hough_line(cnp.ndarray img, int threshold=10, Hough transform for line detection", in IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 1999. """ - if img.ndim != 2: - raise ValueError('The input image must be 2D.') - - if theta is None: - theta = PI_2 - np.arange(180) / 180.0 * 2 * PI_2 - cdef Py_ssize_t height = img.shape[0] cdef Py_ssize_t width = img.shape[1] diff --git a/skimage/transform/hough_transform.py b/skimage/transform/hough_transform.py index 92f2c79a..f9a8d2e5 100644 --- a/skimage/transform/hough_transform.py +++ b/skimage/transform/hough_transform.py @@ -1,7 +1,67 @@ import numpy as np from scipy import ndimage as ndi -from .. import measure, morphology -from ._hough_transform import _hough_circle +from .. import measure +from ._hough_transform import (_hough_circle, + hough_ellipse as _hough_ellipse, + hough_line as _hough_line, + probabilistic_hough_line as _prob_hough_line) + + +# Wrapper for Cython allows function signature introspection +def hough_line(img, theta=None): + """Perform a straight line Hough transform. + + Parameters + ---------- + img : (M, N) ndarray + Input image with nonzero values representing edges. + theta : 1D ndarray of double + Angles at which to compute the transform, in radians. + Defaults to -pi/2 .. pi/2 + + Returns + ------- + H : 2-D ndarray of uint64 + Hough transform accumulator. + theta : ndarray + Angles at which the transform was computed, in radians. + distances : ndarray + Distance values. + + Notes + ----- + The origin is the top left corner of the original image. + X and Y axis are horizontal and vertical edges respectively. + The distance is the minimal algebraic distance from the origin + to the detected line. + + Examples + -------- + Generate a test image: + + >>> img = np.zeros((100, 150), dtype=bool) + >>> img[30, :] = 1 + >>> img[:, 65] = 1 + >>> img[35:45, 35:50] = 1 + >>> for i in range(90): + ... img[i, i] = 1 + >>> img += np.random.random(img.shape) > 0.95 + + Apply the Hough transform: + + >>> out, angles, d = hough_line(img) + + .. plot:: hough_tf.py + + """ + if img.ndim != 2: + raise ValueError('The input image `img` must be 2D.') + + if theta is None: + # These values are approximations of pi/2 + theta = np.linspace(-np.pi / 2, np.pi / 2, 180) + + return _hough_line(img, theta=theta) def hough_line_peaks(hspace, angles, dists, min_distance=9, min_angle=10, @@ -73,7 +133,10 @@ def hough_line_peaks(hspace, angles, dists, min_distance=9, min_angle=10, label_hspace = measure.label(hspace_t) props = measure.regionprops(label_hspace, hspace_max) - props = sorted(props, key= lambda x: x.max_intensity)[::-1] + + # Sort the list of peaks by intensity, not left-right, so larger peaks + # in Hough space cannot be arbitrarily suppressed by smaller neighbors + props = sorted(props, key=lambda x: x.max_intensity)[::-1] coords = np.array([np.round(p.centroid) for p in props], dtype=int) hspace_peaks = [] @@ -126,6 +189,48 @@ def hough_line_peaks(hspace, angles, dists, min_distance=9, min_angle=10, return hspace_peaks, angle_peaks, dist_peaks +# Wrapper for Cython allows function signature introspection +def probabilistic_hough_line(img, threshold=10, line_length=50, line_gap=10, + theta=None): + """Return lines from a progressive probabilistic line Hough transform. + + Parameters + ---------- + img : (M, N) ndarray + Input image with nonzero values representing edges. + threshold : int, optional (default 10) + Threshold + line_length : int, optional (default 50) + Minimum accepted length of detected lines. + Increase the parameter to extract longer lines. + line_gap : int, optional, (default 10) + Maximum gap between pixels to still form a line. + Increase the parameter to merge broken lines more aggresively. + theta : 1D ndarray, dtype=double, optional, default (-pi/2 .. pi/2) + Angles at which to compute the transform, in radians. + + Returns + ------- + lines : list + List of lines identified, lines in format ((x0, y0), (x1, y0)), + indicating line start and end. + + References + ---------- + .. [1] C. Galamhos, J. Matas and J. Kittler, "Progressive probabilistic + Hough transform for line detection", in IEEE Computer Society + Conference on Computer Vision and Pattern Recognition, 1999. + """ + if img.ndim != 2: + raise ValueError('The input image `img` must be 2D.') + + if theta is None: + theta = np.pi / 2 - np.arange(180) / 180.0 * np.pi + + return _prob_hough_line(img, threshold=threshold, line_length=line_length, + line_gap=line_gap, theta=theta) + + def hough_circle(image, radius, normalize=True, full_output=False): """Perform a circular Hough transform. @@ -169,3 +274,56 @@ def hough_circle(image, radius, normalize=True, full_output=False): radius = np.atleast_1d(np.asarray(radius)) return _hough_circle(image, radius.astype(np.intp), normalize=normalize, full_output=full_output) + + +# Wrapper for Cython allows function signature introspection +def hough_ellipse(img, threshold=4, accuracy=1, min_size=4, max_size=None): + """Perform an elliptical Hough transform. + + Parameters + ---------- + img : (M, N) ndarray + Input image with nonzero values representing edges. + threshold: int, optional (default 4) + Accumulator threshold value. + accuracy : double, optional (default 1) + Bin size on the minor axis used in the accumulator. + min_size : int, optional (default 4) + Minimal major axis length. + max_size : int, optional + Maximal minor axis length. (default None) + If None, the value is set to the half of the smaller + image dimension. + + Returns + ------- + result : ndarray with fields [(accumulator, y0, x0, a, b, orientation)] + Where ``(yc, xc)`` is the center, ``(a, b)`` the major and minor + axes, respectively. The `orientation` value follows + `skimage.draw.ellipse_perimeter` convention. + + Examples + -------- + >>> from skimage.transform import hough_ellipse + >>> from skimage.draw import ellipse_perimeter + >>> img = np.zeros((25, 25), dtype=np.uint8) + >>> rr, cc = ellipse_perimeter(10, 10, 6, 8) + >>> img[cc, rr] = 1 + >>> result = hough_ellipse(img, threshold=8) + >>> result.tolist() + [(10, 10.0, 10.0, 8.0, 6.0, 0.0)] + + Notes + ----- + The accuracy must be chosen to produce a peak in the accumulator + distribution. In other words, a flat accumulator distribution with low + values may be caused by a too low bin size. + + References + ---------- + .. [1] Xie, Yonghong, and Qiang Ji. "A new efficient ellipse detection + method." Pattern Recognition, 2002. Proceedings. 16th International + Conference on. Vol. 2. IEEE, 2002 + """ + return _hough_ellipse(img, threshold=threshold, accuracy=accuracy, + min_size=min_size, max_size=max_size) diff --git a/skimage/transform/tests/test_hough_transform.py b/skimage/transform/tests/test_hough_transform.py index fce7013c..0babc769 100644 --- a/skimage/transform/tests/test_hough_transform.py +++ b/skimage/transform/tests/test_hough_transform.py @@ -1,5 +1,5 @@ import numpy as np -from numpy.testing import assert_almost_equal, assert_equal +from numpy.testing import assert_almost_equal, assert_equal, assert_raises import skimage.transform as tf from skimage.draw import line, circle_perimeter, ellipse_perimeter @@ -7,15 +7,6 @@ from skimage._shared._warnings import expected_warnings from skimage._shared.testing import test_parallel -def append_desc(func, description): - """Append the test function ``func`` and append - ``description`` to its name. - """ - func.description = func.__module__ + '.' + func.__name__ + description - - return func - - @test_parallel() def test_hough_line(): # Generate a test image @@ -42,12 +33,21 @@ def test_hough_line_angles(): assert_equal(len(angles), 10) +def test_hough_line_bad_input(): + img = np.zeros(100) + img[10] = 1 + + # Expected error, img must be 2D + assert_raises(ValueError, tf.hough_line, img, np.linspace(0, 360, 10)) + + def test_probabilistic_hough(): # Generate a test image img = np.zeros((100, 100), dtype=int) for i in range(25, 75): img[100 - i, i] = 100 img[i, i] = 100 + # decrease default theta sampling because similar orientations may confuse # as mentioned in article of Galambos et al theta = np.linspace(0, np.pi, 45) @@ -59,9 +59,21 @@ def test_probabilistic_hough(): line = list(line) line.sort(key=lambda x: x[0]) sorted_lines.append(line) + assert([(25, 75), (74, 26)] in sorted_lines) assert([(25, 25), (74, 74)] in sorted_lines) + # Execute with default theta + tf.probabilistic_hough_line(img, line_length=10, line_gap=3) + + +def test_probabilistic_hough_bad_input(): + img = np.zeros(100) + img[10] = 1 + + # Expected error, img must be 2D + assert_raises(ValueError, tf.probabilistic_hough_line, img) + def test_hough_line_peaks(): img = np.zeros((100, 150), dtype=int) @@ -78,6 +90,22 @@ def test_hough_line_peaks(): assert_almost_equal(theta[0], 1.41, 1) +def test_hough_line_peaks_ordered(): + # Regression test per PR #1421 + testim = np.zeros((256, 64), dtype=np.bool) + + testim[50:100, 20] = True + testim[85:200, 25] = True + testim[15:35, 50] = True + testim[1:-1, 58] = True + + hough_space, angles, dists = tf.hough_line(testim) + + with expected_warnings(['`background`']): + hspace, _, _ = tf.hough_line_peaks(hough_space, angles, dists) + assert hspace[0] > hspace[1] + + def test_hough_line_peaks_dist(): img = np.zeros((100, 100), dtype=np.bool_) img[:, 30] = True