Merge branch 'master' of git://github.com/scikit-image/scikit-image

This commit is contained in:
François Boulogne
2015-12-06 16:50:48 -05:00
80 changed files with 1273 additions and 284 deletions
+3
View File
@@ -215,3 +215,6 @@
- Jim Fienup, Alexander Iacchetta
In-depth review of sub-pixel shift registration
- Damian Eads
Structuring elements in morphology module.
-3
View File
@@ -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``.
+4
View File
@@ -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
+3 -5
View File
@@ -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
+3 -2
View File
@@ -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
@@ -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)
+3 -1
View File
@@ -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
+72 -38
View File
@@ -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
+3 -3
View File
@@ -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
+5 -1
View File
@@ -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:
+1 -1
View File
@@ -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')
@@ -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)
+1 -1
View File
@@ -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()
+2 -2
View File
@@ -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')
+49 -8
View File
@@ -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()
+8 -1
View File
@@ -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')
+3 -1
View File
@@ -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()
+15 -13
View File
@@ -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()
"""
+4 -2
View File
@@ -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')
+1 -1
View File
@@ -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')
+28 -25
View File
@@ -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()
+8 -1
View File
@@ -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')
+1 -1
View File
@@ -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),
+8 -1
View File
@@ -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')
+3 -3
View File
@@ -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()
+1 -1
View File
@@ -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')
+1 -1
View File
@@ -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')
+6 -1
View File
@@ -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')
+1 -1
View File
@@ -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')
+1 -1
View File
@@ -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)
+2 -2
View File
@@ -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)
@@ -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,
+9 -12
View File
@@ -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()
+40
View File
@@ -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()
+5 -2
View File
@@ -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
+8 -2
View File
@@ -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()
+1 -1
View File
@@ -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()
+9 -8
View File
@@ -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
"""
+1 -1
View File
@@ -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)
+1 -1
View File
@@ -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')
+1 -1
View File
@@ -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())
+1 -1
View File
@@ -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')
+4 -1
View File
@@ -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()
@@ -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()
+10 -5
View File
@@ -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()
+1 -1
View File
@@ -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')
+3 -2
View File
@@ -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' +
-2
View File
@@ -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.")
+2
View File
@@ -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)
Binary file not shown.
+6
View File
@@ -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
+1 -1
View File
@@ -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
+1 -1
View File
@@ -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
+2 -1
View File
@@ -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',
+2 -2
View File
@@ -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
+34
View File
@@ -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)
+4 -3
View File
@@ -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',
+45 -2
View File
@@ -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):
+31 -1
View File
@@ -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] = <dtype_t_out>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] = <dtype_t_out>round(exp(mean / pop)-1)
else:
out[0] = <dtype_t_out>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,
+33 -2
View File
@@ -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)
+28 -1
View File
@@ -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
+2 -1
View File
@@ -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',
+1 -1
View File
@@ -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__
+5 -1
View File
@@ -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, :]
+2 -2
View File
@@ -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)
+170
View File
@@ -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.
+80 -26
View File
@@ -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`']):
+4 -7
View File
@@ -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__":
+3 -2
View File
@@ -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']
+12 -12
View File
@@ -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
+94 -5
View File
@@ -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
+4 -8
View File
@@ -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
@@ -139,5 +139,6 @@ def test_object():
assert_raises(ValueError, convex_hull_object, image, 7)
if __name__ == "__main__":
np.testing.run_module_suite()
+100 -2
View File
@@ -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()
+17 -4
View File
@@ -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':
@@ -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],
+3 -3
View File
@@ -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
-15
View File
@@ -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]
+161 -3
View File
@@ -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)
+38 -10
View File
@@ -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