diff --git a/.mailmap b/.mailmap new file mode 100644 index 00000000..903ad2df --- /dev/null +++ b/.mailmap @@ -0,0 +1,18 @@ +K.-Michael Aye +Nelson Brown +Luis Pedro Coelho +Marianne Corvellec +Riaan van den Dool +Emmanuelle Gouillart +Thouis (Ray) Jones +Gregory R. Lee +Andreas Mueller +Juan Nunez-Iglesias +Nicolas Pinto +Johannes Schönberger +Tim Sheerman-Chase +Matthew Trentacoste +James Turner +Stefan van der Walt +John Wiggins +Tony S Yu diff --git a/CONTRIBUTING.txt b/CONTRIBUTING.txt index aefde504..29dba55d 100644 --- a/CONTRIBUTING.txt +++ b/CONTRIBUTING.txt @@ -225,6 +225,73 @@ Every time Travis is triggered, it also calls on `Coveralls `_ to inspect the current test overage. +Building docs +------------- + +To build docs, run ``make`` from the ``docs`` directory. ``make help`` lists +all targets. + +Requirements +~~~~~~~~~~~~ + +Sphinx (>= 1.3) and Latex is needed to build doc. + +**Sphinx:** + +.. code:: sh + + pip install sphinx + +**Latex Ubuntu:** + +.. code:: sh + + sudo apt-get install -qq texlive texlive-latex-extra dvipng + +**Latex Mac:** + +Install the full `MacTex `__ installation or +install the smaller +`BasicTex `__ and add *ucs* +and *dvipng* packages: + +.. code:: sh + + sudo tlmgr install ucs dvipng + +Fixing Warnings +~~~~~~~~~~~~~~~ + +- "citation not found: R###" There is probably an underscore after a + reference 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 the docstrings + +- Make sure to use pre-sphinxification paths to images (not the + \_images directory) + +Auto-generating dev docs +~~~~~~~~~~~~~~~~~~~~~~~~ + +This set of instructions was used to create +scikit-image/tools/deploy-docs.sh + +- Go to Github account settings -> personal access tokens +- Create a new token with access rights ``public_repo`` and + ``user:email only`` +- Install the travis command line tool: ``gem install travis``. On OSX, + you can get gem via ``brew install ruby``. +- Take then token generated by Github and run + ``travis encrypt GH_TOKEN=`` from inside a scikit-image repo +- Paste the output into the secure: field of ``.travis.yml``. +- The decrypted GH\_TOKEN env var will be available for travis scripts + +https://help.github.com/articles/creating-an-access-token-for-command-line-use/ +http://docs.travis-ci.com/user/encryption-keys/ + + Bugs ---- diff --git a/DEPENDS.txt b/DEPENDS.txt index 2130e21d..2c2aeb7d 100644 --- a/DEPENDS.txt +++ b/DEPENDS.txt @@ -56,3 +56,10 @@ Testing requirements A Python Unit Testing Framework * `Coverage.py `__ A tool that generates a unit test code coverage report + + +Documentation requirements +-------------------------- + +`sphinx >= 1.3 `_ is required to build the +documentation. diff --git a/TODO.txt b/TODO.txt index 02bd9ead..bed2793a 100644 --- a/TODO.txt +++ b/TODO.txt @@ -12,10 +12,13 @@ Version 0.14 add an alias LineModel = LineModelND. While the deprecated LineModel has for parameters `(dist, theta)`, LineModelND has the more general parameters `(origin, direction)`. +* Remove deprecated old syntax support for ``skimage.transform.integrate``. Version 0.13 ------------ +* Require Python 2.7+, remove warning in `__init__.py` and 2.6 hack in + `doc/release/contribs.py`. * Remove deprecated `None` defaults for `skimage.exposure.rescale_intensity` * Remove deprecated `skimage.filters.canny` import in `filters/__init__.py` file (canny is now in `skimage.feature.canny`). @@ -36,8 +39,6 @@ Version 0.12 ------------ * Change `label` to mark background as 0, not -1, which is consistent with SciPy's labelling. -* Remove `skimage.morphology.label` from `skimage.morphology.__init__`--it now - lives in `skimage.measure.label`. * Remove deprecated `reverse_map` parameter of `skimage.transform.warp` * Change deprecated `enforce_connectivity=False` on skimage.segmentation.slic and set it to True as default diff --git a/doc/README.md b/doc/README.md deleted file mode 100644 index d9fec89a..00000000 --- a/doc/README.md +++ /dev/null @@ -1,53 +0,0 @@ -# Building docs # -To build docs, run `make` in this directory. `make help` lists all targets. - -## Requirements ## -Sphinx and Latex is needed to build doc. - -**Sphinx:** -```sh -pip install sphinx -``` - -**Latex Ubuntu:** -```sh -sudo apt-get install -qq texlive texlive-latex-extra dvipng -``` - -**Latex Mac:** - -Install the full [MacTex](http://www.tug.org/mactex/) installation or install the smaller [BasicTex](http://www.tug.org/mactex/morepackages.html) and add *ucs* and *dvipng* packages: -```sh -sudo tlmgr install ucs dvipng -``` - - -## Fixing Warnings ## - -- "citation not found: R###" - There is probably an underscore after a reference - 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 - the docstrings - -- Make sure to use pre-sphinxification paths to images - (not the _images directory) - - -## Auto-generating dev docs ## - -This set of instructions was used to create scikit-image/tools/deploy-docs.sh - -- Go to Github account settings -> personal access tokens -- Create a new token with access rights `public_repo` and `user:email only` -- Install the travis command line tool: `gem install travis`. On OSX, you can get gem via `brew install ruby`. -- Take then token generated by Github and run `travis encrypt GH_TOKEN=` from inside a scikit-image repo -- Paste the output into the secure: field of `.travis.yml`. -- The decrypted GH_TOKEN env var will be available for travis scripts - -https://help.github.com/articles/creating-an-access-token-for-command-line-use/ -http://docs.travis-ci.com/user/encryption-keys/ diff --git a/doc/examples/color_exposure/README.txt b/doc/examples/color_exposure/README.txt new file mode 100644 index 00000000..cfaa1da3 --- /dev/null +++ b/doc/examples/color_exposure/README.txt @@ -0,0 +1,2 @@ +Manipulating exposure and color channels +---------------------------------------- diff --git a/doc/examples/plot_adapt_rgb.py b/doc/examples/color_exposure/plot_adapt_rgb.py similarity index 100% rename from doc/examples/plot_adapt_rgb.py rename to doc/examples/color_exposure/plot_adapt_rgb.py diff --git a/doc/examples/plot_equalize.py b/doc/examples/color_exposure/plot_equalize.py similarity index 98% rename from doc/examples/plot_equalize.py rename to doc/examples/color_exposure/plot_equalize.py index 2289e270..079b8b6d 100644 --- a/doc/examples/plot_equalize.py +++ b/doc/examples/color_exposure/plot_equalize.py @@ -99,5 +99,5 @@ ax_cdf.set_ylabel('Fraction of total intensity') ax_cdf.set_yticks(np.linspace(0, 1, 5)) # prevent overlap of y-axis labels -fig.subplots_adjust(wspace=0.4) +fig.tight_layout() plt.show() diff --git a/doc/examples/plot_ihc_color_separation.py b/doc/examples/color_exposure/plot_ihc_color_separation.py similarity index 95% rename from doc/examples/plot_ihc_color_separation.py rename to doc/examples/color_exposure/plot_ihc_color_separation.py index 6191a800..2db2b552 100644 --- a/doc/examples/plot_ihc_color_separation.py +++ b/doc/examples/color_exposure/plot_ihc_color_separation.py @@ -26,7 +26,8 @@ from skimage.color import rgb2hed ihc_rgb = data.immunohistochemistry() ihc_hed = rgb2hed(ihc_rgb) -fig, axes = plt.subplots(2, 2, figsize=(7, 6), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +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) @@ -44,7 +45,7 @@ ax3.set_title("DAB") for ax in axes.ravel(): ax.axis('off') -fig.subplots_adjust(hspace=0.3) +fig.tight_layout() """ diff --git a/doc/examples/plot_local_equalize.py b/doc/examples/color_exposure/plot_local_equalize.py similarity index 84% rename from doc/examples/plot_local_equalize.py rename to doc/examples/color_exposure/plot_local_equalize.py index 194bfe6a..f1d1fb21 100644 --- a/doc/examples/plot_local_equalize.py +++ b/doc/examples/color_exposure/plot_local_equalize.py @@ -74,12 +74,14 @@ img_eq = rank.equalize(img, selem=selem) # Display results 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) +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') @@ -94,5 +96,5 @@ ax_cdf.set_ylabel('Fraction of total intensity') # prevent overlap of y-axis labels -fig.subplots_adjust(wspace=0.4) +fig.tight_layout() plt.show() diff --git a/doc/examples/plot_log_gamma.py b/doc/examples/color_exposure/plot_log_gamma.py similarity index 92% rename from doc/examples/plot_log_gamma.py rename to doc/examples/color_exposure/plot_log_gamma.py index 70d5881c..2efa5193 100644 --- a/doc/examples/plot_log_gamma.py +++ b/doc/examples/color_exposure/plot_log_gamma.py @@ -55,10 +55,12 @@ logarithmic_corrected = exposure.adjust_log(img, 1) # Display results fig = plt.figure(figsize=(8, 5)) -axes = np.zeros((2,3), dtype=np.object) +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[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) @@ -80,5 +82,5 @@ ax_cdf.set_ylabel('Fraction of total intensity') ax_cdf.set_yticks(np.linspace(0, 1, 5)) # prevent overlap of y-axis labels -fig.subplots_adjust(wspace=0.4) +fig.tight_layout() plt.show() diff --git a/doc/examples/plot_regional_maxima.py b/doc/examples/color_exposure/plot_regional_maxima.py similarity index 100% rename from doc/examples/plot_regional_maxima.py rename to doc/examples/color_exposure/plot_regional_maxima.py diff --git a/doc/examples/plot_tinting_grayscale_images.py b/doc/examples/color_exposure/plot_tinting_grayscale_images.py similarity index 100% rename from doc/examples/plot_tinting_grayscale_images.py rename to doc/examples/color_exposure/plot_tinting_grayscale_images.py diff --git a/doc/examples/edges/README.txt b/doc/examples/edges/README.txt new file mode 100644 index 00000000..ff6d7e2c --- /dev/null +++ b/doc/examples/edges/README.txt @@ -0,0 +1,2 @@ +Edges and lines +--------------- diff --git a/doc/examples/edges/plot_active_contours.py b/doc/examples/edges/plot_active_contours.py new file mode 100644 index 00000000..f6f2c1f6 --- /dev/null +++ b/doc/examples/edges/plot_active_contours.py @@ -0,0 +1,99 @@ +""" +==================== +Active Contour Model +==================== + +The active contour model is a method to fit open or closed splines to lines or +edges in an image. It works by minimising an energy that is in part defined by +the image and part by the spline's shape: length and smoothness. The +minimization is done implicitly in the shape energy and explicitly in the +image energy. + +In the following two examples the active contour model is used (1) to segment +the face of a person from the rest of an image by fitting a closed curve +to the edges of the face and (2) to find the darkest curve between two fixed +points while obeying smoothness considerations. Typically it is a good idea to +smooth images a bit before analyzing, as done in the following examples. + +.. [1] *Snakes: Active contour models*. Kass, M.; Witkin, A.; Terzopoulos, D. + International Journal of Computer Vision 1 (4): 321 (1988). + +We initialize a circle around the astronaut's face and use the default boundary +condition ``bc='periodic'`` to fit a closed curve. The default parameters +``w_line=0, w_edge=1`` will make the curve search towards edges, such as the +boundaries of the face. +""" + +import numpy as np +import matplotlib.pyplot as plt +from skimage.color import rgb2gray +from skimage import data +from skimage.filters import gaussian_filter +from skimage.segmentation import active_contour + +# Test scipy version, since active contour is only possible +# with recent scipy version +import scipy +scipy_version = list(map(int, scipy.__version__.split('.'))) +new_scipy = scipy_version[0] > 0 or \ + (scipy_version[0] == 0 and scipy_version[1] >= 14) + +img = data.astronaut() +img = rgb2gray(img) + +s = np.linspace(0, 2*np.pi, 400) +x = 220 + 100*np.cos(s) +y = 100 + 100*np.sin(s) +init = np.array([x, y]).T + +if not new_scipy: + print('You are using an old version of scipy. ' + 'Active contours is implemented for scipy versions ' + '0.14.0 and above.') + +if new_scipy: + snake = active_contour(gaussian_filter(img, 3), + init, alpha=0.015, beta=10, gamma=0.001) + + fig = plt.figure(figsize=(7, 7)) + ax = fig.add_subplot(111) + plt.gray() + ax.imshow(img) + ax.plot(init[:, 0], init[:, 1], '--r', lw=3) + ax.plot(snake[:, 0], snake[:, 1], '-b', lw=3) + ax.set_xticks([]), ax.set_yticks([]) + ax.axis([0, img.shape[1], img.shape[0], 0]) + +""" +.. image:: PLOT2RST.current_figure + +Here we initialize a straight line between two points, `(5, 136)` and +`(424, 50)`, and require that the spline has its end points there by giving +the boundary condition `bc='fixed'`. We furthermore make the algorithm search +for dark lines by giving a negative `w_line` value. +""" + +img = data.text() + +x = np.linspace(5, 424, 100) +y = np.linspace(136, 50, 100) +init = np.array([x, y]).T + +if new_scipy: + snake = active_contour(gaussian_filter(img, 1), init, bc='fixed', + alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1) + + fig = plt.figure(figsize=(9, 5)) + ax = fig.add_subplot(111) + plt.gray() + ax.imshow(img) + ax.plot(init[:, 0], init[:, 1], '--r', lw=3) + ax.plot(snake[:, 0], snake[:, 1], '-b', lw=3) + ax.set_xticks([]), ax.set_yticks([]) + ax.axis([0, img.shape[1], img.shape[0], 0]) + +plt.show() + +""" +.. image:: PLOT2RST.current_figure +""" diff --git a/doc/examples/plot_canny.py b/doc/examples/edges/plot_canny.py similarity index 91% rename from doc/examples/plot_canny.py rename to doc/examples/edges/plot_canny.py index 06d2a05d..123ea874 100644 --- a/doc/examples/plot_canny.py +++ b/doc/examples/edges/plot_canny.py @@ -35,7 +35,8 @@ 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), sharex=True, sharey=True) +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') @@ -49,7 +50,6 @@ ax3.imshow(edges2, cmap=plt.cm.gray) ax3.axis('off') ax3.set_title('Canny filter, $\sigma=3$', fontsize=20) -fig.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, - bottom=0.02, left=0.02, right=0.98) +fig.tight_layout() plt.show() diff --git a/doc/examples/plot_circular_elliptical_hough_transform.py b/doc/examples/edges/plot_circular_elliptical_hough_transform.py similarity index 100% rename from doc/examples/plot_circular_elliptical_hough_transform.py rename to doc/examples/edges/plot_circular_elliptical_hough_transform.py diff --git a/doc/examples/plot_contours.py b/doc/examples/edges/plot_contours.py similarity index 100% rename from doc/examples/plot_contours.py rename to doc/examples/edges/plot_contours.py diff --git a/doc/examples/plot_convex_hull.py b/doc/examples/edges/plot_convex_hull.py similarity index 100% rename from doc/examples/plot_convex_hull.py rename to doc/examples/edges/plot_convex_hull.py diff --git a/doc/examples/plot_edge_filter.py b/doc/examples/edges/plot_edge_filter.py similarity index 100% rename from doc/examples/plot_edge_filter.py rename to doc/examples/edges/plot_edge_filter.py diff --git a/doc/examples/plot_line_hough_transform.py b/doc/examples/edges/plot_line_hough_transform.py similarity index 72% rename from doc/examples/plot_line_hough_transform.py rename to doc/examples/edges/plot_line_hough_transform.py index 15464768..c4d7cf95 100644 --- a/doc/examples/plot_line_hough_transform.py +++ b/doc/examples/edges/plot_line_hough_transform.py @@ -1,4 +1,4 @@ -r""" +""" ============================= Straight line Hough transform ============================= @@ -6,7 +6,7 @@ Straight line Hough transform The Hough transform in its simplest form is a `method to detect straight lines `__. -In the following example, we construct an image with a line intersection. We +In the following example, we construct an image with a line intersection. We then use the Hough transform to explore a parameter space for straight lines that may run through the image. @@ -53,9 +53,9 @@ References .. [2] Duda, R. O. and P. E. Hart, "Use of the Hough Transformation to Detect Lines and Curves in Pictures," Comm. ACM, Vol. 15, pp. 11-15 (January, 1972) - """ +from matplotlib import cm from skimage.transform import (hough_line, hough_line_peaks, probabilistic_hough_line) from skimage.feature import canny @@ -64,70 +64,71 @@ from skimage import data import numpy as np import matplotlib.pyplot as plt -# Construct test image - +# Constructing test image. image = np.zeros((100, 100)) - - -# Classic straight-line Hough transform - idx = np.arange(25, 75) image[idx[::-1], idx] = 255 image[idx, idx] = 255 +# Classic straight-line Hough transform. h, theta, d = hough_line(image) -fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8,4)) +# Generating figure 1. +fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(12, 6)) +plt.tight_layout() -ax1.imshow(image, cmap=plt.cm.gray) -ax1.set_title('Input image') -ax1.set_axis_off() +ax0.imshow(image, cmap=cm.gray) +ax0.set_title('Input image') +ax0.set_axis_off() -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) -ax2.set_title('Hough transform') -ax2.set_xlabel('Angles (degrees)') -ax2.set_ylabel('Distance (pixels)') -ax2.axis('image') +ax1.imshow(np.log(1 + h), extent=[np.rad2deg(theta[-1]), np.rad2deg(theta[0]), + d[-1], d[0]], cmap=cm.gray, aspect=1/1.5) +ax1.set_title('Hough transform') +ax1.set_xlabel('Angles (degrees)') +ax1.set_ylabel('Distance (pixels)') +ax1.axis('image') -ax3.imshow(image, cmap=plt.cm.gray) -rows, cols = image.shape +ax2.imshow(image, cmap=cm.gray) +row1, col1 = 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) - 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 + y1 = (dist - col1 * np.cos(angle)) / np.sin(angle) + ax2.plot((0, col1), (y0, y1), '-r') +ax2.axis((0, col1, row1, 0)) +ax2.set_title('Detected lines') +ax2.set_axis_off() +# Line finding using the Probabilistic Hough Transform. image = data.camera() edges = canny(image, 2, 1, 25) lines = probabilistic_hough_line(edges, threshold=10, line_length=5, line_gap=3) -fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8,4), sharex=True, sharey=True) +# Generating figure 2. +fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(16, 6), sharex=True, + sharey=True) +plt.tight_layout() -ax1.imshow(image, cmap=plt.cm.gray) -ax1.set_title('Input image') +ax0.imshow(image, cmap=cm.gray) +ax0.set_title('Input image') +ax0.set_axis_off() +ax0.set_adjustable('box-forced') + +ax1.imshow(edges, cmap=cm.gray) +ax1.set_title('Canny edges') ax1.set_axis_off() ax1.set_adjustable('box-forced') -ax2.imshow(edges, cmap=plt.cm.gray) -ax2.set_title('Canny edges') +ax2.imshow(edges * 0) +for line in lines: + p0, p1 = line + ax2.plot((p0[0], p1[0]), (p0[1], p1[1])) + +row2, col2 = image.shape +ax2.axis((0, col2, row2, 0)) + +ax2.set_title('Probabilistic Hough') ax2.set_axis_off() ax2.set_adjustable('box-forced') -ax3.imshow(edges * 0) - -for line in lines: - p0, p1 = line - ax3.plot((p0[0], p1[0]), (p0[1], p1[1])) - -ax3.set_title('Probabilistic Hough') -ax3.set_axis_off() -ax3.set_adjustable('box-forced') plt.show() diff --git a/doc/examples/plot_marching_cubes.py b/doc/examples/edges/plot_marching_cubes.py similarity index 100% rename from doc/examples/plot_marching_cubes.py rename to doc/examples/edges/plot_marching_cubes.py diff --git a/doc/examples/plot_medial_transform.py b/doc/examples/edges/plot_medial_transform.py similarity index 94% rename from doc/examples/plot_medial_transform.py rename to doc/examples/edges/plot_medial_transform.py index 8b9775cd..cc754268 100644 --- a/doc/examples/plot_medial_transform.py +++ b/doc/examples/edges/plot_medial_transform.py @@ -54,12 +54,13 @@ 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), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +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') ax2.contour(data, [0.5], colors='w') ax2.axis('off') -fig.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1) +fig.tight_layout() plt.show() diff --git a/doc/examples/plot_polygon.py b/doc/examples/edges/plot_polygon.py similarity index 100% rename from doc/examples/plot_polygon.py rename to doc/examples/edges/plot_polygon.py diff --git a/doc/examples/plot_shapes.py b/doc/examples/edges/plot_shapes.py similarity index 96% rename from doc/examples/plot_shapes.py rename to doc/examples/edges/plot_shapes.py index 34d1108c..0fd565f1 100644 --- a/doc/examples/plot_shapes.py +++ b/doc/examples/edges/plot_shapes.py @@ -5,16 +5,16 @@ Shapes This example shows how to draw several different shapes: - - line - - Bezier curve - - polygon - - circle - - ellipse +- line +- Bezier curve +- polygon +- circle +- ellipse Anti-aliased drawing for: - - line - - circle +- line +- circle """ import math diff --git a/doc/examples/plot_skeleton.py b/doc/examples/edges/plot_skeleton.py similarity index 87% rename from doc/examples/plot_skeleton.py rename to doc/examples/edges/plot_skeleton.py index 2bba3567..a406e297 100644 --- a/doc/examples/plot_skeleton.py +++ b/doc/examples/edges/plot_skeleton.py @@ -47,7 +47,9 @@ image[circle2] = 0 skeleton = skeletonize(image) # display results -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(8, 4.5), + sharex=True, sharey=True, + subplot_kw={'adjustable': 'box-forced'}) ax1.imshow(image, cmap=plt.cm.gray) ax1.axis('off') @@ -57,7 +59,6 @@ ax2.imshow(skeleton, cmap=plt.cm.gray) ax2.axis('off') ax2.set_title('skeleton', fontsize=20) -fig.subplots_adjust(wspace=0.02, hspace=0.02, top=0.98, - bottom=0.02, left=0.02, right=0.98) +fig.tight_layout() plt.show() diff --git a/doc/examples/features_detection/README.txt b/doc/examples/features_detection/README.txt new file mode 100644 index 00000000..82d1ab54 --- /dev/null +++ b/doc/examples/features_detection/README.txt @@ -0,0 +1,2 @@ +Detection of features and objects +--------------------------------- diff --git a/doc/examples/plot_blob.py b/doc/examples/features_detection/plot_blob.py similarity index 92% rename from doc/examples/plot_blob.py rename to doc/examples/features_detection/plot_blob.py index bb5dd089..c17cdf7d 100644 --- a/doc/examples/plot_blob.py +++ b/doc/examples/features_detection/plot_blob.py @@ -34,19 +34,20 @@ independent of the size of blobs as internally the implementation uses box filters instead of convolutions. Bright on dark as well as dark on bright blobs are detected. The downside is that small blobs (<3px) are not detected accurately. See :py:meth:`skimage.feature.blob_doh` for usage. - """ -from matplotlib import pyplot as plt from skimage import data from skimage.feature import blob_dog, blob_log, blob_doh from math import sqrt from skimage.color import rgb2gray +import matplotlib.pyplot as plt + image = data.hubble_deep_field()[0:500, 0:500] image_gray = rgb2gray(image) blobs_log = blob_log(image_gray, max_sigma=30, num_sigma=10, threshold=.1) + # Compute radii in the 3rd column. blobs_log[:, 2] = blobs_log[:, 2] * sqrt(2) @@ -61,14 +62,17 @@ titles = ['Laplacian of Gaussian', 'Difference of Gaussian', 'Determinant of Hessian'] sequence = zip(blobs_list, colors, titles) +fig, axes = plt.subplots(1, 3, figsize=(14, 4), sharex=True, sharey=True, + subplot_kw={'adjustable': 'box-forced'}) +plt.tight_layout() -fig,axes = plt.subplots(1, 3, sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) axes = axes.ravel() for blobs, color, title in sequence: ax = axes[0] axes = axes[1:] ax.set_title(title) ax.imshow(image, interpolation='nearest') + ax.set_axis_off() for blob in blobs: y, x, r = blob c = plt.Circle((x, y), r, color=color, linewidth=2, fill=False) diff --git a/doc/examples/plot_brief.py b/doc/examples/features_detection/plot_brief.py similarity index 100% rename from doc/examples/plot_brief.py rename to doc/examples/features_detection/plot_brief.py diff --git a/doc/examples/plot_censure.py b/doc/examples/features_detection/plot_censure.py similarity index 100% rename from doc/examples/plot_censure.py rename to doc/examples/features_detection/plot_censure.py diff --git a/doc/examples/plot_corner.py b/doc/examples/features_detection/plot_corner.py similarity index 100% rename from doc/examples/plot_corner.py rename to doc/examples/features_detection/plot_corner.py diff --git a/doc/examples/plot_daisy.py b/doc/examples/features_detection/plot_daisy.py similarity index 100% rename from doc/examples/plot_daisy.py rename to doc/examples/features_detection/plot_daisy.py diff --git a/doc/examples/plot_gabor.py b/doc/examples/features_detection/plot_gabor.py similarity index 100% rename from doc/examples/plot_gabor.py rename to doc/examples/features_detection/plot_gabor.py diff --git a/doc/examples/plot_gabors_from_astronaut.py b/doc/examples/features_detection/plot_gabors_from_astronaut.py similarity index 99% rename from doc/examples/plot_gabors_from_astronaut.py rename to doc/examples/features_detection/plot_gabors_from_astronaut.py index ed429203..a8cb50c2 100644 --- a/doc/examples/plot_gabors_from_astronaut.py +++ b/doc/examples/features_detection/plot_gabors_from_astronaut.py @@ -87,5 +87,5 @@ ax3.set_title("K-means filterbank (codebook)\non LGN-like DoG image") for ax in axes.ravel(): ax.axis('off') -fig.subplots_adjust(hspace=0.3) +fig.tight_layout() plt.show() diff --git a/doc/examples/plot_glcm.py b/doc/examples/features_detection/plot_glcm.py similarity index 100% rename from doc/examples/plot_glcm.py rename to doc/examples/features_detection/plot_glcm.py diff --git a/doc/examples/plot_hog.py b/doc/examples/features_detection/plot_hog.py similarity index 95% rename from doc/examples/plot_hog.py rename to doc/examples/features_detection/plot_hog.py index ab71bb36..e7b524cb 100644 --- a/doc/examples/plot_hog.py +++ b/doc/examples/features_detection/plot_hog.py @@ -15,11 +15,11 @@ Algorithm overview Compute a Histogram of Oriented Gradients (HOG) by - 1. (optional) global image normalisation - 2. computing the gradient image in x and y - 3. computing gradient histograms - 4. normalising across blocks - 5. flattening into a feature vector +1. (optional) global image normalisation +2. computing the gradient image in x and y +3. computing gradient histograms +4. normalising across blocks +5. flattening into a feature vector The first stage applies an optional global image normalisation equalisation that is designed to reduce the influence of illumination diff --git a/doc/examples/plot_holes_and_peaks.py b/doc/examples/features_detection/plot_holes_and_peaks.py similarity index 100% rename from doc/examples/plot_holes_and_peaks.py rename to doc/examples/features_detection/plot_holes_and_peaks.py diff --git a/doc/examples/plot_local_binary_pattern.py b/doc/examples/features_detection/plot_local_binary_pattern.py similarity index 100% rename from doc/examples/plot_local_binary_pattern.py rename to doc/examples/features_detection/plot_local_binary_pattern.py diff --git a/doc/examples/plot_multiblock_local_binary_pattern.py b/doc/examples/features_detection/plot_multiblock_local_binary_pattern.py similarity index 100% rename from doc/examples/plot_multiblock_local_binary_pattern.py rename to doc/examples/features_detection/plot_multiblock_local_binary_pattern.py diff --git a/doc/examples/plot_orb.py b/doc/examples/features_detection/plot_orb.py similarity index 100% rename from doc/examples/plot_orb.py rename to doc/examples/features_detection/plot_orb.py diff --git a/doc/examples/plot_template.py b/doc/examples/features_detection/plot_template.py similarity index 100% rename from doc/examples/plot_template.py rename to doc/examples/features_detection/plot_template.py diff --git a/doc/examples/plot_windowed_histogram.py b/doc/examples/features_detection/plot_windowed_histogram.py similarity index 100% rename from doc/examples/plot_windowed_histogram.py rename to doc/examples/features_detection/plot_windowed_histogram.py diff --git a/doc/examples/filters/README.txt b/doc/examples/filters/README.txt new file mode 100644 index 00000000..e7e55450 --- /dev/null +++ b/doc/examples/filters/README.txt @@ -0,0 +1,2 @@ +Filtering and restoration +------------------------- diff --git a/doc/examples/plot_denoise.py b/doc/examples/filters/plot_denoise.py similarity index 92% rename from doc/examples/plot_denoise.py rename to doc/examples/filters/plot_denoise.py index f39591d2..e66d2478 100644 --- a/doc/examples/plot_denoise.py +++ b/doc/examples/filters/plot_denoise.py @@ -38,7 +38,8 @@ 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), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5), sharex=True, + sharey=True, subplot_kw={'adjustable': 'box-forced'}) plt.gray() @@ -62,7 +63,6 @@ ax[1, 2].imshow(astro) ax[1, 2].axis('off') ax[1, 2].set_title('original') -fig.subplots_adjust(wspace=0.02, hspace=0.2, - top=0.9, bottom=0.05, left=0, right=1) +fig.tight_layout() plt.show() diff --git a/doc/examples/plot_entropy.py b/doc/examples/filters/plot_entropy.py similarity index 100% rename from doc/examples/plot_entropy.py rename to doc/examples/filters/plot_entropy.py diff --git a/doc/examples/plot_nonlocal_means.py b/doc/examples/filters/plot_nonlocal_means.py similarity index 88% rename from doc/examples/plot_nonlocal_means.py rename to doc/examples/filters/plot_nonlocal_means.py index 6082ee9c..d5c21b0e 100644 --- a/doc/examples/plot_nonlocal_means.py +++ b/doc/examples/filters/plot_nonlocal_means.py @@ -26,7 +26,8 @@ noisy = np.clip(noisy, 0, 1) denoise = denoise_nl_means(noisy, 7, 9, 0.08) -fig, ax = plt.subplots(ncols=2, figsize=(8, 4), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +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') @@ -35,7 +36,6 @@ ax[1].imshow(denoise) ax[1].axis('off') ax[1].set_title('non-local means') -fig.subplots_adjust(wspace=0.02, hspace=0.2, - top=0.9, bottom=0.05, left=0, right=1) +fig.tight_layout() plt.show() diff --git a/doc/examples/plot_phase_unwrap.py b/doc/examples/filters/plot_phase_unwrap.py similarity index 100% rename from doc/examples/plot_phase_unwrap.py rename to doc/examples/filters/plot_phase_unwrap.py diff --git a/doc/examples/plot_rank_mean.py b/doc/examples/filters/plot_rank_mean.py similarity index 71% rename from doc/examples/plot_rank_mean.py rename to doc/examples/filters/plot_rank_mean.py index 51f69788..a6fd75d6 100644 --- a/doc/examples/plot_rank_mean.py +++ b/doc/examples/filters/plot_rank_mean.py @@ -5,12 +5,12 @@ Mean filters This example compares the following mean filters of the rank filter package: - * **local mean**: all pixels belonging to the structuring element to compute - average gray level. - * **percentile mean**: only use values between percentiles p0 and p1 - (here 10% and 90%). - * **bilateral mean**: only use pixels of the structuring element having a gray - level situated inside g-s0 and g+s1 (here g-500 and g+500) +* **local mean**: all pixels belonging to the structuring element to compute + average gray level. +* **percentile mean**: only use values between percentiles p0 and p1 + (here 10% and 90%). +* **bilateral mean**: only use pixels of the structuring element having a gray + level situated inside g-s0 and g+s1 (here g-500 and g+500) Percentile and usual mean give here similar results, these filters smooth the complete image (background and details). Bilateral mean exhibits a high @@ -34,7 +34,8 @@ bilateral_result = rank.mean_bilateral(image, selem=selem, s0=500, s1=500) normal_result = rank.mean(image, selem=selem) -fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(8, 10), sharex=True, sharey=True) +fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(8, 10), + sharex=True, sharey=True) ax = axes.ravel() titles = ['Original', 'Percentile mean', 'Bilateral mean', 'Local mean'] diff --git a/doc/examples/plot_restoration.py b/doc/examples/filters/plot_restoration.py similarity index 88% rename from doc/examples/plot_restoration.py rename to doc/examples/filters/plot_restoration.py index 221a53b7..a22a68a8 100644 --- a/doc/examples/plot_restoration.py +++ b/doc/examples/filters/plot_restoration.py @@ -42,7 +42,9 @@ 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), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 5), + sharex=True, sharey=True, + subplot_kw={'adjustable': 'box-forced'}) plt.gray() @@ -54,7 +56,6 @@ ax[1].imshow(deconvolved) ax[1].axis('off') ax[1].set_title('Self tuned restoration') -fig.subplots_adjust(wspace=0.02, hspace=0.2, - top=0.9, bottom=0.05, left=0, right=1) +fig.tight_layout() plt.show() diff --git a/doc/examples/numpy_operations/README.txt b/doc/examples/numpy_operations/README.txt new file mode 100644 index 00000000..265f8a5f --- /dev/null +++ b/doc/examples/numpy_operations/README.txt @@ -0,0 +1,2 @@ +Operations on NumPy arrays +-------------------------- diff --git a/doc/examples/plot_camera_numpy.py b/doc/examples/numpy_operations/plot_camera_numpy.py similarity index 100% rename from doc/examples/plot_camera_numpy.py rename to doc/examples/numpy_operations/plot_camera_numpy.py diff --git a/doc/examples/plot_view_as_blocks.py b/doc/examples/numpy_operations/plot_view_as_blocks.py similarity index 95% rename from doc/examples/plot_view_as_blocks.py rename to doc/examples/numpy_operations/plot_view_as_blocks.py index 1a85bc41..53797aac 100644 --- a/doc/examples/plot_view_as_blocks.py +++ b/doc/examples/numpy_operations/plot_view_as_blocks.py @@ -49,8 +49,9 @@ 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, extent=(0, 128, 128, 0), interpolation='nearest', 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") @@ -65,5 +66,5 @@ ax3.set_title("Block view with\n local median pooling") ax3.imshow(median_view, interpolation='nearest', cmap=cm.Greys_r) ax3.set_axis_off() -fig.subplots_adjust(hspace=0.4, wspace=0.4) +fig.tight_layout() plt.show() diff --git a/doc/examples/segmentation/README.txt b/doc/examples/segmentation/README.txt new file mode 100644 index 00000000..c1d1fbeb --- /dev/null +++ b/doc/examples/segmentation/README.txt @@ -0,0 +1,2 @@ +Segmentation of objects +----------------------- diff --git a/doc/examples/plot_join_segmentations.py b/doc/examples/segmentation/plot_join_segmentations.py similarity index 93% rename from doc/examples/plot_join_segmentations.py rename to doc/examples/segmentation/plot_join_segmentations.py index 0ad7e57c..e1065a7d 100644 --- a/doc/examples/plot_join_segmentations.py +++ b/doc/examples/segmentation/plot_join_segmentations.py @@ -40,7 +40,8 @@ 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), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +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') @@ -58,5 +59,5 @@ axes[3].set_title('Join') for ax in axes: ax.axis('off') -fig.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1) +fig.tight_layout() plt.show() diff --git a/doc/examples/plot_label.py b/doc/examples/segmentation/plot_label.py similarity index 100% rename from doc/examples/plot_label.py rename to doc/examples/segmentation/plot_label.py diff --git a/doc/examples/plot_local_otsu.py b/doc/examples/segmentation/plot_local_otsu.py similarity index 62% rename from doc/examples/plot_local_otsu.py rename to doc/examples/segmentation/plot_local_otsu.py index d853c25a..8d7475d7 100644 --- a/doc/examples/plot_local_otsu.py +++ b/doc/examples/segmentation/plot_local_otsu.py @@ -10,23 +10,21 @@ structuring element. The example compares the local threshold with the global threshold. -.. note: local is much slower than global thresholding +.. Note: local is much slower than global thresholding .. [1] http://en.wikipedia.org/wiki/Otsu's_method """ -import matplotlib -import matplotlib.pyplot as plt from skimage import data from skimage.morphology import disk from skimage.filters import threshold_otsu, rank from skimage.util import img_as_ubyte +import matplotlib +import matplotlib.pyplot as plt matplotlib.rcParams['font.size'] = 9 - - img = img_as_ubyte(data.page()) radius = 15 @@ -36,26 +34,26 @@ local_otsu = rank.otsu(img, selem) threshold_global_otsu = threshold_otsu(img) global_otsu = img >= threshold_global_otsu +fig, ax = plt.subplots(2, 2, figsize=(8, 5), sharex=True, sharey=True, + subplot_kw={'adjustable': 'box-forced'}) +ax0, ax1, ax2, ax3 = ax.ravel() -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(ax0.imshow(img, cmap=plt.cm.gray), + ax=ax0, orientation='horizontal') +ax0.set_title('Original') +ax0.axis('off') -fig.colorbar(ax1.imshow(img, cmap=plt.cm.gray), +fig.colorbar(ax1.imshow(local_otsu, cmap=plt.cm.gray), ax=ax1, orientation='horizontal') -ax1.set_title('Original') +ax1.set_title('Local Otsu (radius=%d)' % radius) ax1.axis('off') -fig.colorbar(ax2.imshow(local_otsu, cmap=plt.cm.gray), - ax=ax2, orientation='horizontal') -ax2.set_title('Local Otsu (radius=%d)' % radius) +ax2.imshow(img >= local_otsu, cmap=plt.cm.gray) +ax2.set_title('Original >= Local Otsu' % threshold_global_otsu) ax2.axis('off') -ax3.imshow(img >= local_otsu, cmap=plt.cm.gray) -ax3.set_title('Original >= Local Otsu' % threshold_global_otsu) +ax3.imshow(global_otsu, cmap=plt.cm.gray) +ax3.set_title('Global Otsu (threshold = %d)' % threshold_global_otsu) ax3.axis('off') -ax4.imshow(global_otsu, cmap=plt.cm.gray) -ax4.set_title('Global Otsu (threshold = %d)' % threshold_global_otsu) -ax4.axis('off') - plt.show() diff --git a/doc/examples/plot_marked_watershed.py b/doc/examples/segmentation/plot_marked_watershed.py similarity index 100% rename from doc/examples/plot_marked_watershed.py rename to doc/examples/segmentation/plot_marked_watershed.py diff --git a/doc/examples/plot_ncut.py b/doc/examples/segmentation/plot_ncut.py similarity index 100% rename from doc/examples/plot_ncut.py rename to doc/examples/segmentation/plot_ncut.py diff --git a/doc/examples/plot_otsu.py b/doc/examples/segmentation/plot_otsu.py similarity index 100% rename from doc/examples/plot_otsu.py rename to doc/examples/segmentation/plot_otsu.py diff --git a/doc/examples/plot_peak_local_max.py b/doc/examples/segmentation/plot_peak_local_max.py similarity index 89% rename from doc/examples/plot_peak_local_max.py rename to doc/examples/segmentation/plot_peak_local_max.py index 8a690c5e..76a38b3b 100644 --- a/doc/examples/plot_peak_local_max.py +++ b/doc/examples/segmentation/plot_peak_local_max.py @@ -25,7 +25,8 @@ 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), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +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') @@ -41,7 +42,6 @@ ax3.plot(coordinates[:, 1], coordinates[:, 0], 'r.') ax3.axis('off') ax3.set_title('Peak local max') -fig.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, - bottom=0.02, left=0.02, right=0.98) +fig.tight_layout() plt.show() diff --git a/doc/examples/plot_rag.py b/doc/examples/segmentation/plot_rag.py similarity index 100% rename from doc/examples/plot_rag.py rename to doc/examples/segmentation/plot_rag.py diff --git a/doc/examples/plot_rag_draw.py b/doc/examples/segmentation/plot_rag_draw.py similarity index 100% rename from doc/examples/plot_rag_draw.py rename to doc/examples/segmentation/plot_rag_draw.py diff --git a/doc/examples/plot_rag_mean_color.py b/doc/examples/segmentation/plot_rag_mean_color.py similarity index 100% rename from doc/examples/plot_rag_mean_color.py rename to doc/examples/segmentation/plot_rag_mean_color.py diff --git a/doc/examples/plot_rag_merge.py b/doc/examples/segmentation/plot_rag_merge.py similarity index 100% rename from doc/examples/plot_rag_merge.py rename to doc/examples/segmentation/plot_rag_merge.py diff --git a/doc/examples/plot_random_walker_segmentation.py b/doc/examples/segmentation/plot_random_walker_segmentation.py similarity index 91% rename from doc/examples/plot_random_walker_segmentation.py rename to doc/examples/segmentation/plot_random_walker_segmentation.py index 81e0bd91..da050f8f 100644 --- a/doc/examples/plot_random_walker_segmentation.py +++ b/doc/examples/segmentation/plot_random_walker_segmentation.py @@ -38,7 +38,8 @@ 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), sharex=True, sharey=True) +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') @@ -52,6 +53,5 @@ 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, - right=1) +fig.tight_layout() plt.show() diff --git a/doc/examples/plot_regionprops.py b/doc/examples/segmentation/plot_regionprops.py similarity index 100% rename from doc/examples/plot_regionprops.py rename to doc/examples/segmentation/plot_regionprops.py diff --git a/doc/examples/plot_segmentations.py b/doc/examples/segmentation/plot_segmentations.py similarity index 96% rename from doc/examples/plot_segmentations.py rename to doc/examples/segmentation/plot_segmentations.py index af601f53..33008069 100644 --- a/doc/examples/plot_segmentations.py +++ b/doc/examples/segmentation/plot_segmentations.py @@ -79,9 +79,10 @@ 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, sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +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) +fig.tight_layout() ax[0].imshow(mark_boundaries(img, segments_fz)) ax[0].set_title("Felzenszwalbs's method") diff --git a/doc/examples/plot_threshold_adaptive.py b/doc/examples/segmentation/plot_threshold_adaptive.py similarity index 100% rename from doc/examples/plot_threshold_adaptive.py rename to doc/examples/segmentation/plot_threshold_adaptive.py diff --git a/doc/examples/plot_watershed.py b/doc/examples/segmentation/plot_watershed.py similarity index 93% rename from doc/examples/plot_watershed.py rename to doc/examples/segmentation/plot_watershed.py index 664c4c77..fb6f8a9b 100644 --- a/doc/examples/plot_watershed.py +++ b/doc/examples/segmentation/plot_watershed.py @@ -48,7 +48,8 @@ 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), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +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') @@ -61,6 +62,5 @@ ax2.set_title('Separated objects') for ax in axes: ax.axis('off') -fig.subplots_adjust(hspace=0.01, wspace=0.01, top=0.9, bottom=0, left=0, - right=1) +fig.tight_layout() plt.show() diff --git a/doc/examples/transform/README.txt b/doc/examples/transform/README.txt new file mode 100644 index 00000000..dcc361fb --- /dev/null +++ b/doc/examples/transform/README.txt @@ -0,0 +1,2 @@ +Geometrical transformations and registration +-------------------------------------------- diff --git a/doc/examples/plot_edge_modes.py b/doc/examples/transform/plot_edge_modes.py similarity index 70% rename from doc/examples/plot_edge_modes.py rename to doc/examples/transform/plot_edge_modes.py index 8689e4a8..9c583006 100644 --- a/doc/examples/plot_edge_modes.py +++ b/doc/examples/transform/plot_edge_modes.py @@ -7,10 +7,10 @@ This example illustrates the different edge modes available during interpolation in routines such as `skimage.transform.rescale` and `skimage.transform.resize`. """ -from skimage._shared.interpolation import extend_image -import skimage.data -import matplotlib.pyplot as plt import numpy as np +import matplotlib.pyplot as plt + +from skimage._shared.interpolation import extend_image img = np.zeros((16, 16)) img[:8, :8] += 1 @@ -20,18 +20,19 @@ img[:1, :1] += 2 img[8, 8] = 4 modes = ['constant', 'edge', 'wrap', 'reflect', 'symmetric'] -fig, axes = plt.subplots(1, 5, figsize=(15, 5)) +fig, axes = plt.subplots(2, 3) +axes = axes.flatten() + for n, mode in enumerate(modes): img_extended = extend_image(img, pad=img.shape[0], mode=mode) axes[n].imshow(img_extended, cmap=plt.cm.gray, interpolation='nearest') - axes[n].plot([15.5, 15.5], [15.5, 31.5], 'y--', linewidth=0.5) - axes[n].plot([31.5, 31.5], [15.5, 31.5], 'y--', linewidth=0.5) - axes[n].plot([15.5, 31.5], [15.5, 15.5], 'y--', linewidth=0.5) - axes[n].plot([15.5, 31.5], [31.5, 31.5], 'y--', linewidth=0.5) - axes[n].set_axis_off() - axes[n].set_aspect('equal') + axes[n].plot([15.5, 15.5, 31.5, 31.5, 15.5], + [15.5, 31.5, 31.5, 15.5, 15.5], 'y--', linewidth=0.5) axes[n].set_title(mode) +for n in range(len(axes)): + axes[n].set_axis_off() + axes[n].set_aspect('equal') + plt.tight_layout() - plt.show() diff --git a/doc/examples/plot_matching.py b/doc/examples/transform/plot_matching.py similarity index 100% rename from doc/examples/plot_matching.py rename to doc/examples/transform/plot_matching.py diff --git a/doc/examples/plot_piecewise_affine.py b/doc/examples/transform/plot_piecewise_affine.py similarity index 100% rename from doc/examples/plot_piecewise_affine.py rename to doc/examples/transform/plot_piecewise_affine.py diff --git a/doc/examples/plot_pyramid.py b/doc/examples/transform/plot_pyramid.py similarity index 100% rename from doc/examples/plot_pyramid.py rename to doc/examples/transform/plot_pyramid.py diff --git a/doc/examples/plot_radon_transform.py b/doc/examples/transform/plot_radon_transform.py similarity index 94% rename from doc/examples/plot_radon_transform.py rename to doc/examples/transform/plot_radon_transform.py index cc327d72..541abf99 100644 --- a/doc/examples/plot_radon_transform.py +++ b/doc/examples/transform/plot_radon_transform.py @@ -31,9 +31,9 @@ Technique (SART). For further information on tomographic reconstruction, see - - AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging", - http://www.slaney.org/pct/pct-toc.html - - http://en.wikipedia.org/wiki/Radon_transform +- AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging", + http://www.slaney.org/pct/pct-toc.html +- http://en.wikipedia.org/wiki/Radon_transform The forward transform ===================== @@ -73,7 +73,7 @@ ax2.set_ylabel("Projection position (pixels)") ax2.imshow(sinogram, cmap=plt.cm.Greys_r, extent=(0, 180, 0, sinogram.shape[0]), aspect='auto') -fig.subplots_adjust(hspace=0.4, wspace=0.5) +fig.tight_layout() plt.show() """ @@ -101,7 +101,9 @@ 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), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +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 +154,8 @@ 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), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +fig, ax = plt.subplots(2, 2, figsize=(8, 8.5), sharex=True, sharey=True, + subplot_kw={'adjustable': 'box-forced'}) ax1, ax2, ax3, ax4 = ax.ravel() ax1.set_title("Reconstruction\nSART") ax1.imshow(reconstruction_sart, cmap=plt.cm.Greys_r) diff --git a/doc/examples/plot_ransac.py b/doc/examples/transform/plot_ransac.py similarity index 100% rename from doc/examples/plot_ransac.py rename to doc/examples/transform/plot_ransac.py diff --git a/doc/examples/plot_ransac3D.py b/doc/examples/transform/plot_ransac3D.py similarity index 100% rename from doc/examples/plot_ransac3D.py rename to doc/examples/transform/plot_ransac3D.py diff --git a/doc/examples/plot_register_translation.py b/doc/examples/transform/plot_register_translation.py similarity index 100% rename from doc/examples/plot_register_translation.py rename to doc/examples/transform/plot_register_translation.py diff --git a/doc/examples/plot_seam_carving.py b/doc/examples/transform/plot_seam_carving.py similarity index 100% rename from doc/examples/plot_seam_carving.py rename to doc/examples/transform/plot_seam_carving.py diff --git a/doc/examples/plot_ssim.py b/doc/examples/transform/plot_ssim.py similarity index 86% rename from doc/examples/plot_ssim.py rename to doc/examples/transform/plot_ssim.py index 112a6971..1a206653 100644 --- a/doc/examples/plot_ssim.py +++ b/doc/examples/transform/plot_ssim.py @@ -19,19 +19,14 @@ but with very different mean structural similarity indices. assessment: From error visibility to structural similarity," IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600-612, Apr. 2004. - """ + import numpy as np -import matplotlib import matplotlib.pyplot as plt from skimage import data, img_as_float from skimage.measure import structural_similarity as ssim - -matplotlib.rcParams['font.size'] = 9 - - img = img_as_float(data.camera()) rows, cols = img.shape @@ -45,7 +40,10 @@ 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), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'}) +fig, (ax0, ax1, ax2) = plt.subplots(nrows=1, ncols=3, figsize=(16, 6), + sharex=True, sharey=True, + subplot_kw={'adjustable': 'box-forced'}) +plt.tight_layout() mse_none = mse(img, img) ssim_none = ssim(img, img, dynamic_range=img.max() - img.min()) @@ -63,13 +61,16 @@ label = 'MSE: %2.f, SSIM: %.2f' ax0.imshow(img, cmap=plt.cm.gray, vmin=0, vmax=1) ax0.set_xlabel(label % (mse_none, ssim_none)) ax0.set_title('Original image') +ax0.axes.get_yaxis().set_visible(False) ax1.imshow(img_noise, cmap=plt.cm.gray, vmin=0, vmax=1) ax1.set_xlabel(label % (mse_noise, ssim_noise)) ax1.set_title('Image with noise') +ax1.axes.get_yaxis().set_visible(False) ax2.imshow(img_const, cmap=plt.cm.gray, vmin=0, vmax=1) ax2.set_xlabel(label % (mse_const, ssim_const)) ax2.set_title('Image plus constant') +ax2.axes.get_yaxis().set_visible(False) plt.show() diff --git a/doc/examples/plot_swirl.py b/doc/examples/transform/plot_swirl.py similarity index 100% rename from doc/examples/plot_swirl.py rename to doc/examples/transform/plot_swirl.py diff --git a/doc/examples/applications/README.txt b/doc/examples/xx_applications/README.txt similarity index 100% rename from doc/examples/applications/README.txt rename to doc/examples/xx_applications/README.txt diff --git a/doc/examples/applications/plot_coins_segmentation.py b/doc/examples/xx_applications/plot_coins_segmentation.py similarity index 100% rename from doc/examples/applications/plot_coins_segmentation.py rename to doc/examples/xx_applications/plot_coins_segmentation.py diff --git a/doc/examples/applications/plot_geometric.py b/doc/examples/xx_applications/plot_geometric.py similarity index 100% rename from doc/examples/applications/plot_geometric.py rename to doc/examples/xx_applications/plot_geometric.py diff --git a/doc/examples/applications/plot_morphology.py b/doc/examples/xx_applications/plot_morphology.py similarity index 100% rename from doc/examples/applications/plot_morphology.py rename to doc/examples/xx_applications/plot_morphology.py diff --git a/doc/examples/applications/plot_rank_filters.py b/doc/examples/xx_applications/plot_rank_filters.py similarity index 100% rename from doc/examples/applications/plot_rank_filters.py rename to doc/examples/xx_applications/plot_rank_filters.py diff --git a/doc/release/contribs.py b/doc/release/contribs.py index c0a34e6e..6fcc0d87 100755 --- a/doc/release/contribs.py +++ b/doc/release/contribs.py @@ -10,41 +10,45 @@ if len(sys.argv) != 2: tag = sys.argv[1] -def call(cmd): - return subprocess.check_output(shlex.split(cmd), universal_newlines=True).split('\n') -tag_date = call("git show --format='%%ci' %s" % tag)[0] -print("Release %s was on %s\n" % (tag, tag_date)) +if not sys.version_info[:2] == (2, 6): -merges = call("git log --since='%s' --merges --format='>>>%%B' --reverse" % tag_date) -merges = [m for m in merges if m.strip()] -merges = '\n'.join(merges).split('>>>') -merges = [m.split('\n')[:2] for m in merges] -merges = [m for m in merges if len(m) == 2 and m[1].strip()] + def call(cmd): + return subprocess.check_output(shlex.split(cmd), universal_newlines=True).split('\n') -num_commits = call("git rev-list %s..HEAD --count" % tag)[0] -print("A total of %s changes have been committed.\n" % num_commits) + tag_date = call("git show --format='%%ci' %s" % tag)[0] + print("Release %s was on %s\n" % (tag, tag_date)) -print("It contained the following %d merges:\n" % len(merges)) -for (merge, message) in merges: - if merge.startswith('Merge pull request #'): - PR = ' (%s)' % merge.split()[3] - else: - PR = '' + merges = call("git log --since='%s' --merges --format='>>>%%B' --reverse" % tag_date) + merges = [m for m in merges if m.strip()] + merges = '\n'.join(merges).split('>>>') + merges = [m.split('\n')[:2] for m in merges] + merges = [m for m in merges if len(m) == 2 and m[1].strip()] - print('- ' + message + PR) + num_commits = call("git rev-list %s..HEAD --count" % tag)[0] + print("A total of %s changes have been committed.\n" % num_commits) + + print("It contained the following %d merges:\n" % len(merges)) + for (merge, message) in merges: + if merge.startswith('Merge pull request #'): + PR = ' (%s)' % merge.split()[3] + else: + PR = '' + + print('- ' + message + PR) -print("\nMade by the following committers [alphabetical by last name]:\n") + print("\nMade by the following committers [alphabetical by last name]:\n") -authors = call("git log --since='%s' --format=%%aN" % tag_date) -authors = [a.strip() for a in authors if a.strip()] + authors = call("git log --since='%s' --format=%%aN" % tag_date) + authors = [a.strip() for a in authors if a.strip()] -def key(author): - author = [v for v in author.split() if v[0] in string.ascii_letters] - return author[-1] + def key(author): + author = [v for v in author.split() if v[0] in string.ascii_letters] + if len(author) > 0: + return author[-1] -authors = sorted(set(authors), key=key) + authors = sorted(set(authors), key=key) -for a in authors: - print('- ' + a) + for a in authors: + print('- ' + a) diff --git a/doc/source/user_guide/numpy_images.txt b/doc/source/user_guide/numpy_images.txt index 46f21afa..f4b3cf42 100644 --- a/doc/source/user_guide/numpy_images.txt +++ b/doc/source/user_guide/numpy_images.txt @@ -81,9 +81,9 @@ disk: :: ... (nrows / 2)**2) >>> camera[outer_disk_mask] = 0 -.. image:: ../auto_examples/images/plot_camera_numpy_1.png +.. image:: ../auto_examples/numpy_operations/images/plot_camera_numpy_1.png :width: 45% - :target: ../auto_examples/plot_camera_numpy.html + :target: ../auto_examples/numpy_operations/plot_camera_numpy.html Boolean arithmetic can be used to define more complex masks: :: diff --git a/doc/source/user_guide/transforming_image_data.txt b/doc/source/user_guide/transforming_image_data.txt index 9dfe3139..472b55dc 100644 --- a/doc/source/user_guide/transforming_image_data.txt +++ b/doc/source/user_guide/transforming_image_data.txt @@ -78,8 +78,8 @@ using an array of labels to encode the regions to be represented with the same color. -.. image: ../auto_examples/images/plot_join_segmentations_1.png - :target: ../auto_examples/plot_join_segmentations.html +.. image: ../auto_examples/segmentation/images/plot_join_segmentations_1.png + :target: ../auto_examples/segmentation/plot_join_segmentations.html :align: center :width: 80% @@ -87,9 +87,9 @@ same color. .. topic:: Examples: - * :ref:`example_plot_tinting_grayscale_images.py` - * :ref:`example_plot_join_segmentations.py` - * :ref:`example_plot_rag_mean_color.py` + * :ref:`example_color_exposure_plot_tinting_grayscale_images.py` + * :ref:`example_segmentation_plot_join_segmentations.py` + * :ref:`example_segmentation_plot_rag_mean_color.py` Contrast and exposure @@ -122,7 +122,7 @@ the image. The histogram of pixel values is computed with :func:`histogram` returns the number of pixels for each value bin, and the centers of the bins. The behavior of :func:`histogram` is therefore -slightly different from the one of :func:`np.histogram`, which returns +slightly different from the one of :func:`numpy.histogram`, which returns the boundaries of the bins. The simplest contrast enhancement :func:`rescale_intensity` consists in @@ -157,16 +157,16 @@ details are enhanced in large regions with poor contrast. As a further refinement, histogram equalization can be performed in subregions of the image with :func:`equalize_adapthist`, in order to correct for exposure gradients across the image. See the example -:ref:`example_plot_equalize.py`. +:ref:`example_color_exposure_plot_equalize.py`. -.. image:: ../auto_examples/images/plot_equalize_1.png - :target: ../auto_examples/plot_equalize.html +.. image:: ../auto_examples/color_exposure/images/plot_equalize_1.png + :target: ../auto_examples/color_exposure/plot_equalize.html :align: center :width: 90% .. topic:: Examples: - * :ref:`example_plot_equalize.py` + * :ref:`example_color_exposure_plot_equalize.py` diff --git a/doc/source/user_guide/tutorial_segmentation.txt b/doc/source/user_guide/tutorial_segmentation.txt index 76af842e..e5382528 100644 --- a/doc/source/user_guide/tutorial_segmentation.txt +++ b/doc/source/user_guide/tutorial_segmentation.txt @@ -11,8 +11,8 @@ the coins cannot be done directly from the histogram of grey values, because the background shares enough grey levels with the coins that a thresholding segmentation is not sufficient. -.. image:: ../auto_examples/applications/images/plot_coins_segmentation_1.png - :target: ../auto_examples/applications/plot_coins_segmentation.html +.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_1.png + :target: ../auto_examples/xx_applications/plot_coins_segmentation.html :align: center :: @@ -26,8 +26,8 @@ Simply thresholding the image leads either to missing significant parts of the coins, or to merging parts of the background with the coins. This is due to the inhomogeneous lighting of the image. -.. image:: ../auto_examples/applications/images/plot_coins_segmentation_2.png - :target: ../auto_examples/applications/plot_coins_segmentation.html +.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_2.png + :target: ../auto_examples/xx_applications/plot_coins_segmentation.html :align: center A first idea is to take advantage of the local contrast, that is, to @@ -53,8 +53,8 @@ boundary of the coins, or inside the coins. >>> from scipy import ndimage as ndi >>> fill_coins = ndi.binary_fill_holes(edges) -.. image:: ../auto_examples/applications/images/plot_coins_segmentation_3.png - :target: ../auto_examples/applications/plot_coins_segmentation.html +.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_3.png + :target: ../auto_examples/xx_applications/plot_coins_segmentation.html :align: center Now that we have contours that delineate the outer boundary of the coins, @@ -62,8 +62,8 @@ we fill the inner part of the coins using the ``ndi.binary_fill_holes`` function, which uses mathematical morphology to fill the holes. -.. image:: ../auto_examples/applications/images/plot_coins_segmentation_4.png - :target: ../auto_examples/applications/plot_coins_segmentation.html +.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_4.png + :target: ../auto_examples/xx_applications/plot_coins_segmentation.html :align: center Most coins are well segmented out of the background. Small objects from @@ -83,8 +83,8 @@ has not been segmented correctly at all. The reason is that the contour that we got from the Canny detector was not completely closed, therefore the filling function did not fill the inner part of the coin. -.. image:: ../auto_examples/applications/images/plot_coins_segmentation_5.png - :target: ../auto_examples/applications/plot_coins_segmentation.html +.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_5.png + :target: ../auto_examples/xx_applications/plot_coins_segmentation.html :align: center Therefore, this segmentation method is not very robust: if we miss a @@ -128,8 +128,8 @@ separate the coins from the background. and here is the corresponding 2-D plot: -.. image:: ../auto_examples/applications/images/plot_coins_segmentation_6.png - :target: ../auto_examples/applications/plot_coins_segmentation.html +.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_6.png + :target: ../auto_examples/xx_applications/plot_coins_segmentation.html :align: center The next step is to find markers of the background and the coins based on the @@ -139,8 +139,8 @@ extreme parts of the histogram of grey values:: >>> markers[coins < 30] = 1 >>> markers[coins > 150] = 2 -.. image:: ../auto_examples/applications/images/plot_coins_segmentation_7.png - :target: ../auto_examples/applications/plot_coins_segmentation.html +.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_7.png + :target: ../auto_examples/xx_applications/plot_coins_segmentation.html :align: center Let us now compute the watershed transform:: @@ -148,8 +148,8 @@ Let us now compute the watershed transform:: >>> from skimage.morphology import watershed >>> segmentation = watershed(elevation_map, markers) -.. image:: ../auto_examples/applications/images/plot_coins_segmentation_8.png - :target: ../auto_examples/applications/plot_coins_segmentation.html +.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_8.png + :target: ../auto_examples/xx_applications/plot_coins_segmentation.html :align: center With this method, the result is satisfying for all coins. Even if the @@ -165,7 +165,7 @@ We can now label all the coins one by one using ``ndi.label``:: >>> labeled_coins, _ = ndi.label(segmentation) -.. image:: ../auto_examples/applications/images/plot_coins_segmentation_9.png - :target: ../auto_examples/applications/plot_coins_segmentation.html +.. image:: ../auto_examples/xx_applications/images/plot_coins_segmentation_9.png + :target: ../auto_examples/xx_applications/plot_coins_segmentation.html :align: center diff --git a/requirements.txt b/requirements.txt index 8b97dd4b..19c541f5 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,5 +3,5 @@ numpy>=1.6.1 scipy>=0.9.0 six>=1.4 networkx>=1.8 -pillow>=1.7.8 +pillow>=2.1.0 dask[array]>=0.5.0 diff --git a/skimage/__init__.py b/skimage/__init__.py index 2b152156..673562a7 100644 --- a/skimage/__init__.py +++ b/skimage/__init__.py @@ -156,4 +156,9 @@ else: _raise_build_error(e) from .util.dtype import * + +if sys.version.startswith('2.6'): + warnings.warn("Python 2.6 is deprecated and will not be supported in scikit-image 0.13+") + + del warnings, functools, osp, imp, sys diff --git a/skimage/data/foo3x5x4indexed.png b/skimage/data/foo3x5x4indexed.png new file mode 100644 index 00000000..cc969d03 Binary files /dev/null and b/skimage/data/foo3x5x4indexed.png differ diff --git a/skimage/data/lena_GRAY_U8_hog.npy b/skimage/data/lena_GRAY_U8_hog.npy new file mode 100644 index 00000000..71ebf61f Binary files /dev/null and b/skimage/data/lena_GRAY_U8_hog.npy differ diff --git a/skimage/draw/draw.py b/skimage/draw/draw.py index 9a816810..c91aca2c 100644 --- a/skimage/draw/draw.py +++ b/skimage/draw/draw.py @@ -5,19 +5,19 @@ from ._draw import _coords_inside_image def _ellipse_in_shape(shape, center, radiuses): """Generate coordinates of points within ellipse bounded by shape.""" - y, x = np.ogrid[0:float(shape[0]), 0:float(shape[1])] - cy, cx = center + r_lim, c_lim = np.ogrid[0:float(shape[0]), 0:float(shape[1])] + r, c = center ry, rx = radiuses - distances = ((y - cy) / ry) ** 2 + ((x - cx) / rx) ** 2 + distances = ((r_lim - r) / ry) ** 2 + ((c_lim - c) / rx) ** 2 return np.nonzero(distances < 1) -def ellipse(cy, cx, yradius, xradius, shape=None): +def ellipse(r, c, yradius, xradius, shape=None): """Generate coordinates of pixels within ellipse. Parameters ---------- - cy, cx : double + r, c : double Centre coordinate of ellipse. yradius, xradius : double Minor and major semi-axes. ``(x/xradius)**2 + (y/yradius)**2 = 1``. @@ -53,7 +53,7 @@ def ellipse(cy, cx, yradius, xradius, shape=None): """ - center = np.array([cy, cx]) + center = np.array([r, c]) radiuses = np.array([yradius, xradius]) # The upper_left and lower_right corners of the @@ -77,12 +77,12 @@ def ellipse(cy, cx, yradius, xradius, shape=None): return rr, cc -def circle(cy, cx, radius, shape=None): +def circle(r, c, radius, shape=None): """Generate coordinates of pixels within circle. Parameters ---------- - cy, cx : double + r, c : double Centre coordinate of circle. radius: double Radius of circle. @@ -122,7 +122,7 @@ def circle(cy, cx, radius, shape=None): """ - return ellipse(cy, cx, radius, radius, shape) + return ellipse(r, c, radius, radius, shape) def set_color(img, coords, color): diff --git a/skimage/feature/_hog.py b/skimage/feature/_hog.py index afb2b4e9..760102e7 100644 --- a/skimage/feature/_hog.py +++ b/skimage/feature/_hog.py @@ -5,7 +5,8 @@ from . import _hoghistogram def hog(image, orientations=9, pixels_per_cell=(8, 8), - cells_per_block=(3, 3), visualise=False, normalise=False): + cells_per_block=(3, 3), visualise=False, normalise=False, + feature_vector=True): """Extract Histogram of Oriented Gradients (HOG) for a given image. Compute a Histogram of Oriented Gradients (HOG) by @@ -31,6 +32,9 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), normalise : bool, optional Apply power law compression to normalise the image before processing. + feature_vector : bool, optional + Return the data as a feature vector by calling .ravel() on the result + just before returning. Returns ------- @@ -127,13 +131,11 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), orientations_arr = np.arange(orientations) dx_arr = radius * np.cos(orientations_arr / orientations * np.pi) dy_arr = radius * np.sin(orientations_arr / orientations * np.pi) - cr2 = cy + cy - cc2 = cx + cx hog_image = np.zeros((sy, sx), dtype=float) for x in range(n_cellsx): for y in range(n_cellsy): for o, dx, dy in zip(orientations_arr, dx_arr, dy_arr): - centre = tuple([y * cr2 // 2, x * cc2 // 2]) + centre = tuple([y * cy + cy // 2, x * cx + cx // 2]) rr, cc = draw.line(int(centre[0] - dx), int(centre[1] + dy), int(centre[0] + dx), @@ -171,8 +173,11 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), overlapping grid of blocks covering the detection window into a combined feature vector for use in the window classifier. """ + + if feature_vector: + normalised_blocks = normalised_blocks.ravel() if visualise: - return normalised_blocks.ravel(), hog_image + return normalised_blocks, hog_image else: - return normalised_blocks.ravel() + return normalised_blocks diff --git a/skimage/feature/_hoghistogram.pyx b/skimage/feature/_hoghistogram.pyx index 3f8899dd..cd34caae 100644 --- a/skimage/feature/_hoghistogram.pyx +++ b/skimage/feature/_hoghistogram.pyx @@ -10,7 +10,9 @@ cdef float cell_hog(double[:, ::1] magnitude, float orientation_start, float orientation_end, int cell_columns, int cell_rows, int column_index, int row_index, - int size_columns, int size_rows) nogil: + int size_columns, int size_rows, + int range_rows_start, int range_rows_stop, + int range_columns_start, int range_columns_stop) nogil: """Calculation of the cell's HOG value Parameters @@ -35,22 +37,23 @@ cdef float cell_hog(double[:, ::1] magnitude, Number of columns. size_rows : int Number of rows. + range_rows_start : int + Start row of cell. + range_rows_stop : int + Stop row of cell. + range_columns_start : int + Start column of cell. + range_columns_stop : int + Stop column of cell Returns ------- total : float The total HOG value. """ - cdef int cell_column, cell_row, cell_row_index, cell_column_index, \ - range_columns_start, range_columns_stop, range_rows_start, \ - range_rows_stop - - range_rows_stop = cell_rows/2 - range_rows_start = -range_rows_stop - range_columns_stop = cell_columns/2 - range_columns_start = -range_columns_stop - + cdef int cell_column, cell_row, cell_row_index, cell_column_index cdef float total = 0. + for cell_row in range(range_rows_start, range_rows_stop): cell_row_index = row_index + cell_row if (cell_row_index < 0 or cell_row_index >= size_rows): @@ -67,7 +70,7 @@ cdef float cell_hog(double[:, ::1] magnitude, total += magnitude[cell_row_index, cell_column_index] - return total + return total / (cell_rows * cell_columns) def hog_histograms(double[:, ::1] gradient_columns, double[:, ::1] gradient_rows, @@ -106,7 +109,9 @@ def hog_histograms(double[:, ::1] gradient_columns, gradient_rows) cdef double[:, ::1] orientation = \ np.arctan2(gradient_rows, gradient_columns) * (180 / np.pi) % 180 - cdef int i, x, y, o, yi, xi, cc, cr, x0, y0 + cdef int i, x, y, o, yi, xi, cc, cr, x0, y0, \ + range_rows_start, range_rows_stop, \ + range_columns_start, range_columns_stop cdef float orientation_start, orientation_end, \ number_of_orientations_per_180 @@ -115,6 +120,10 @@ def hog_histograms(double[:, ::1] gradient_columns, cc = cell_rows * number_of_cells_rows cr = cell_columns * number_of_cells_columns number_of_orientations_per_180 = 180. / number_of_orientations + range_rows_stop = cell_rows/2 + range_rows_start = -range_rows_stop + range_columns_stop = cell_columns/2 + range_columns_start = -range_columns_stop with nogil: # compute orientations integral images @@ -134,7 +143,8 @@ def hog_histograms(double[:, ::1] gradient_columns, while x < cr: orientation_histogram[yi, xi, i] = cell_hog(magnitude, orientation, orientation_start, orientation_end, - cell_columns, cell_rows, x, y, size_columns, size_rows) + cell_columns, cell_rows, x, y, size_columns, size_rows, + range_rows_start, range_rows_stop, range_columns_start, range_columns_stop) xi += 1 x += cell_columns diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx index b90ab6dd..57d8ec81 100644 --- a/skimage/feature/_texture.pyx +++ b/skimage/feature/_texture.pyx @@ -271,8 +271,8 @@ def _local_binary_pattern(double[:, ::1] image, # Values represent offsets of neighbour rectangles relative to central one. # It has order starting from top left and going clockwise. cdef: - Py_ssize_t[::1] mlbp_r_offsets = np.asarray([-1, -1, -1, 0, 1, 1, 1, 0]) - Py_ssize_t[::1] mlbp_c_offsets = np.asarray([-1, 0, 1, 1, 1, 0, -1, -1]) + Py_ssize_t[::1] mlbp_r_offsets = np.asarray([-1, -1, -1, 0, 1, 1, 1, 0], dtype=np.intp) + Py_ssize_t[::1] mlbp_c_offsets = np.asarray([-1, 0, 1, 1, 1, 0, -1, -1], dtype=np.intp) def _multiblock_lbp(float[:, ::1] int_image, diff --git a/skimage/feature/tests/test_hog.py b/skimage/feature/tests/test_hog.py index 6112d726..3afc0d89 100644 --- a/skimage/feature/tests/test_hog.py +++ b/skimage/feature/tests/test_hog.py @@ -1,5 +1,7 @@ +import os import numpy as np from scipy import ndimage as ndi +import skimage as si from skimage import data from skimage import feature from skimage import img_as_float @@ -9,7 +11,7 @@ from numpy.testing import (assert_raises, ) -def test_histogram_of_oriented_gradients(): +def test_histogram_of_oriented_gradients_output_size(): img = img_as_float(data.astronaut()[:256, :].mean(axis=2)) fd = feature.hog(img, orientations=9, pixels_per_cell=(8, 8), @@ -18,6 +20,17 @@ def test_histogram_of_oriented_gradients(): assert len(fd) == 9 * (256 // 8) * (512 // 8) +def test_histogram_of_oriented_gradients_output_correctness(): + img = np.load(os.path.join(si.data_dir, 'lena_GRAY_U8.npy')) + correct_output = np.load(os.path.join(si.data_dir, 'lena_GRAY_U8_hog.npy')) + + output = feature.hog(img, orientations=9, pixels_per_cell=(8, 8), + cells_per_block=(3, 3), feature_vector=True, + normalise=False, visualise=False) + + assert_almost_equal(output, correct_output) + + def test_hog_image_size_cell_size_mismatch(): image = data.camera()[:150, :200] fd = feature.hog(image, orientations=9, pixels_per_cell=(8, 8), diff --git a/skimage/future/graph/rag.py b/skimage/future/graph/rag.py index 2975a027..6a45447b 100644 --- a/skimage/future/graph/rag.py +++ b/skimage/future/graph/rag.py @@ -1,16 +1,6 @@ -try: - import networkx as nx -except ImportError: - msg = "Graph functions require networkx, which is not installed" - - class nx: - class Graph: - def __init__(self, *args, **kwargs): - raise ImportError(msg) - import warnings - warnings.warn(msg) - +import networkx as nx import numpy as np +from numpy.lib.stride_tricks import as_strided from scipy import ndimage as ndi import math from ... import draw, measure, segmentation, util, color @@ -51,21 +41,79 @@ def min_weight(graph, src, dst, n): return min(w1, w2) +def _add_edge_filter(values, graph): + """Create edge in `graph` between central element of `values` and the rest. + + Add an edge between the middle element in `values` and + all other elements of `values` into `graph`. ``values[len(values) // 2]`` + is expected to be the central value of the footprint used. + + Parameters + ---------- + values : array + The array to process. + graph : RAG + The graph to add edges in. + + Returns + ------- + 0 : float + Always returns 0. The return value is required so that `generic_filter` + can put it in the output array, but it is ignored by this filter. + """ + values = values.astype(int) + center = values[len(values) // 2] + for value in values: + if value != center and not graph.has_edge(center, value): + graph.add_edge(center, value) + return 0. + + class RAG(nx.Graph): """ The Region Adjacency Graph (RAG) of an image, subclasses `networx.Graph `_ + + Parameters + ---------- + label_image : array of int + An initial segmentation, with each region labeled as a different + integer. Every unique value in ``label_image`` will correspond to + a node in the graph. + connectivity : int in {1, ..., ``label_image.ndim``}, optional + The connectivity between pixels in ``label_image``. For a 2D image, + a connectivity of 1 corresponds to immediate neighbors up, down, + left, and right, while a connectivity of 2 also includes diagonal + neighbors. See `scipy.ndimage.generate_binary_structure`. + data : networkx Graph specification, optional + Initial or additional edges to pass to the NetworkX Graph + constructor. See `networkx.Graph`. Valid edge specifications + include edge list (list of tuples), NumPy arrays, and SciPy + sparse matrices. + **attr : keyword arguments, optional + Additional attributes to add to the graph. """ - def __init__(self, data=None, **attr): + def __init__(self, label_image=None, connectivity=1, data=None, **attr): super(RAG, self).__init__(data, **attr) - try: - self.max_id = max(self.nodes_iter()) - except ValueError: - # Empty sequence + if self.number_of_nodes() == 0: self.max_id = 0 + else: + self.max_id = max(self.nodes_iter()) + + if label_image is not None: + fp = ndi.generate_binary_structure(label_image.ndim, connectivity) + ndi.generic_filter( + label_image, + function=_add_edge_filter, + footprint=fp, + mode='nearest', + output=as_strided(np.empty((1,), dtype=np.float_), + shape=label_image.shape, + strides=((0,) * label_image.ndim)), + extra_arguments=(self,)) def merge_nodes(self, src, dst, weight_func=min_weight, in_place=True, extra_arguments=[], extra_keywords={}): @@ -172,36 +220,6 @@ class RAG(nx.Graph): super(RAG, self).add_node(n) -def _add_edge_filter(values, graph): - """Create edge in `g` between the first element of `values` and the rest. - - Add an edge between the first element in `values` and - all other elements of `values` in the graph `g`. `values[0]` - is expected to be the central value of the footprint used. - - Parameters - ---------- - values : array - The array to process. - graph : RAG - The graph to add edges in. - - Returns - ------- - 0 : int - Always returns 0. The return value is required so that `generic_filter` - can put it in the output array. - - """ - values = values.astype(int) - current = values[0] - for value in values[1:]: - if value != current: - graph.add_edge(current, value) - - return 0 - - def rag_mean_color(image, labels, connectivity=2, mode='distance', sigma=255.0): """Compute the Region Adjacency Graph using mean colors. @@ -224,7 +242,7 @@ def rag_mean_color(image, labels, connectivity=2, mode='distance', Pixels with a squared distance less than `connectivity` from each other are considered adjacent. It can range from 1 to `labels.ndim`. Its behavior is the same as `connectivity` parameter in - `scipy.ndimage.filters.generate_binary_structure`. + `scipy.ndimage.generate_binary_structure`. mode : {'distance', 'similarity'}, optional The strategy to assign edge weights. @@ -263,35 +281,7 @@ def rag_mean_color(image, labels, connectivity=2, mode='distance', http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.11.5274 """ - graph = RAG() - - # The footprint is constructed in such a way that the first - # element in the array being passed to _add_edge_filter is - # the central value. - fp = ndi.generate_binary_structure(labels.ndim, connectivity) - for d in range(fp.ndim): - fp = fp.swapaxes(0, d) - fp[0, ...] = 0 - fp = fp.swapaxes(0, d) - - # For example - # if labels.ndim = 2 and connectivity = 1 - # fp = [[0,0,0], - # [0,1,1], - # [0,1,0]] - # - # if labels.ndim = 2 and connectivity = 2 - # fp = [[0,0,0], - # [0,1,1], - # [0,1,1]] - - ndi.generic_filter( - labels, - function=_add_edge_filter, - footprint=fp, - mode='nearest', - output=np.zeros(labels.shape, dtype=np.uint8), - extra_arguments=(graph,)) + graph = RAG(labels, connectivity=connectivity) for n in graph: graph.node[n].update({'labels': [n], diff --git a/skimage/future/graph/tests/test_rag.py b/skimage/future/graph/tests/test_rag.py index 35c22218..e5882e20 100644 --- a/skimage/future/graph/tests/test_rag.py +++ b/skimage/future/graph/tests/test_rag.py @@ -178,3 +178,21 @@ def test_ncut_stable_subgraph(): new_labels, _, _ = segmentation.relabel_sequential(new_labels) assert new_labels.max() == 0 + + +def test_generic_rag_2d(): + labels = np.array([[1, 2], [3, 4]], dtype=np.uint8) + g = graph.RAG(labels) + assert g.has_edge(1, 2) and g.has_edge(2, 4) and not g.has_edge(1, 4) + h = graph.RAG(labels, connectivity=2) + assert h.has_edge(1, 2) and h.has_edge(1, 4) and h.has_edge(2, 3) + + +def test_generic_rag_3d(): + labels = np.arange(8, dtype=np.uint8).reshape((2, 2, 2)) + g = graph.RAG(labels) + assert g.has_edge(0, 1) and g.has_edge(1, 3) and not g.has_edge(0, 3) + h = graph.RAG(labels, connectivity=2) + assert h.has_edge(0, 1) and h.has_edge(0, 3) and not h.has_edge(0, 7) + k = graph.RAG(labels, connectivity=3) + assert k.has_edge(0, 1) and k.has_edge(1, 2) and k.has_edge(2, 5) diff --git a/skimage/io/__init__.py b/skimage/io/__init__.py index d0ec1f43..21cc5d62 100644 --- a/skimage/io/__init__.py +++ b/skimage/io/__init__.py @@ -58,5 +58,5 @@ def _update_doc(doc): return doc - -__doc__ = _update_doc(__doc__) +if __doc__ is not None: + __doc__ = _update_doc(__doc__) diff --git a/skimage/io/_plugins/pil_plugin.py b/skimage/io/_plugins/pil_plugin.py index 7cb8e5c1..59811dd2 100644 --- a/skimage/io/_plugins/pil_plugin.py +++ b/skimage/io/_plugins/pil_plugin.py @@ -85,7 +85,10 @@ def pil_to_ndarray(im, dtype=None, img_num=None): if grayscale: frame = im.convert('L') else: - frame = im.convert('RGB') + if im.format == 'PNG' and 'transparency' in im.info: + frame = im.convert('RGBA') + else: + frame = im.convert('RGB') elif im.mode == '1': frame = im.convert('L') diff --git a/skimage/io/tests/test_mpl_imshow.py b/skimage/io/tests/test_mpl_imshow.py index 9b814f57..e369273a 100644 --- a/skimage/io/tests/test_mpl_imshow.py +++ b/skimage/io/tests/test_mpl_imshow.py @@ -87,8 +87,7 @@ def test_outside_standard_range(): def test_nonstandard_type(): plt.figure() - with expected_warnings(["Non-standard image type", - "Low image dynamic range"]): + with expected_warnings(["Low image dynamic range"]): ax_im = io.imshow(im64) assert ax_im.get_clim() == (im64.min(), im64.max()) assert n_subplots(ax_im) == 2 diff --git a/skimage/io/tests/test_pil.py b/skimage/io/tests/test_pil.py index ad3870a6..ee962906 100644 --- a/skimage/io/tests/test_pil.py +++ b/skimage/io/tests/test_pil.py @@ -84,6 +84,28 @@ def test_imread_palette(): assert img.ndim == 3 +def test_imread_index_png_with_alpha(): + # The file `foo3x5x4indexed.png` was created with this array + # (3x5 is (height)x(width)): + data = np.array([[[127, 0, 255, 255], + [127, 0, 255, 255], + [127, 0, 255, 255], + [127, 0, 255, 255], + [127, 0, 255, 255]], + [[192, 192, 255, 0], + [192, 192, 255, 0], + [0, 0, 255, 0], + [0, 0, 255, 0], + [0, 0, 255, 0]], + [[0, 31, 255, 255], + [0, 31, 255, 255], + [0, 31, 255, 255], + [0, 31, 255, 255], + [0, 31, 255, 255]]], dtype=np.uint8) + img = imread(os.path.join(data_dir, 'foo3x5x4indexed.png')) + assert_array_equal(img, data) + + def test_palette_is_gray(): gray = Image.open(os.path.join(data_dir, 'palette_gray.png')) assert _palette_is_grayscale(gray) diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py index 96fb1102..cd0e1b60 100644 --- a/skimage/measure/_marching_cubes.py +++ b/skimage/measure/_marching_cubes.py @@ -1,8 +1,10 @@ import numpy as np +import scipy.ndimage as ndi from . import _marching_cubes_cy -def marching_cubes(volume, level, spacing=(1., 1., 1.)): +def marching_cubes(volume, level, spacing=(1., 1., 1.), + gradient_direction='descent'): """ Marching cubes algorithm to find iso-valued surfaces in 3d volumetric data @@ -15,6 +17,12 @@ def marching_cubes(volume, level, spacing=(1., 1., 1.)): spacing : length-3 tuple of floats Voxel spacing in spatial dimensions corresponding to numpy array indexing dimensions (M, N, P) as in `volume`. + gradient_direction : string + Controls if the mesh was generated from an isosurface with gradient + descent toward objects of interest (the default), or the opposite. + The two options are: + * descent : Object was greater than exterior + * ascent : Exterior was greater than object Returns ------- @@ -68,14 +76,10 @@ def marching_cubes(volume, level, spacing=(1., 1., 1.)): lexicographical order) coordinate in the contour. This is a side-effect of how the input array is traversed, but can be relied upon. - The generated mesh does not guarantee coherent orientation because of how - symmetry is used in the algorithm. If this is required, e.g. due to a - particular visualization package or for generating 3D printing STL files, - the utility ``skimage.measure.correct_mesh_orientation`` is available to - fix this in post-processing. + The generated mesh guarantees coherent orientation as of version 0.12. To quantify the area of an isosurface generated by this algorithm, pass - the outputs directly into `skimage.measure.mesh_surface_area`. + outputs directly into `skimage.measure.mesh_surface_area`. Regarding visualization of algorithm output, the ``mayavi`` package is recommended. To contour a volume named `myvolume` about the level 0.0:: @@ -122,8 +126,18 @@ def marching_cubes(volume, level, spacing=(1., 1., 1.)): # Returns a true mesh with no degenerate faces. verts, faces = _marching_cubes_cy.unpack_unique_verts(raw_faces) + verts = np.asarray(verts) + faces = np.asarray(faces) + + # Calculate gradient of `volume`, then interpolate to vertices in `verts` + grad_x, grad_y, grad_z = np.gradient(volume) + + # Fancy indexing to define two vector arrays from triangle vertices + faces = _correct_mesh_orientation(volume, verts[faces], faces, spacing, + gradient_direction) + # Adjust for non-isotropic spacing in `verts` at time of return - return np.asarray(verts) * np.r_[spacing], np.asarray(faces) + return verts * np.r_[spacing], faces def mesh_surface_area(verts, faces): @@ -187,7 +201,7 @@ def correct_mesh_orientation(volume, verts, faces, spacing=(1., 1., 1.), indexing dimensions (M, N, P) as in `volume`. gradient_direction : string Controls if the mesh was generated from an isosurface with gradient - ascent toward objects of interest (the default), or the opposite. + descent toward objects of interest (the default), or the opposite. The two options are: * descent : Object was greater than exterior * ascent : Exterior was greater than object @@ -225,17 +239,86 @@ def correct_mesh_orientation(volume, verts, faces, spacing=(1., 1., 1.), skimage.measure.mesh_surface_area """ - import scipy.ndimage as ndi + import warnings + warnings.warn( + DeprecationWarning("`correct_mesh_orientation` is deprecated for " + "removal as `marching_cubes` now guarantess " + "correct mesh orientation.")) - # Calculate gradient of `volume`, then interpolate to vertices in `verts` - grad_x, grad_y, grad_z = np.gradient(volume) + verts = verts.copy() + verts[:, 0] /= spacing[0] + verts[:, 1] /= spacing[1] + verts[:, 2] /= spacing[2] # 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] - + + return _correct_mesh_orientation(volume, actual_verts, faces, spacing, + gradient_direction) + + +def _correct_mesh_orientation(volume, actual_verts, faces, + spacing=(1., 1., 1.), + gradient_direction='descent'): + """ + Correct orientations of mesh faces. + + Parameters + ---------- + volume : (M, N, P) array of doubles + Input data volume to find isosurfaces. Will be cast to `np.float64`. + actual_verts : (F, 3, 3) array of floats + Array with (face, vertex, coords) index coordinates. + faces : (F, 3) array of ints + List of length-3 lists of integers, referencing vertex coordinates as + provided in `verts`. + spacing : length-3 tuple of floats + Voxel spacing in spatial dimensions corresponding to numpy array + indexing dimensions (M, N, P) as in `volume`. + gradient_direction : string + Controls if the mesh was generated from an isosurface with gradient + descent toward objects of interest (the default), or the opposite. + The two options are: + * descent : Object was greater than exterior + * ascent : Exterior was greater than object + + Returns + ------- + faces_corrected (F, 3) array of ints + Corrected list of faces referencing vertex coordinates in `verts`. + + Notes + ----- + Certain applications and mesh processing algorithms require all faces + to be oriented in a consistent way. Generally, this means a normal vector + points "out" of the meshed shapes. This algorithm corrects the output from + `skimage.measure.marching_cubes` by flipping the orientation of + mis-oriented faces. + + Because marching cubes could be used to find isosurfaces either on + gradient descent (where the desired object has greater values than the + exterior) or ascent (where the desired object has lower values than the + exterior), the ``gradient_direction`` kwarg allows the user to inform this + algorithm which is correct. If the resulting mesh appears to be oriented + completely incorrectly, try changing this option. + + The arguments expected by this function are the exact outputs from + `skimage.measure.marching_cubes` except `actual_verts`, which is an + uncorrected version of the fancy indexing operation `verts[faces]`. + Only `faces` is corrected and returned as the vertices do not change, + only the order in which they are referenced. + + This algorithm assumes ``faces`` provided are exclusively triangles. + + See Also + -------- + skimage.measure.marching_cubes + skimage.measure.mesh_surface_area + + """ + # Calculate gradient of `volume`, then interpolate to vertices in `verts` + grad_x, grad_y, grad_z = np.gradient(volume) + a = actual_verts[:, 0, :] - actual_verts[:, 1, :] b = actual_verts[:, 0, :] - actual_verts[:, 2, :] diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index f38fda22..297a6390 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -1,4 +1,5 @@ # coding: utf-8 +from __future__ import division from math import sqrt, atan2, pi as PI import numpy as np from scipy import ndimage as ndi @@ -7,6 +8,9 @@ from ._label import label from . import _moments +from functools import wraps +from collections import defaultdict + __all__ = ['regionprops', 'perimeter'] @@ -14,20 +18,21 @@ STREL_4 = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]], dtype=np.uint8) STREL_8 = np.ones((3, 3), dtype=np.uint8) +STREL_26_3D = np.ones((3, 3, 3), dtype=np.uint8) PROPS = { 'Area': 'area', 'BoundingBox': 'bbox', 'CentralMoments': 'moments_central', 'Centroid': 'centroid', 'ConvexArea': 'convex_area', -# 'ConvexHull', + # 'ConvexHull', 'ConvexImage': 'convex_image', 'Coordinates': 'coords', 'Eccentricity': 'eccentricity', 'EquivDiameter': 'equivalent_diameter', 'EulerNumber': 'euler_number', 'Extent': 'extent', -# 'Extrema', + # 'Extrema', 'FilledArea': 'filled_area', 'FilledImage': 'filled_image', 'HuMoments': 'moments_hu', @@ -42,10 +47,10 @@ PROPS = { 'NormalizedMoments': 'moments_normalized', 'Orientation': 'orientation', 'Perimeter': 'perimeter', -# 'PixelIdxList', -# 'PixelList', + # 'PixelIdxList', + # 'PixelList', 'Solidity': 'solidity', -# 'SubarrayIdx' + # 'SubarrayIdx' 'WeightedCentralMoments': 'weighted_moments_central', 'WeightedCentroid': 'weighted_centroid', 'WeightedHuMoments': 'weighted_moments_hu', @@ -53,131 +58,121 @@ PROPS = { 'WeightedNormalizedMoments': 'weighted_moments_normalized' } -PROP_VALS = PROPS.values() +PROP_VALS = set(PROPS.values()) -class _cached_property(object): - """Decorator to use a function as a cached property. +def _cached(f): + @wraps(f) + def wrapper(obj): + cache = obj._cache + prop = f.__name__ - The function is only called the first time and each successive call returns - the cached result of the first call. + if not ((prop in cache) and obj._cache_active): + cache[prop] = f(obj) - class Foo(object): + return cache[prop] - @_cached_property - def foo(self): - return "Cached" + return wrapper - class Foo(object): - def __init__(self): - self._cache_active = False - - @_cached_property - def foo(self): - return "Not cached" - - Adapted from . - - """ - - def __init__(self, func, name=None, doc=None): - self.__name__ = name or func.__name__ - self.__module__ = func.__module__ - self.__doc__ = doc or func.__doc__ - self.func = func - - def __get__(self, obj, type=None): - if obj is None: - return self - - # call every time, if cache is not active - if not obj.__dict__.get('_cache_active', True): - return self.func(obj) - - # try to retrieve from cache or call and store result in cache - try: - value = obj.__dict__[self.__name__] - except KeyError: - value = self.func(obj) - obj.__dict__[self.__name__] = value - return value +def only2d(method): + @wraps(method) + def func2d(self, *args, **kwargs): + if self._ndim > 2: + raise NotImplementedError('Property %s is not implemented for ' + '3D images' % method.__name__) + return method(self, *args, **kwargs) + return func2d class _RegionProperties(object): + """Please refer to `skimage.measure.regionprops` for more information + on the available region properties. + """ def __init__(self, slice, label, label_image, intensity_image, cache_active): + + if intensity_image is not None: + if not intensity_image.shape == label_image.shape: + raise ValueError('Label and intensity image must have the same shape.') + self.label = label + self._slice = slice self._label_image = label_image self._intensity_image = intensity_image + self._cache_active = cache_active + self._cache = {} + self._ndim = label_image.ndim - @property + @_cached def area(self): - return self.moments[0, 0] + return np.sum(self.image) - @property def bbox(self): - return (self._slice[0].start, self._slice[1].start, - self._slice[0].stop, self._slice[1].stop) + return tuple([self._slice[i].start for i in range(self._ndim)] + + [self._slice[i].stop for i in range(self._ndim)]) - @property def centroid(self): - row, col = self.local_centroid - return row + self._slice[0].start, col + self._slice[1].start + centroid_coords = self.local_centroid + return tuple(centroid_coords + + np.array([self._slice[i].start + for i in range(self._ndim)])) - @property + @only2d def convex_area(self): return np.sum(self.convex_image) - @_cached_property + @_cached + @only2d def convex_image(self): from ..morphology.convex_hull import convex_hull_image return convex_hull_image(self.image) - @property def coords(self): - rr, cc = np.nonzero(self.image) - return np.vstack((rr + self._slice[0].start, - cc + self._slice[1].start)).T + indices = np.nonzero(self.image) + return np.vstack([indices[i] + self._slice[i].start + for i in range(self._ndim)]).T - @property + @only2d def eccentricity(self): l1, l2 = self.inertia_tensor_eigvals if l1 == 0: return 0 return sqrt(1 - l2 / l1) - @property def equivalent_diameter(self): - return sqrt(4 * self.moments[0, 0] / PI) + if self._ndim == 2: + return sqrt(4 * self.area / PI) + elif self._ndim == 3: + return (6 * self.area / PI) ** (1. / 3) - @property + @only2d def euler_number(self): euler_array = self.filled_image != self.image - _, num = label(euler_array, neighbors=8, return_num=True, background=-1) + _, num = label(euler_array, neighbors=8, return_num=True, + background=0) return -num + 1 - @property def extent(self): - rows, cols = self.image.shape - return self.moments[0, 0] / (rows * cols) + return self.area / self.image.size - @property def filled_area(self): return np.sum(self.filled_image) - @_cached_property + @_cached def filled_image(self): - return ndi.binary_fill_holes(self.image, STREL_8) + structure = STREL_8 if self._ndim == 2 else STREL_26_3D + return ndi.binary_fill_holes(self.image, structure) - @_cached_property + @_cached def image(self): return self._label_image[self._slice] == self.label - @_cached_property + @_cached + @only2d def inertia_tensor(self): mu = self.moments_central a = mu[2, 0] / mu[0, 0] @@ -185,7 +180,8 @@ class _RegionProperties(object): c = mu[0, 2] / mu[0, 0] return np.array([[a, b], [b, c]]) - @_cached_property + @_cached + @only2d def inertia_tensor_eigvals(self): a, b, b, c = self.inertia_tensor.flat # eigen values of inertia tensor @@ -193,64 +189,63 @@ class _RegionProperties(object): l2 = (a + c) / 2 - sqrt(4 * b ** 2 + (a - c) ** 2) / 2 return l1, l2 - @_cached_property + @_cached def intensity_image(self): if self._intensity_image is None: raise AttributeError('No intensity image specified.') return self._intensity_image[self._slice] * self.image - @property def _intensity_image_double(self): return self.intensity_image.astype(np.double) - @property + @only2d def local_centroid(self): m = self.moments row = m[0, 1] / m[0, 0] col = m[1, 0] / m[0, 0] return row, col - @property def max_intensity(self): return np.max(self.intensity_image[self.image]) - @property def mean_intensity(self): return np.mean(self.intensity_image[self.image]) - @property def min_intensity(self): return np.min(self.intensity_image[self.image]) - @property + @only2d def major_axis_length(self): l1, _ = self.inertia_tensor_eigvals return 4 * sqrt(l1) - @property + @only2d def minor_axis_length(self): _, l2 = self.inertia_tensor_eigvals return 4 * sqrt(l2) - @_cached_property + @_cached + @only2d def moments(self): return _moments.moments(self.image.astype(np.uint8), 3) - @_cached_property + @_cached + @only2d def moments_central(self): row, col = self.local_centroid return _moments.moments_central(self.image.astype(np.uint8), row, col, 3) - @property + @only2d def moments_hu(self): return _moments.moments_hu(self.moments_normalized) - @_cached_property + @_cached + @only2d def moments_normalized(self): return _moments.moments_normalized(self.moments_central, 3) - @property + @only2d def orientation(self): a, b, b, c = self.inertia_tensor.flat b = -b @@ -262,46 +257,65 @@ class _RegionProperties(object): else: return - 0.5 * atan2(2 * b, (a - c)) - @property + @only2d def perimeter(self): return perimeter(self.image, 4) - @property + @only2d def solidity(self): return self.moments[0, 0] / np.sum(self.convex_image) - @property + @only2d def weighted_centroid(self): row, col = self.weighted_local_centroid return row + self._slice[0].start, col + self._slice[1].start - @property + @only2d def weighted_local_centroid(self): m = self.weighted_moments row = m[0, 1] / m[0, 0] col = m[1, 0] / m[0, 0] return row, col - @_cached_property + @_cached + @only2d def weighted_moments(self): - return _moments.moments_central(self._intensity_image_double, 0, 0, 3) + return _moments.moments_central(self._intensity_image_double(), 0, 0, 3) - @_cached_property + @_cached + @only2d def weighted_moments_central(self): row, col = self.weighted_local_centroid - return _moments.moments_central(self._intensity_image_double, + return _moments.moments_central(self._intensity_image_double(), row, col, 3) - @property + @only2d def weighted_moments_hu(self): return _moments.moments_hu(self.weighted_moments_normalized) - @_cached_property + @_cached + @only2d def weighted_moments_normalized(self): return _moments.moments_normalized(self.weighted_moments_central, 3) def __iter__(self): - return iter(PROPS.values()) + props = PROP_VALS + + if self._intensity_image is None: + unavailable_props = ('intensity_image', + 'max_intensity', + 'mean_intensity', + 'min_intensity', + 'weighted_moments', + 'weighted_moments_central', + 'weighted_centroid', + 'weighted_local_centroid', + 'weighted_moments_hu', + 'weighted_moments_normalized') + + props = props.difference(unavailable_props) + + return iter(sorted(props)) def __getitem__(self, key): value = getattr(self, key, None) @@ -316,7 +330,7 @@ class _RegionProperties(object): for key in PROP_VALS: try: - #so that NaNs are equal + # so that NaNs are equal np.testing.assert_equal(getattr(self, key, None), getattr(other, key, None)) except AssertionError: @@ -333,7 +347,8 @@ def regionprops(label_image, intensity_image=None, cache=True): label_image : (N, M) ndarray Labeled input image. Labels with value 0 are ignored. intensity_image : (N, M) ndarray, optional - Intensity image with same size as labeled image. Default is None. + Intensity (i.e., input) image with same size as labeled image. + Default is None. cache : bool, optional Determine whether to cache calculated properties. The computation is much faster for cached properties, whereas the memory consumption @@ -352,7 +367,7 @@ def regionprops(label_image, intensity_image=None, cache=True): **area** : int Number of pixels of region. **bbox** : tuple - Bounding box ``(min_row, min_col, max_row, max_col)`` + Bounding box ``(min_row, min_col, max_row, max_col)`` **centroid** : array Centroid coordinate tuple ``(row, col)``. **convex_area** : int @@ -370,7 +385,7 @@ def regionprops(label_image, intensity_image=None, cache=True): **equivalent_diameter** : float The diameter of a circle with the same area as the region. **euler_number** : int - Euler number of region. Computed as number of objects (= 1) + Euler characteristic of region. Computed as number of objects (= 1) subtracted by number of holes (8-connectivity). **extent** : float Ratio of pixels in the region to pixels in the total bounding box. @@ -386,8 +401,13 @@ def regionprops(label_image, intensity_image=None, cache=True): Inertia tensor of the region for the rotation around its mass. **inertia_tensor_eigvals** : tuple The two eigen values of the inertia tensor in decreasing order. + **intensity_image** : ndarray + Image inside region bounding box. **label** : int The label in the labeled input image. + **local_centroid** : array + Centroid coordinate tuple ``(row, col)``, relative to region bounding + box. **major_axis_length** : float The length of the major axis of the ellipse that has the same normalized second central moments as the region. @@ -433,6 +453,9 @@ def regionprops(label_image, intensity_image=None, cache=True): **weighted_centroid** : array Centroid coordinate tuple ``(row, col)`` weighted with intensity image. + **weighted_local_centroid** : array + Centroid coordinate tuple ``(row, col)``, relative to region bounding + box, weighted with intensity image. **weighted_moments** : (3, 3) ndarray Spatial moments of intensity image up to 3rd order:: @@ -457,7 +480,12 @@ def regionprops(label_image, intensity_image=None, cache=True): wnu_ji = wmu_ji / wm_00^[(i+j)/2 + 1] - where `wm_00` is the zeroth spatial moment (intensity-weighted area). + where ``wm_00`` is the zeroth spatial moment (intensity-weighted area). + + Each region also supports iteration, so that you can do:: + + for prop in region: + print(prop, region[prop]) References ---------- @@ -488,8 +516,8 @@ def regionprops(label_image, intensity_image=None, cache=True): label_image = np.squeeze(label_image) - if label_image.ndim != 2: - raise TypeError('Only 2-D images supported.') + if label_image.ndim not in (2, 3): + raise TypeError('Only 2-D and 3-D images supported.') if not np.issubdtype(label_image.dtype, np.integer): raise TypeError('Label image must be of integral type.') @@ -544,7 +572,6 @@ def perimeter(image, neighbourhood=4): perimeter_weights[[21, 33]] = sqrt(2) perimeter_weights[[13, 23]] = (1 + sqrt(2)) / 2 - perimeter_image = ndi.convolve(border_image, np.array([[10, 2, 10], [ 2, 1, 2], [10, 2, 10]]), @@ -557,3 +584,32 @@ def perimeter(image, neighbourhood=4): perimeter_histogram = np.bincount(perimeter_image.ravel(), minlength=50) total_perimeter = np.dot(perimeter_histogram, perimeter_weights) return total_perimeter + + +def _parse_docs(): + import re + import textwrap + + doc = regionprops.__doc__ + matches = re.finditer('\*\*(\w+)\*\* \:.*?\n(.*?)(?=\n [\*\S]+)', + doc, flags=re.DOTALL) + prop_doc = dict((m.group(1), textwrap.dedent(m.group(2))) for m in matches) + + return prop_doc + + +def _install_properties_docs(): + prop_doc = _parse_docs() + + for p in [member for member in dir(_RegionProperties) + if not member.startswith('_')]: + try: + getattr(_RegionProperties, p).__doc__ = prop_doc[p] + except AttributeError: + # For Python 2.x + getattr(_RegionProperties, p).im_func.__doc__ = prop_doc[p] + + setattr(_RegionProperties, p, property(getattr(_RegionProperties, p))) + + +_install_properties_docs() diff --git a/skimage/measure/tests/test_marching_cubes.py b/skimage/measure/tests/test_marching_cubes.py index 60adbb0d..0e3d77ec 100644 --- a/skimage/measure/tests/test_marching_cubes.py +++ b/skimage/measure/tests/test_marching_cubes.py @@ -39,7 +39,24 @@ def test_invalid_input(): def test_correct_mesh_orientation(): sphere_small = ellipsoid(1, 1, 1, levelset=True) - verts, faces = marching_cubes(sphere_small, 0.) + + # Mesh with incorrectly oriented faces which was previously returned from + # `marching_cubes`, before it guaranteed correct mesh orientation + verts = np.array([[1., 2., 2.], + [2., 2., 1.], + [2., 1., 2.], + [2., 2., 3.], + [2., 3., 2.], + [3., 2., 2.]]) + + faces = np.array([[0, 1, 2], + [2, 0, 3], + [1, 0, 4], + [4, 0, 3], + [1, 2, 5], + [2, 3, 5], + [1, 4, 5], + [5, 4, 3]]) # Correct mesh orientation - descent corrected_faces1 = correct_mesh_orientation(sphere_small, verts, faces, diff --git a/skimage/measure/tests/test_regionprops.py b/skimage/measure/tests/test_regionprops.py index 6865818e..f724577d 100644 --- a/skimage/measure/tests/test_regionprops.py +++ b/skimage/measure/tests/test_regionprops.py @@ -3,8 +3,8 @@ from numpy.testing import assert_array_equal, assert_almost_equal, \ import numpy as np import math -from skimage.measure._regionprops import regionprops, PROPS, perimeter -from skimage._shared._warnings import expected_warnings +from skimage.measure._regionprops import (regionprops, PROPS, perimeter, + _parse_docs) SAMPLE = np.array( @@ -22,6 +22,10 @@ SAMPLE = np.array( INTENSITY_SAMPLE = SAMPLE.copy() INTENSITY_SAMPLE[1, 9:11] = 2 +SAMPLE_3D = np.zeros((6, 6, 6), dtype=np.uint8) +SAMPLE_3D[1:3, 1:3, 1:3] = 1 +SAMPLE_3D[3, 2, 2] = 1 +INTENSITY_SAMPLE_3D = SAMPLE_3D.copy() def test_all_props(): region = regionprops(SAMPLE, INTENSITY_SAMPLE)[0] @@ -29,6 +33,14 @@ def test_all_props(): assert_almost_equal(region[prop], getattr(region, PROPS[prop])) +def test_all_props_3d(): + region = regionprops(SAMPLE_3D, INTENSITY_SAMPLE_3D)[0] + for prop in PROPS: + try: + assert_almost_equal(region[prop], getattr(region, PROPS[prop])) + except NotImplementedError: + pass + def test_dtype(): regionprops(np.zeros((10, 10), dtype=np.int)) regionprops(np.zeros((10, 10), dtype=np.uint)) @@ -42,12 +54,15 @@ def test_ndim(): regionprops(np.zeros((10, 10), dtype=np.int)) regionprops(np.zeros((10, 10, 1), dtype=np.int)) regionprops(np.zeros((10, 10, 1, 1), dtype=np.int)) - assert_raises(TypeError, regionprops, np.zeros((10, 10, 2), dtype=np.int)) + regionprops(np.zeros((10, 10, 10), dtype=np.int)) + assert_raises(TypeError, regionprops, np.zeros((10, 10, 10, 2), dtype=np.int)) def test_area(): area = regionprops(SAMPLE)[0].area assert area == np.sum(SAMPLE) + area = regionprops(SAMPLE_3D)[0].area + assert area == np.sum(SAMPLE_3D) def test_bbox(): @@ -59,6 +74,9 @@ def test_bbox(): bbox = regionprops(SAMPLE_mod)[0].bbox assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1]-1)) + bbox = regionprops(SAMPLE_3D)[0].bbox + assert_array_almost_equal(bbox, (1, 1, 1, 4, 3, 3)) + def test_moments_central(): mu = regionprops(SAMPLE)[0].moments_central @@ -110,6 +128,11 @@ def test_coordinates(): prop_coords = regionprops(sample)[0].coords assert_array_equal(prop_coords, coords) + sample = np.zeros((6, 6, 6), dtype=np.int8) + coords = np.array([[1, 1, 1], [1, 2, 1], [1, 3, 1]]) + sample[coords[:, 0], coords[:, 1], coords[:, 2]] = 1 + prop_coords = regionprops(sample)[0].coords + assert_array_equal(prop_coords, coords) def test_eccentricity(): eps = regionprops(SAMPLE)[0].eccentricity @@ -129,12 +152,12 @@ def test_equiv_diameter(): def test_euler_number(): en = regionprops(SAMPLE)[0].euler_number - assert en == 0 + assert en == 1 SAMPLE_mod = SAMPLE.copy() SAMPLE_mod[7, -3] = 0 en = regionprops(SAMPLE_mod)[0].euler_number - assert en == -1 + assert en == 0 def test_extent(): @@ -145,13 +168,13 @@ def test_extent(): def test_moments_hu(): hu = regionprops(SAMPLE)[0].moments_hu ref = np.array([ - 3.27117627e-01, - 2.63869194e-02, - 2.35390060e-02, - 1.23151193e-03, - 1.38882330e-06, + 3.27117627e-01, + 2.63869194e-02, + 2.35390060e-02, + 1.23151193e-03, + 1.38882330e-06, -2.72586158e-05, - 6.48350653e-06 + 6.48350653e-06 ]) # bug in OpenCV caused in Central Moments calculation? assert_array_almost_equal(hu, ref) @@ -161,11 +184,17 @@ def test_image(): img = regionprops(SAMPLE)[0].image assert_array_equal(img, SAMPLE) + img = regionprops(SAMPLE_3D)[0].image + assert_array_equal(img, SAMPLE_3D[1:4, 1:3, 1:3]) + def test_label(): label = regionprops(SAMPLE)[0].label assert_array_equal(label, 1) + label = regionprops(SAMPLE_3D)[0].label + assert_array_equal(label, 1) + def test_filled_area(): area = regionprops(SAMPLE)[0].filled_area @@ -355,11 +384,18 @@ def test_pure_background(): def test_invalid(): ps = regionprops(SAMPLE) + def get_intensity_image(): ps[0].intensity_image + assert_raises(AttributeError, get_intensity_image) +def test_invalid_size(): + wrong_intensity_sample = np.array([[1], [1]]) + assert_raises(ValueError, regionprops, SAMPLE, wrong_intensity_sample) + + def test_equals(): arr = np.zeros((100, 100), dtype=np.int) arr[0:25, 0:25] = 1 @@ -376,6 +412,47 @@ def test_equals(): assert_equal(r1 != r3, True, "Different regionprops are equal") +def test_iterate_all_props(): + region = regionprops(SAMPLE)[0] + p0 = dict((p, region[p]) for p in region) + + region = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[0] + p1 = dict((p, region[p]) for p in region) + + assert len(p0) < len(p1) + + +def test_cache(): + region = regionprops(SAMPLE)[0] + f0 = region.filled_image + region._label_image[:10] = 1 + f1 = region.filled_image + + # Changed underlying image, but cache keeps result the same + assert_array_equal(f0, f1) + + # Now invalidate cache + region._cache_active = False + f1 = region.filled_image + + assert np.any(f0 != f1) + + +def test_docstrings_and_props(): + region = regionprops(SAMPLE)[0] + + docs = _parse_docs() + props = [m for m in dir(region) if not m.startswith('_')] + + nr_docs_parsed = len(docs) + nr_props = len(props) + assert_equal(nr_docs_parsed, nr_props) + + ds = docs['weighted_moments_normalized'] + assert 'iteration' not in ds + assert len(ds.split('\n')) > 3 + + if __name__ == "__main__": from numpy.testing import run_module_suite run_module_suite() diff --git a/skimage/morphology/__init__.py b/skimage/morphology/__init__.py index 159a4ac2..a313f4ca 100644 --- a/skimage/morphology/__init__.py +++ b/skimage/morphology/__init__.py @@ -11,8 +11,6 @@ from .greyreconstruct import reconstruction from .misc import remove_small_objects, remove_small_holes from ..measure._label import label -from .._shared.utils import deprecated as _deprecated -label = _deprecated('skimage.measure.label')(label) __all__ = ['binary_erosion', diff --git a/skimage/morphology/convex_hull.py b/skimage/morphology/convex_hull.py index bb9ca734..31b92a47 100644 --- a/skimage/morphology/convex_hull.py +++ b/skimage/morphology/convex_hull.py @@ -1,11 +1,12 @@ -__all__ = ['convex_hull_image', 'convex_hull_object'] - +"""Convex Hull.""" import numpy as np from ..measure._pnpoly import grid_points_in_poly from ._convex_hull import possible_hull from ..measure._label import label from ..util import unique_rows +__all__ = ['convex_hull_image', 'convex_hull_object'] + try: from scipy.spatial import Delaunay except ImportError: @@ -33,6 +34,8 @@ def convex_hull_image(image): .. [1] http://blogs.mathworks.com/steve/2011/10/04/binary-image-convex-hull-algorithm-notes/ """ + if image.ndim > 2: + raise ValueError("Input must be a 2D image") if Delaunay is None: raise ImportError("Could not import scipy.spatial.Delaunay, " @@ -85,7 +88,7 @@ def convex_hull_object(image, neighbors=8): Parameters ---------- - image : ndarray + image : (M, N) array Binary input image. neighbors : {4, 8}, int Whether to use 4- or 8-connectivity. @@ -104,6 +107,8 @@ def convex_hull_object(image, neighbors=8): convex_hull_image separately on each object. """ + if image.ndim > 2: + raise ValueError("Input must be a 2D image") if neighbors != 4 and neighbors != 8: raise ValueError('Neighbors must be either 4 or 8.') diff --git a/skimage/morphology/tests/test_convex_hull.py b/skimage/morphology/tests/test_convex_hull.py index c0cc954e..1a41b6c8 100644 --- a/skimage/morphology/tests/test_convex_hull.py +++ b/skimage/morphology/tests/test_convex_hull.py @@ -31,20 +31,24 @@ def test_basic(): assert_array_equal(convex_hull_image(image), expected) + # Test that an error is raised on passing a 3D image: + image3d = np.empty((5, 5, 5)) + assert_raises(ValueError, convex_hull_image, image3d) + @skipif(not scipy_spatial) def test_qhull_offset_example(): - nonzeros = (([1367, 1368, 1368, 1368, 1369, 1369, 1369, 1369, 1369, 1370, 1370, - 1370, 1370, 1370, 1370, 1370, 1371, 1371, 1371, 1371, 1371, 1371, - 1371, 1371, 1371, 1372, 1372, 1372, 1372, 1372, 1372, 1372, 1372, - 1372, 1373, 1373, 1373, 1373, 1373, 1373, 1373, 1373, 1373, 1374, - 1374, 1374, 1374, 1374, 1374, 1374, 1375, 1375, 1375, 1375, 1375, - 1376, 1376, 1376, 1377]), - ([151, 150, 151, 152, 149, 150, 151, 152, 153, 148, 149, 150, 151, - 152, 153, 154, 147, 148, 149, 150, 151, 152, 153, 154, 155, 146, - 147, 148, 149, 150, 151, 152, 153, 154, 146, 147, 148, 149, 150, - 151, 152, 153, 154, 147, 148, 149, 150, 151, 152, 153, 148, 149, - 150, 151, 152, 149, 150, 151, 150])) + nonzeros = (([1367, 1368, 1368, 1368, 1369, 1369, 1369, 1369, 1369, 1370, + 1370, 1370, 1370, 1370, 1370, 1370, 1371, 1371, 1371, 1371, + 1371, 1371, 1371, 1371, 1371, 1372, 1372, 1372, 1372, 1372, + 1372, 1372, 1372, 1372, 1373, 1373, 1373, 1373, 1373, 1373, + 1373, 1373, 1373, 1374, 1374, 1374, 1374, 1374, 1374, 1374, + 1375, 1375, 1375, 1375, 1375, 1376, 1376, 1376, 1377]), + ([151, 150, 151, 152, 149, 150, 151, 152, 153, 148, 149, 150, + 151, 152, 153, 154, 147, 148, 149, 150, 151, 152, 153, 154, + 155, 146, 147, 148, 149, 150, 151, 152, 153, 154, 146, 147, + 148, 149, 150, 151, 152, 153, 154, 147, 148, 149, 150, 151, + 152, 153, 148, 149, 150, 151, 152, 149, 150, 151, 150])) image = np.zeros((1392, 1040), dtype=bool) image[nonzeros] = True expected = image.copy() @@ -139,6 +143,9 @@ def test_object(): assert_raises(ValueError, convex_hull_object, image, 7) + # Test that an error is raised on passing a 3D image: + image3d = np.empty((5, 5, 5)) + assert_raises(ValueError, convex_hull_object, image3d) if __name__ == "__main__": np.testing.run_module_suite() diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index 16fd07df..97fab246 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -22,8 +22,8 @@ def denoise_bilateral(image, win_size=5, sigma_range=None, sigma_spatial=1, Parameters ---------- - image : ndarray - Input image. + image : ndarray, shape (M, N[, 3]) + Input image, 2D grayscale or RGB. win_size : int Window size for filtering. sigma_range : float @@ -31,7 +31,7 @@ def denoise_bilateral(image, win_size=5, sigma_range=None, sigma_spatial=1, similarity). A larger value results in averaging of pixels with larger radiometric differences. Note, that the image will be converted using the `img_as_float` function and thus the standard deviation is in - respect to the range ``[0, 1]``. If the value is ``None`` the standard + respect to the range ``[0, 1]``. If the value is ``None`` the standard deviation of the ``image`` will be used. sigma_spatial : float Standard deviation for range distance. A larger value results in @@ -116,7 +116,7 @@ def denoise_tv_bregman(image, weight, max_iter=100, eps=1e-3, isotropic=True): return _denoise_tv_bregman(image, weight, max_iter, eps, isotropic) -def _denoise_tv_chambolle_3d(im, weight=100, eps=2.e-4, n_iter_max=200): +def _denoise_tv_chambolle_3d(im, weight=0.2, eps=2.e-4, n_iter_max=200): """Perform total-variation denoising on 3D images. Parameters @@ -189,7 +189,7 @@ def _denoise_tv_chambolle_3d(im, weight=100, eps=2.e-4, n_iter_max=200): return out -def _denoise_tv_chambolle_2d(im, weight=50, eps=2.e-4, n_iter_max=200): +def _denoise_tv_chambolle_2d(im, weight=0.2, eps=2.e-4, n_iter_max=200): """Perform total-variation denoising on 2D images. Parameters @@ -265,7 +265,7 @@ def _denoise_tv_chambolle_2d(im, weight=50, eps=2.e-4, n_iter_max=200): return out -def denoise_tv_chambolle(im, weight=50, eps=2.e-4, n_iter_max=200, +def denoise_tv_chambolle(im, weight=0.2, eps=2.e-4, n_iter_max=200, multichannel=False): """Perform total-variation denoising on 2D and 3D images. diff --git a/skimage/restoration/_denoise_cy.pyx b/skimage/restoration/_denoise_cy.pyx index b679b488..3c71b2e4 100644 --- a/skimage/restoration/_denoise_cy.pyx +++ b/skimage/restoration/_denoise_cy.pyx @@ -6,7 +6,6 @@ cimport numpy as cnp import numpy as np from libc.math cimport exp, fabs, sqrt -from libc.stdlib cimport malloc, free from libc.float cimport DBL_MAX from .._shared.interpolation cimport get_pixel3d from ..util import img_as_float @@ -16,10 +15,10 @@ cdef inline double _gaussian_weight(double sigma, double value): return exp(-0.5 * (value / sigma)**2) -cdef double* _compute_color_lut(Py_ssize_t bins, double sigma, double max_value): +cdef double[:] _compute_color_lut(Py_ssize_t bins, double sigma, double max_value): cdef: - double* color_lut = malloc(bins * sizeof(double)) + double[:] color_lut = np.empty(bins, dtype=np.double) Py_ssize_t b for b in range(bins): @@ -28,10 +27,10 @@ cdef double* _compute_color_lut(Py_ssize_t bins, double sigma, double max_value) return color_lut -cdef double* _compute_range_lut(Py_ssize_t win_size, double sigma): +cdef double[:] _compute_range_lut(Py_ssize_t win_size, double sigma): cdef: - double* range_lut = malloc(win_size**2 * sizeof(double)) + double[:] range_lut = np.empty(win_size**2, dtype=np.double) Py_ssize_t kr, kc Py_ssize_t window_ext = (win_size - 1) / 2 double dist @@ -74,16 +73,16 @@ def _denoise_bilateral(image, Py_ssize_t win_size, sigma_range, double[:, :, ::1] cimage double[:, :, ::1] out - double* color_lut - double* range_lut + double[:] color_lut + double[:] range_lut Py_ssize_t r, c, d, wr, wc, kr, kc, rr, cc, pixel_addr, color_lut_bin double value, weight, dist, total_weight, csigma_range, color_weight, \ range_weight double dist_scale - double* values - double* centres - double* total_values + double[:] values + double[:] centres + double[:] total_values if sigma_range is None: csigma_range = image.std() @@ -95,20 +94,20 @@ def _denoise_bilateral(image, Py_ssize_t win_size, sigma_range, if max_value == 0.0: raise ValueError("The maximum value found in the image was 0.") + if mode not in ('constant', 'wrap', 'symmetric', 'reflect', 'edge'): + raise ValueError("Invalid mode specified. Please use `constant`, " + "`edge`, `wrap`, `symmetric` or `reflect`.") + cdef char cmode = ord(mode[0].upper()) + cimage = np.ascontiguousarray(image) out = np.zeros((rows, cols, dims), dtype=np.double) color_lut = _compute_color_lut(bins, csigma_range, max_value) range_lut = _compute_range_lut(win_size, sigma_spatial) dist_scale = bins / dims / max_value - values = malloc(dims * sizeof(double)) - centres = malloc(dims * sizeof(double)) - total_values = malloc(dims * sizeof(double)) - - if mode not in ('constant', 'wrap', 'symmetric', 'reflect', 'edge'): - raise ValueError("Invalid mode specified. Please use `constant`, " - "`edge`, `wrap`, `symmetric` or `reflect`.") - cdef char cmode = ord(mode[0].upper()) + values = np.empty(dims, dtype=np.double) + centres = np.empty(dims, dtype=np.double) + total_values = np.empty(dims, dtype=np.double) for r in range(rows): for c in range(cols): @@ -146,12 +145,6 @@ def _denoise_bilateral(image, Py_ssize_t win_size, sigma_range, for d in range(dims): out[r, c, d] = total_values[d] / total_weight - free(color_lut) - free(range_lut) - free(values) - free(centres) - free(total_values) - return np.squeeze(np.asarray(out)) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index 089914fc..7d7ef45b 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -359,6 +359,11 @@ cdef inline float _integral_to_distance_2d(IMGDTYPE [:, ::] integral, """ References ---------- + J. Darbon, A. Cunha, T.F. Chan, S. Osher, and G.J. Jensen, Fast + nonlocal filtering applied to electron cryomicroscopy, in 5th IEEE + International Symposium on Biomedical Imaging: From Nano to Macro, + 2008, pp. 1331-1334. + Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326. @@ -381,6 +386,11 @@ cdef inline float _integral_to_distance_3d(IMGDTYPE [:, :, ::] integral, """ References ---------- + J. Darbon, A. Cunha, T.F. Chan, S. Osher, and G.J. Jensen, Fast + nonlocal filtering applied to electron cryomicroscopy, in 5th IEEE + International Symposium on Biomedical Imaging: From Nano to Macro, + 2008, pp. 1331-1334. + Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326. @@ -523,6 +533,16 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): ------- result : ndarray Denoised image, of same shape as input image. + + References + ---------- + J. Darbon, A. Cunha, T.F. Chan, S. Osher, and G.J. Jensen, Fast + nonlocal filtering applied to electron cryomicroscopy, in 5th IEEE + International Symposium on Biomedical Imaging: From Nano to Macro, + 2008, pp. 1331-1334. + + Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means + Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326. """ if s % 2 == 0: s += 1 # odd value for symmetric patch @@ -623,6 +643,16 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): ------- result : ndarray Denoised image, of same shape as input image. + + References + ---------- + J. Darbon, A. Cunha, T.F. Chan, S. Osher, and G.J. Jensen, Fast + nonlocal filtering applied to electron cryomicroscopy, in 5th IEEE + International Symposium on Biomedical Imaging: From Nano to Macro, + 2008, pp. 1331-1334. + + Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means + Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326. """ if s % 2 == 0: s += 1 # odd value for symmetric patch diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py index 632d6788..35629e9b 100644 --- a/skimage/restoration/deconvolution.py +++ b/skimage/restoration/deconvolution.py @@ -7,7 +7,7 @@ from __future__ import division import numpy as np import numpy.random as npr -from scipy.signal import convolve2d +from scipy.signal import fftconvolve, convolve from . import uft @@ -336,7 +336,7 @@ def richardson_lucy(image, psf, iterations=50, clip=True): Parameters ---------- image : ndarray - Input degraded image. + Input degraded image (can be N dimensional). psf : ndarray The point spread function. iterations : int @@ -365,13 +365,29 @@ def richardson_lucy(image, psf, iterations=50, clip=True): ---------- .. [1] http://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution """ + # compute the times for direct convolution and the fft method. The fft is of + # complexity O(N log(N)) for each dimension and the direct method does + # straight arithmetic (and is O(n*k) to add n elements k times) + direct_time = np.prod(image.shape + psf.shape) + fft_time = np.sum([n*np.log(n) for n in image.shape + psf.shape]) + + # see whether the fourier transform convolution method or the direct + # convolution method is faster (discussed in scikit-image PR #1792) + time_ratio = 40.032 * fft_time / direct_time + + if time_ratio <= 1 or len(image.shape) > 2: + convolve_method = fftconvolve + else: + convolve_method = convolve + image = image.astype(np.float) psf = psf.astype(np.float) im_deconv = 0.5 * np.ones(image.shape) psf_mirror = psf[::-1, ::-1] + for _ in range(iterations): - relative_blur = image / convolve2d(im_deconv, psf, 'same') - im_deconv *= convolve2d(relative_blur, psf_mirror, 'same') + relative_blur = image / convolve_method(im_deconv, psf, 'same') + im_deconv *= convolve_method(relative_blur, psf_mirror, 'same') if clip: im_deconv[im_deconv > 1] = 1 diff --git a/skimage/restoration/non_local_means.py b/skimage/restoration/non_local_means.py index c7bb22b9..c2d56681 100644 --- a/skimage/restoration/non_local_means.py +++ b/skimage/restoration/non_local_means.py @@ -69,6 +69,8 @@ def denoise_nl_means(image, patch_size=7, patch_distance=11, h=0.1, shift, that reduces the number of operations [1]_. Therefore, this algorithm executes faster than the classic algorith (``fast_mode=False``), at the expense of using twice as much memory. + This implementation has been proven to be more efficient compared to + other alternatives, see e.g. [3]_. Compared to the classic algorithm, all pixels of a patch contribute to the distance to another patch with the same weight, no matter @@ -84,9 +86,14 @@ def denoise_nl_means(image, patch_size=7, patch_distance=11, h=0.1, References ---------- .. [1] Buades, A., Coll, B., & Morel, J. M. (2005, June). A non-local - algorithm for image denoising. In CVPR 2005, Vol. 2, pp. 60-65, IEEE. + algorithm for image denoising. In CVPR 2005, Vol. 2, pp. 60-65, IEEE. - .. [2] Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means + .. [2] J. Darbon, A. Cunha, T.F. Chan, S. Osher, and G.J. Jensen, Fast + nonlocal filtering applied to electron cryomicroscopy, in 5th IEEE + International Symposium on Biomedical Imaging: From Nano to Macro, + 2008, pp. 1331-1334. + + .. [3] Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326. Examples diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index ef70c11b..fb199afd 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -21,7 +21,7 @@ def test_denoise_tv_chambolle_2d(): # clip noise so that it does not exceed allowed range for float images. img = np.clip(img, 0, 1) # denoise - denoised_astro = restoration.denoise_tv_chambolle(img, weight=60.0) + denoised_astro = restoration.denoise_tv_chambolle(img, weight=0.25) # which dtype? assert denoised_astro.dtype in [np.float, np.float32, np.float64] from scipy import ndimage as ndi @@ -30,12 +30,12 @@ def test_denoise_tv_chambolle_2d(): # test if the total variation has decreased assert grad_denoised.dtype == np.float assert (np.sqrt((grad_denoised**2).sum()) - < np.sqrt((grad**2).sum()) / 2) + < np.sqrt((grad**2).sum())) def test_denoise_tv_chambolle_multichannel(): - denoised0 = restoration.denoise_tv_chambolle(astro[..., 0], weight=60.0) - denoised = restoration.denoise_tv_chambolle(astro, weight=60.0, + denoised0 = restoration.denoise_tv_chambolle(astro[..., 0], weight=0.25) + denoised = restoration.denoise_tv_chambolle(astro, weight=0.25, multichannel=True) assert_equal(denoised[..., 0], denoised0) @@ -46,7 +46,7 @@ def test_denoise_tv_chambolle_float_result_range(): int_astro = np.multiply(img, 255).astype(np.uint8) assert np.max(int_astro) > 1 denoised_int_astro = restoration.denoise_tv_chambolle(int_astro, - weight=60.0) + weight=0.25) # test if the value range of output float data is within [0.0:1.0] assert denoised_int_astro.dtype == np.float assert np.max(denoised_int_astro) <= 1.0 @@ -62,7 +62,7 @@ def test_denoise_tv_chambolle_3d(): mask += 20 * np.random.rand(*mask.shape) mask[mask < 0] = 0 mask[mask > 255] = 255 - res = restoration.denoise_tv_chambolle(mask.astype(np.uint8), weight=100) + res = restoration.denoise_tv_chambolle(mask.astype(np.uint8), weight=0.4) assert res.dtype == np.float assert res.std() * 255 < mask.std() diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index f79fb482..63f7d66d 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -1,4 +1,5 @@ from .random_walker_segmentation import random_walker +from .active_contour_model import active_contour from ._felzenszwalb import felzenszwalb from .slic_superpixels import slic from ._quickshift import quickshift @@ -8,6 +9,7 @@ from ._join import join_segmentations, relabel_from_one, relabel_sequential __all__ = ['random_walker', + 'active_contour', 'felzenszwalb', 'slic', 'quickshift', diff --git a/skimage/segmentation/_clear_border.py b/skimage/segmentation/_clear_border.py index 266e17e2..bc972754 100644 --- a/skimage/segmentation/_clear_border.py +++ b/skimage/segmentation/_clear_border.py @@ -1,37 +1,41 @@ import numpy as np -from scipy.ndimage import label +#from scipy.ndimage import label +from ..measure import label -def clear_border(image, buffer_size=0, bgval=0): - """Clear objects connected to image border. +def clear_border(labels, buffer_size=0, bgval=0, in_place=False): + """Clear objects connected to the label image border. - The changes will be applied to the input image. + The changes will be applied directly to the input. Parameters ---------- - image : (N, M) array - Binary image. + labels : (N, M) array of int + Label or binary image. buffer_size : int, optional - Define additional buffer around image border. + The width of the border examined. By default, only objects + that touch the outside of the image are removed. bgval : float or int, optional - Value for cleared objects. + Cleared objects are set to this value. + in_place : bool, optional + Whether or not to manipulate the labels array in-place. Returns ------- - image : (N, M) array + labels : (N, M) array Cleared binary image. Examples -------- >>> import numpy as np >>> from skimage.segmentation import clear_border - >>> image = np.array([[0, 0, 0, 0, 0, 0, 0, 1, 0], - ... [0, 0, 0, 0, 1, 0, 0, 0, 0], - ... [1, 0, 0, 1, 0, 1, 0, 0, 0], - ... [0, 0, 1, 1, 1, 1, 1, 0, 0], - ... [0, 1, 1, 1, 1, 1, 1, 1, 0], - ... [0, 0, 0, 0, 0, 0, 0, 0, 0]]) - >>> clear_border(image) + >>> labels = np.array([[0, 0, 0, 0, 0, 0, 0, 1, 0], + ... [0, 0, 0, 0, 1, 0, 0, 0, 0], + ... [1, 0, 0, 1, 0, 1, 0, 0, 0], + ... [0, 0, 1, 1, 1, 1, 1, 0, 0], + ... [0, 1, 1, 1, 1, 1, 1, 1, 0], + ... [0, 0, 0, 0, 0, 0, 0, 0, 0]]) + >>> clear_border(labels) array([[0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 1, 0, 0, 0], @@ -40,6 +44,7 @@ def clear_border(image, buffer_size=0, bgval=0): [0, 0, 0, 0, 0, 0, 0, 0, 0]]) """ + image = labels rows, cols = image.shape if buffer_size >= rows or buffer_size >= cols: @@ -53,7 +58,10 @@ def clear_border(image, buffer_size=0, bgval=0): borders[:, :ext] = True borders[:, - ext:] = True - labels, number = label(image) + # Re-label, in case we are dealing with a binary image + # and to get consistent labeling + labels = label(image, background=0) + 1 + number = np.max(labels) + 1 # determine all objects that are connected to borders borders_indices = np.unique(labels[borders]) @@ -63,6 +71,9 @@ def clear_border(image, buffer_size=0, bgval=0): # create mask for pixels to clear mask = label_mask[labels.ravel()].reshape(labels.shape) + if not in_place: + image = image.copy() + # clear border pixels image[mask] = bgval diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx index 19fc1627..bde81749 100644 --- a/skimage/segmentation/_slic.pyx +++ b/skimage/segmentation/_slic.pyx @@ -76,7 +76,8 @@ def _slic_cython(double[:, :, :, ::1] image_zyx, # approximate grid size for desired n_segments cdef Py_ssize_t step_z, step_y, step_x slices = regular_grid((depth, height, width), n_segments) - step_z, step_y, step_x = [int(s.step) for s in slices] + step_z, step_y, step_x = [int(s.step if s.step is not None else 1) + for s in slices] cdef Py_ssize_t[:, :, ::1] nearest_segments \ = np.empty((depth, height, width), dtype=np.intp) diff --git a/skimage/segmentation/active_contour_model.py b/skimage/segmentation/active_contour_model.py new file mode 100644 index 00000000..b2f89b3b --- /dev/null +++ b/skimage/segmentation/active_contour_model.py @@ -0,0 +1,234 @@ +import numpy as np +from skimage import img_as_float +import scipy +import scipy.linalg +from scipy.interpolate import RectBivariateSpline, interp2d +from skimage.filters import sobel + + +def active_contour(image, snake, alpha=0.01, beta=0.1, + w_line=0, w_edge=1, gamma=0.01, + bc='periodic', max_px_move=1.0, + max_iterations=2500, convergence=0.1): + """Active contour model. + + Active contours by fitting snakes to features of images. Supports single + and multichannel 2D images. Snakes can be periodic (for segmentation) or + have fixed and/or free ends. + + Parameters + ---------- + image : (N, M) or (N, M, 3) ndarray + Input image. + snake : (N, 2) ndarray + Initialisation coordinates of snake. For periodic snakes, it should + not include duplicate endpoints. + alpha : float, optional + Snake length shape parameter. Higher values makes snake contract + faster. + beta : float, optional + Snake smoothness shape parameter. Higher values makes snake smoother. + w_line : float, optional + Controls attraction to brightness. Use negative values to attract to + dark regions. + w_edge : float, optional + Controls attraction to edges. Use negative values to repel snake from + edges. + gamma : float, optional + Explicit time stepping parameter. + bc : {'periodic', 'free', 'fixed'}, optional + Boundary conditions for worm. 'periodic' attaches the two ends of the + snake, 'fixed' holds the end-points in place, and'free' allows free + movement of the ends. 'fixed' and 'free' can be combined by parsing + 'fixed-free', 'free-fixed'. Parsing 'fixed-fixed' or 'free-free' + yields same behaviour as 'fixed' and 'free', respectively. + max_px_move : float, optional + Maximum pixel distance to move per iteration. + max_iterations : int, optional + Maximum iterations to optimize snake shape. + convergence: float, optional + Convergence criteria. + + Returns + ------- + snake : (N, 2) ndarray + Optimised snake, same shape as input parameter. + + References + ---------- + .. [1] Kass, M.; Witkin, A.; Terzopoulos, D. "Snakes: Active contour + models". International Journal of Computer Vision 1 (4): 321 + (1988). + + Examples + -------- + >>> from skimage.draw import circle_perimeter + >>> from skimage.filters import gaussian_filter + + Create and smooth image: + + >>> img = np.zeros((100, 100)) + >>> rr, cc = circle_perimeter(35, 45, 25) + >>> img[rr, cc] = 1 + >>> img = gaussian_filter(img, 2) + + Initiliaze spline: + + >>> s = np.linspace(0, 2*np.pi,100) + >>> init = 50*np.array([np.cos(s), np.sin(s)]).T+50 + + Fit spline to image: + + >>> snake = active_contour(img, init, w_edge=0, w_line=1) #doctest: +SKIP + >>> dist = np.sqrt((45-snake[:, 0])**2 +(35-snake[:, 1])**2) #doctest: +SKIP + >>> int(np.mean(dist)) #doctest: +SKIP + 25 + + """ + scipy_version = list(map(int, scipy.__version__.split('.'))) + new_scipy = scipy_version[0] > 0 or \ + (scipy_version[0] == 0 and scipy_version[1] >= 14) + if not new_scipy: + raise NotImplementedError('You are using an old version of scipy. ' + 'Active contours is implemented for scipy versions ' + '0.14.0 and above.') + + max_iterations = int(max_iterations) + if max_iterations <= 0: + raise ValueError("max_iterations should be >0.") + convergence_order = 10 + valid_bcs = ['periodic', 'free', 'fixed', 'free-fixed', + 'fixed-free', 'fixed-fixed', 'free-free'] + if bc not in valid_bcs: + raise ValueError("Invalid boundary condition.\n" + + "Should be one of: "+", ".join(valid_bcs)+'.') + img = img_as_float(image) + RGB = img.ndim == 3 + + # Find edges using sobel: + if w_edge != 0: + if RGB: + edge = [sobel(img[:, :, 0]), sobel(img[:, :, 1]), + sobel(img[:, :, 2])] + else: + edge = [sobel(img)] + for i in range(3 if RGB else 1): + edge[i][0, :] = edge[i][1, :] + edge[i][-1, :] = edge[i][-2, :] + edge[i][:, 0] = edge[i][:, 1] + edge[i][:, -1] = edge[i][:, -2] + else: + edge = [0] + + # Superimpose intensity and edge images: + if RGB: + img = w_line*np.sum(img, axis=2) \ + + w_edge*sum(edge) + else: + img = w_line*img + w_edge*edge[0] + + # Interpolate for smoothness: + if new_scipy: + intp = RectBivariateSpline(np.arange(img.shape[1]), + np.arange(img.shape[0]), + img.T, kx=2, ky=2, s=0) + else: + intp = np.vectorize(interp2d(np.arange(img.shape[1]), + np.arange(img.shape[0]), img, kind='cubic', + copy=False, bounds_error=False, fill_value=0)) + + x, y = snake[:, 0].copy(), snake[:, 1].copy() + xsave = np.empty((convergence_order, len(x))) + ysave = np.empty((convergence_order, len(x))) + + # Build snake shape matrix for Euler equation + n = len(x) + a = np.roll(np.eye(n), -1, axis=0) + \ + np.roll(np.eye(n), -1, axis=1) - \ + 2*np.eye(n) # second order derivative, central difference + b = np.roll(np.eye(n), -2, axis=0) + \ + np.roll(np.eye(n), -2, axis=1) - \ + 4*np.roll(np.eye(n), -1, axis=0) - \ + 4*np.roll(np.eye(n), -1, axis=1) + \ + 6*np.eye(n) # fourth order derivative, central difference + A = -alpha*a + beta*b + + # Impose boundary conditions different from periodic: + sfixed = False + if bc.startswith('fixed'): + A[0, :] = 0 + A[1, :] = 0 + A[1, :3] = [1, -2, 1] + sfixed = True + efixed = False + if bc.endswith('fixed'): + A[-1, :] = 0 + A[-2, :] = 0 + A[-2, -3:] = [1, -2, 1] + efixed = True + sfree = False + if bc.startswith('free'): + A[0, :] = 0 + A[0, :3] = [1, -2, 1] + A[1, :] = 0 + A[1, :4] = [-1, 3, -3, 1] + sfree = True + efree = False + if bc.endswith('free'): + A[-1, :] = 0 + A[-1, -3:] = [1, -2, 1] + A[-2, :] = 0 + A[-2, -4:] = [-1, 3, -3, 1] + efree = True + + # Only one inversion is needed for implicit spline energy minimization: + inv = scipy.linalg.inv(A+gamma*np.eye(n)) + + # Explicit time stepping for image energy minimization: + for i in range(max_iterations): + if new_scipy: + fx = intp(x, y, dx=1, grid=False) + fy = intp(x, y, dy=1, grid=False) + else: + fx = intp(x, y, dx=1) + fy = intp(x, y, dy=1) + if sfixed: + fx[0] = 0 + fy[0] = 0 + if efixed: + fx[-1] = 0 + fy[-1] = 0 + if sfree: + fx[0] *= 2 + fy[0] *= 2 + if efree: + fx[-1] *= 2 + fy[-1] *= 2 + xn = np.dot(inv, gamma*x + fx) + yn = np.dot(inv, gamma*y + fy) + + # Movements are capped to max_px_move per iteration: + dx = max_px_move*np.tanh(xn-x) + dy = max_px_move*np.tanh(yn-y) + if sfixed: + dx[0] = 0 + dy[0] = 0 + if efixed: + dx[-1] = 0 + dy[-1] = 0 + x[:] += dx + y[:] += dy + + # Convergence criteria needs to compare to a number of previous + # configurations since oscillations can occur. + j = i % (convergence_order+1) + if j < convergence_order: + xsave[j, :] = x + ysave[j, :] = y + else: + dist = np.min(np.max(np.abs(xsave-x[None, :]) + + np.abs(ysave-y[None, :]), 1)) + if dist < convergence: + break + + return np.array([x, y]).T diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py index 3eca9401..adc21a65 100644 --- a/skimage/segmentation/slic_superpixels.py +++ b/skimage/segmentation/slic_superpixels.py @@ -155,7 +155,8 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=0, # initialize cluster centroids for desired number of segments grid_z, grid_y, grid_x = np.mgrid[:depth, :height, :width] slices = regular_grid(image.shape[:3], n_segments) - step_z, step_y, step_x = [int(s.step) for s in slices] + step_z, step_y, step_x = [int(s.step if s.step is not None else 1) + for s in slices] segments_z = grid_z[slices] segments_y = grid_y[slices] segments_x = grid_x[slices] diff --git a/skimage/segmentation/tests/test_active_contour_model.py b/skimage/segmentation/tests/test_active_contour_model.py new file mode 100644 index 00000000..48a86d45 --- /dev/null +++ b/skimage/segmentation/tests/test_active_contour_model.py @@ -0,0 +1,116 @@ +import numpy as np +from skimage import data +from skimage.color import rgb2gray +from skimage.filters import gaussian_filter +from skimage.segmentation import active_contour +from numpy.testing import assert_equal, assert_allclose, assert_raises +from numpy.testing.decorators import skipif +import scipy + +scipy_version = list(map(int, scipy.__version__.split('.'))) +new_scipy = scipy_version[0] > 0 or \ + (scipy_version[0] == 0 and scipy_version[1] >= 14) + +@skipif(not new_scipy) +def test_periodic_reference(): + img = data.astronaut() + img = rgb2gray(img) + s = np.linspace(0, 2*np.pi, 400) + x = 220 + 100*np.cos(s) + y = 100 + 100*np.sin(s) + init = np.array([x, y]).T + snake = active_contour(gaussian_filter(img, 3), init, + alpha=0.015, beta=10, w_line=0, w_edge=1, gamma=0.001) + refx = [299, 298, 298, 298, 298, 297, 297, 296, 296, 295] + refy = [98, 99, 100, 101, 102, 103, 104, 105, 106, 108] + assert_equal(np.array(snake[:10, 0], dtype=np.int32), refx) + assert_equal(np.array(snake[:10, 1], dtype=np.int32), refy) + +@skipif(not new_scipy) +def test_fixed_reference(): + img = data.text() + x = np.linspace(5, 424, 100) + y = np.linspace(136, 50, 100) + init = np.array([x, y]).T + snake = active_contour(gaussian_filter(img, 1), init, bc='fixed', + alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1) + refx = [5, 9, 13, 17, 21, 25, 30, 34, 38, 42] + refy = [136, 135, 134, 133, 132, 131, 129, 128, 127, 125] + assert_equal(np.array(snake[:10, 0], dtype=np.int32), refx) + assert_equal(np.array(snake[:10, 1], dtype=np.int32), refy) + +@skipif(not new_scipy) +def test_free_reference(): + img = data.text() + x = np.linspace(5, 424, 100) + y = np.linspace(70, 40, 100) + init = np.array([x, y]).T + snake = active_contour(gaussian_filter(img, 3), init, bc='free', + alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1) + refx = [10, 13, 16, 19, 23, 26, 29, 32, 36, 39] + refy = [76, 76, 75, 74, 73, 72, 71, 70, 69, 69] + assert_equal(np.array(snake[:10, 0], dtype=np.int32), refx) + assert_equal(np.array(snake[:10, 1], dtype=np.int32), refy) + +@skipif(not new_scipy) +def test_RGB(): + img = gaussian_filter(data.text(), 1) + imgR = np.zeros((img.shape[0], img.shape[1], 3)) + imgG = np.zeros((img.shape[0], img.shape[1], 3)) + imgRGB = np.zeros((img.shape[0], img.shape[1], 3)) + imgR[:, :, 0] = img + imgG[:, :, 1] = img + imgRGB[:, :, :] = img[:, :, None] + x = np.linspace(5, 424, 100) + y = np.linspace(136, 50, 100) + init = np.array([x, y]).T + snake = active_contour(imgR, init, bc='fixed', + alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1) + refx = [5, 9, 13, 17, 21, 25, 30, 34, 38, 42] + refy = [136, 135, 134, 133, 132, 131, 129, 128, 127, 125] + assert_equal(np.array(snake[:10, 0], dtype=np.int32), refx) + assert_equal(np.array(snake[:10, 1], dtype=np.int32), refy) + snake = active_contour(imgG, init, bc='fixed', + alpha=0.1, beta=1.0, w_line=-5, w_edge=0, gamma=0.1) + assert_equal(np.array(snake[:10, 0], dtype=np.int32), refx) + assert_equal(np.array(snake[:10, 1], dtype=np.int32), refy) + snake = active_contour(imgRGB, init, bc='fixed', + alpha=0.1, beta=1.0, w_line=-5/3., w_edge=0, gamma=0.1) + assert_equal(np.array(snake[:10, 0], dtype=np.int32), refx) + assert_equal(np.array(snake[:10, 1], dtype=np.int32), refy) + +@skipif(not new_scipy) +def test_end_points(): + img = data.astronaut() + img = rgb2gray(img) + s = np.linspace(0, 2*np.pi, 400) + x = 220 + 100*np.cos(s) + y = 100 + 100*np.sin(s) + init = np.array([x, y]).T + snake = active_contour(gaussian_filter(img, 3), init, + bc='periodic', alpha=0.015, beta=10, w_line=0, w_edge=1, + gamma=0.001, max_iterations=100) + assert np.sum(np.abs(snake[0, :]-snake[-1, :])) < 2 + snake = active_contour(gaussian_filter(img, 3), init, + bc='free', alpha=0.015, beta=10, w_line=0, w_edge=1, + gamma=0.001, max_iterations=100) + assert np.sum(np.abs(snake[0, :]-snake[-1, :])) > 2 + snake = active_contour(gaussian_filter(img, 3), init, + bc='fixed', alpha=0.015, beta=10, w_line=0, w_edge=1, + gamma=0.001, max_iterations=100) + assert_allclose(snake[0, :], [x[0], y[0]], atol=1e-5) + +@skipif(not new_scipy) +def test_bad_input(): + img = np.zeros((10, 10)) + x = np.linspace(5, 424, 100) + y = np.linspace(136, 50, 100) + init = np.array([x, y]).T + assert_raises(ValueError, active_contour, img, init, + bc='wrong') + assert_raises(ValueError, active_contour, img, init, + max_iterations=-15) + + +if __name__ == "__main__": + np.testing.run_module_suite() diff --git a/skimage/segmentation/tests/test_clear_border.py b/skimage/segmentation/tests/test_clear_border.py index 5d6852cf..b626ea2b 100644 --- a/skimage/segmentation/tests/test_clear_border.py +++ b/skimage/segmentation/tests/test_clear_border.py @@ -1,5 +1,5 @@ import numpy as np -from numpy.testing import assert_array_equal +from numpy.testing import assert_array_equal, assert_ from skimage.segmentation import clear_border @@ -24,9 +24,40 @@ def test_clear_border(): assert_array_equal(result, np.zeros(result.shape)) # test background value - result = clear_border(image.copy(), 1, 2) + result = clear_border(image.copy(), buffer_size=1, bgval=2) assert_array_equal(result, 2 * np.ones_like(image)) +def test_clear_border_non_binary(): + image = np.array([[1, 2, 3, 1, 2], + [3, 3, 5, 4, 2], + [3, 4, 5, 4, 2], + [3, 3, 2, 1, 2]]) + + result = clear_border(image) + expected = np.array([[0, 0, 0, 0, 0], + [0, 0, 5, 4, 0], + [0, 4, 5, 4, 0], + [0, 0, 0, 0, 0]]) + + assert_array_equal(result, expected) + assert_(not np.all(image == result)) + + +def test_clear_border_non_binary_inplace(): + image = np.array([[1, 2, 3, 1, 2], + [3, 3, 5, 4, 2], + [3, 4, 5, 4, 2], + [3, 3, 2, 1, 2]]) + + result = clear_border(image, in_place=True) + expected = np.array([[0, 0, 0, 0, 0], + [0, 0, 5, 4, 0], + [0, 4, 5, 4, 0], + [0, 0, 0, 0, 0]]) + + assert_array_equal(result, expected) + assert_array_equal(image, result) + if __name__ == "__main__": np.testing.run_module_suite() diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py index 73629103..83424543 100644 --- a/skimage/segmentation/tests/test_slic.py +++ b/skimage/segmentation/tests/test_slic.py @@ -195,6 +195,20 @@ def test_slic_zero(): assert_equal(seg[10:, 10:], 3) +def test_more_segments_than_pixels(): + rnd = np.random.RandomState(0) + img = np.zeros((20, 21)) + img[:10, :10] = 0.33 + img[10:, :10] = 0.67 + img[10:, 10:] = 1.00 + img += 0.0033 * rnd.normal(size=img.shape) + img[img > 1] = 1 + img[img < 0] = 0 + seg = slic(img, sigma=0, n_segments=500, compactness=1, + multichannel=False, convert2lab=False) + assert np.all(seg.ravel() == np.arange(seg.size)) + + if __name__ == '__main__': from numpy import testing diff --git a/skimage/transform/_seam_carving.pyx b/skimage/transform/_seam_carving.pyx index 3499514c..73306bb3 100644 --- a/skimage/transform/_seam_carving.pyx +++ b/skimage/transform/_seam_carving.pyx @@ -180,7 +180,7 @@ def _seam_carve_v(img, energy_map, iters, border): cdef Py_ssize_t seams_left = iters cdef Py_ssize_t seams_removed cdef Py_ssize_t seam_idx - cdef Py_ssize_t[::1] seam_buffer = np.zeros(rows, dtype=np.int) + cdef Py_ssize_t[::1] seam_buffer = np.zeros(rows, dtype=np.intp) cdef cnp.double_t[:, :, ::1] image = img cdef cnp.int8_t[:, ::1] track_img = np.zeros(img.shape[0:2], dtype=np.int8) diff --git a/skimage/transform/integral.py b/skimage/transform/integral.py index 3503a652..91b07ac7 100644 --- a/skimage/transform/integral.py +++ b/skimage/transform/integral.py @@ -1,7 +1,8 @@ import numpy as np +import collections +import warnings - -def integral_image(x): +def integral_image(img): """Integral image / summed area table. The integral image contains the sum of all elements above and to the @@ -13,13 +14,13 @@ def integral_image(x): Parameters ---------- - x : ndarray + img : ndarray Input image. Returns ------- S : ndarray - Integral image / summed area table. + Integral image/summed area table of same shape as input image. References ---------- @@ -27,44 +28,121 @@ def integral_image(x): ACM SIGGRAPH Computer Graphics, vol. 18, 1984, pp. 207-212. """ - return x.cumsum(1).cumsum(0) + S = img + for i in range(img.ndim): + S = S.cumsum(axis=i) + return S -def integrate(ii, r0, c0, r1, c1): +def integrate(ii, start, end, *args): """Use an integral image to integrate over a given window. Parameters ---------- ii : ndarray Integral image. - r0, c0 : int or ndarray - Top-left corner(s) of block to be summed. - r1, c1 : int or ndarray - Bottom-right corner(s) of block to be summed. + start : List of tuples, each tuple of length equal to dimension of `ii` + Coordinates of top left corner of window(s). + Each tuple in the list contains the starting row, col, ... index + i.e `[(row_win1, col_win1, ...), (row_win2, col_win2,...), ...]`. + end : List of tuples, each tuple of length equal to dimension of `ii` + Coordinates of bottom right corner of window(s). + Each tuple in the list containing the end row, col, ... index i.e + `[(row_win1, col_win1, ...), (row_win2, col_win2, ...), ...]`. + args: optional + For backward compatibility with versions prior to 0.12. + The earlier function signature was `integrate(ii, r0, c0, r1, c1)`, + where `r0`, `c0` are int(lists) specifying start coordinates + of window(s) to be integrated and `r1`, `c1` the end coordinates. Returns ------- S : scalar or ndarray Integral (sum) over the given window(s). + + Examples + -------- + >>> arr = np.ones((5, 6), dtype=np.float) + >>> ii = integral_image(arr) + >>> integrate(ii, (1, 0), (1, 2)) # sum from (1, 0) to (1, 2) + array([ 3.]) + >>> integrate(ii, [(3, 3)], [(4, 5)]) # sum from (3, 3) to (4, 5) + array([ 6.]) + >>> # sum from (1, 0) to (1, 2) and from (3, 3) to (4, 5) + >>> integrate(ii, [(1, 0), (3, 3)], [(1, 2), (4, 5)]) + array([ 3., 6.]) """ - if np.isscalar(r0): - r0, c0, r1, c1 = [np.asarray([x]) for x in (r0, c0, r1, c1)] + rows = 1 + # handle input from new input format + if len(args) == 0: + start = np.atleast_2d(np.array(start)) + end = np.atleast_2d(np.array(end)) + rows = start.shape[0] + # handle deprecated input format + else: + warnings.warn("The syntax 'integrate(ii, r0, c0, r1, c1)' is " + "deprecated, and will be phased out in release 0.14. " + "The new syntax is " + "'integrate(ii, (r0, c0), (r1, c1))'.") + if isinstance(start, collections.Iterable): + rows = len(start) + args = (start, end) + args + start = np.array(args[:int(len(args)/2)]).T + end = np.array(args[int(len(args)/2):]).T - S = np.zeros(r0.shape, ii.dtype) + total_shape = ii.shape + total_shape = np.tile(total_shape, [rows, 1]) - S += ii[r1, c1] + # convert negative indices into equivalent positive indices + start_negatives = start < 0 + end_negatives = end < 0 + start = (start + total_shape) * start_negatives + \ + start * ~(start_negatives) + end = (end + total_shape) * end_negatives + \ + end * ~(end_negatives) - good = (r0 >= 1) & (c0 >= 1) - S[good] += ii[r0[good] - 1, c0[good] - 1] + if np.any((end - start) < 0): + raise IndexError('end coordinates must be greater or equal to start') - good = r0 >= 1 - S[good] -= ii[r0[good] - 1, c1[good]] + # bit_perm is the total number of terms in the expression + # of S. For example, in the case of a 4x4 2D image + # sum of image from (1,1) to (2,2) is given by + # S = + ii[2, 2] + # - ii[0, 2] - ii[2, 0] + # + ii[0, 0] + # The total terms = 4 = 2 ** 2(dims) - good = c0 >= 1 - S[good] -= ii[r1[good], c0[good] - 1] + S = np.zeros(rows) + bit_perm = 2 ** ii.ndim + width = len(bin(bit_perm - 1)[2:]) - if S.size == 1: - return np.asscalar(S) + # Sum of a (hyper)cube, from an integral image is computed using + # values at the corners of the cube. The corners of cube are + # selected using binary numbers as described in the following example. + # In a 3D cube there are 8 corners. The corners are selected using + # binary numbers 000 to 111. Each number is called a permutation, where + # perm(000) means, select end corner where none of the coordinates + # is replaced, i.e ii[end_row, end_col, end_depth]. Similarly, perm(001) + # means replace last coordinate by start - 1, i.e + # ii[end_row, end_col, start_depth - 1], and so on. + # Sign of even permutations is positive, while those of odd is negative. + # If 'start_coord - 1' is -ve it is labeled bad and not considered in + # the final sum. + for i in range(bit_perm): # for all permutations + # boolean permutation array eg [True, False] for '10' + binary = bin(i)[2:].zfill(width) + bool_mask = [bit == '1' for bit in binary] + + sign = (-1)**sum(bool_mask) # determine sign of permutation + + bad = [np.any(((start[r] - 1) * bool_mask) < 0) + for r in range(rows)] # find out bad start rows + + corner_points = (end * (np.invert(bool_mask))) + \ + ((start - 1) * bool_mask) # find corner for each row + + S += [sign * ii[tuple(corner_points[r])] if(not bad[r]) else 0 + for r in range(rows)] # add only good rows return S diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index e49496c1..f71b7566 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -20,6 +20,7 @@ from scipy.interpolate import interp1d from ._warps_cy import _warp_fast from ._radon_transform import sart_projection_update from .. import util +from warnings import warn __all__ = ["radon", "iradon", "iradon_sart"] @@ -49,11 +50,6 @@ def radon(image, theta=None, circle=False): at the pixel index ``radon_image.shape[0] // 2`` along the 0th dimension of ``radon_image``. - Raises - ------ - ValueError - If called with ``circle=True`` and ``image != 0`` outside the inscribed - circle """ if image.ndim != 2: raise ValueError('The input image must be 2-D') @@ -67,8 +63,8 @@ def radon(image, theta=None, circle=False): + (c1 - image.shape[1] // 2) ** 2) reconstruction_circle = reconstruction_circle <= radius ** 2 if not np.all(reconstruction_circle | (image == 0)): - raise ValueError('Image must be zero outside the reconstruction' - ' circle') + warn('Radon transform: image must be zero outside the ' + 'reconstruction circle') # Crop image to make it square slices = [] for d in (0, 1): diff --git a/skimage/transform/tests/test_integral.py b/skimage/transform/tests/test_integral.py index eb92a40b..e05ed397 100644 --- a/skimage/transform/tests/test_integral.py +++ b/skimage/transform/tests/test_integral.py @@ -17,15 +17,16 @@ def test_validity(): def test_basic(): - assert_equal(x[12:24, 10:20].sum(), integrate(s, 12, 10, 23, 19)) - assert_equal(x[:20, :20].sum(), integrate(s, 0, 0, 19, 19)) - assert_equal(x[:20, 10:20].sum(), integrate(s, 0, 10, 19, 19)) - assert_equal(x[10:20, :20].sum(), integrate(s, 10, 0, 19, 19)) + assert_equal(x[12:24, 10:20].sum(), integrate(s, (12, 10), (23, 19))) + assert_equal(x[:20, :20].sum(), integrate(s, (0, 0), (19, 19))) + assert_equal(x[:20, 10:20].sum(), integrate(s, (0, 10), (19, 19))) + assert_equal(x[10:20, :20].sum(), integrate(s, (10, 0), (19, 19))) def test_single(): - assert_equal(x[0, 0], integrate(s, 0, 0, 0, 0)) - assert_equal(x[10, 10], integrate(s, 10, 10, 10, 10)) + assert_equal(x[0, 0], integrate(s, (0, 0), (0, 0))) + assert_equal(x[10, 10], integrate(s, (10, 10), (10, 10))) + def test_vectorized_integrate(): r0 = np.array([12, 0, 0, 10, 0, 10, 30]) @@ -40,7 +41,10 @@ def test_vectorized_integrate(): x[0,0], x[10, 10], x[30:, 31:].sum()]) - assert_equal(expected, integrate(s, r0, c0, r1, c1)) + start_pts = [(r0[i], c0[i]) for i in range(len(r0))] + end_pts = [(r1[i], c1[i]) for i in range(len(r0))] + assert_equal(expected, integrate(s, r0, c0, r1, c1)) # test deprecated + assert_equal(expected, integrate(s, start_pts, end_pts)) if __name__ == '__main__': diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py index 6df8b7f8..63bfa31f 100644 --- a/skimage/transform/tests/test_radon_transform.py +++ b/skimage/transform/tests/test_radon_transform.py @@ -9,7 +9,7 @@ from skimage.transform import radon, iradon, iradon_sart, rescale from skimage.io import imread from skimage import data_dir from skimage._shared.testing import test_parallel - +from skimage._shared._warnings import expected_warnings PHANTOM = imread(os.path.join(data_dir, "phantom.png"), as_grey=True)[::2, ::2] @@ -215,7 +215,8 @@ def _random_circle(shape): def test_radon_circle(): a = np.ones((10, 10)) - assert_raises(ValueError, radon, a, circle=True) + with expected_warnings(['reconstruction circle']): + radon(a, circle=True) # Synthetic data, circular symmetry shape = (61, 79) diff --git a/skimage/util/dtype.py b/skimage/util/dtype.py index 1c594534..077b14be 100644 --- a/skimage/util/dtype.py +++ b/skimage/util/dtype.py @@ -21,8 +21,8 @@ dtype_range = {np.bool_: (False, True), integer_types = (np.uint8, np.uint16, np.int8, np.int16) _supported_types = (np.bool_, np.bool8, - np.uint8, np.uint16, np.uint32, - np.int8, np.int16, np.int32, + np.uint8, np.uint16, np.uint32, np.uint64, + np.int8, np.int16, np.int32, np.int64, np.float32, np.float64) if np.__version__ >= "1.6.0": @@ -125,7 +125,18 @@ def convert(image, dtype, force_copy=False, uniform=False): # Numbers can be represented exactly only if m is a multiple of n # Output array is of same kind as input. kind = a.dtype.kind - if n == m: + if n > m and a.max() < 2 ** m: + mnew = int(np.ceil(m / 2) * 2) + if mnew > m: + dtype = "int%s" % mnew + else: + dtype = "uint%s" % mnew + n = int(np.ceil(n / 2) * 2) + msg = ("Downcasting %s to %s without scaling because max " + "value %s fits in %s" % (a.dtype, dtype, a.max(), dtype)) + warn(msg) + return a.astype(_dtype2(kind, m)) + elif n == m: return a.copy() if copy else a elif n > m: # downscale with precision loss diff --git a/skimage/util/tests/test_dtype.py b/skimage/util/tests/test_dtype.py index 612c43e6..b8e77543 100644 --- a/skimage/util/tests/test_dtype.py +++ b/skimage/util/tests/test_dtype.py @@ -29,7 +29,7 @@ def test_range(): (img_as_float, np.float64), (img_as_uint, np.uint16), (img_as_ubyte, np.ubyte)]: - + with expected_warnings(['precision loss|sign loss|\A\Z']): y = f(x) @@ -62,7 +62,7 @@ def test_range_extra_dtypes(): for dtype_in, dt in dtype_pairs: imin, imax = dtype_range_extra[dtype_in] x = np.linspace(imin, imax, 10).astype(dtype_in) - + with expected_warnings(['precision loss|sign loss|\A\Z']): y = convert(x, dt) @@ -72,9 +72,12 @@ def test_range_extra_dtypes(): y, omin, omax, np.dtype(dt)) -def test_unsupported_dtype(): +def test_downcast(): x = np.arange(10).astype(np.uint64) - assert_raises(ValueError, img_as_int, x) + with expected_warnings('Downcasting'): + y = img_as_int(x) + assert np.allclose(y, x.astype(np.int16)) + assert y.dtype == np.int16, y.dtype def test_float_out_of_range(): diff --git a/skimage/viewer/tests/test_tools.py b/skimage/viewer/tests/test_tools.py index b79a40e2..64eb7fac 100644 --- a/skimage/viewer/tests/test_tools.py +++ b/skimage/viewer/tests/test_tools.py @@ -147,7 +147,7 @@ def test_rect_tool(): do_event(viewer, 'mouse_press', xdata=100, ydata=100) do_event(viewer, 'move', xdata=120, ydata=120) do_event(viewer, 'mouse_release') - assert_equal(tool.geometry, [120, 150, 120, 150]) + #assert_equal(tool.geometry, [120, 150, 120, 150]) # create a new line do_event(viewer, 'mouse_press', xdata=10, ydata=10) diff --git a/tools/mailmap.py b/tools/mailmap.py new file mode 100755 index 00000000..30f2463c --- /dev/null +++ b/tools/mailmap.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python +# Requires package 'editdistance' + +# A mailmap file is used (by GitHub and other tools) to associate multiple +# commit emails with one user. This helps to count number of commits, +# contributors, etc. + +import subprocess +import shlex +import numpy as np +from collections import defaultdict + +from editdistance import eval as dist + +threshold = 5 + +def call(cmd): + return subprocess.check_output(shlex.split(cmd), universal_newlines=True).split('\n') + + +def _clean_email(email): + if not '@' in email: + return + + name, domain = email.split('@') + name = name.split('+', 1)[0] + + return '{}@{}'.format(name, domain).lower() + + +call("rm -f .mailmap") +authors = call("git log --format='%aN::%aE'") + +names, emails = [], [] + +for (name, email) in (author.split('::') for author in authors if author.strip()): + if email not in emails: + names.append(name) + emails.append(email) + +N = len(names) +D = np.zeros((N, N)) + np.infty + +for i in range(1, N): + for j in range(i): + D[i, j] = dist(names[i], names[j]) + +for i in range(N): + dupes, = np.where(D[:, i] < threshold) + for j in dupes: + names[j] = names[i] + +mailmap = defaultdict(set) +for (name, email) in zip(names, emails): + email = _clean_email(email) + if email: + mailmap[name].add(email) + +for key, value in list(mailmap.items()): + if len(value) < 2 or (len(key.split()) < 2): + mailmap.pop(key) + +entries = [] +for name, emails in mailmap.items(): + entries.append([name]) + entries[-1].extend(['<{}>'.format(email) for email in emails]) + +entries = sorted(entries, key=lambda x: x[0].split()[-1]) +for entry in entries: + print(' '.join(entry)) diff --git a/tools/travis_before_install.sh b/tools/travis_before_install.sh index 159fc2d3..b4baeaf4 100755 --- a/tools/travis_before_install.sh +++ b/tools/travis_before_install.sh @@ -33,16 +33,15 @@ retry () { echo "cython>=0.21" >> requirements.txt # require networkx 1.9.1 on 2.6, as 2.6 support was dropped in 1.10 +# require matplotlib 1.4.3 on 2.6, as 2.6 support was dropped in 1.5 if [[ $TRAVIS_PYTHON_VERSION == 2.6* ]]; then sed -i 's/networkx.*/networkx==1.9.1/g' requirements.txt + sed -i 's/matplotlib.*/matplotlib==1.4.3/g' requirements.txt fi # test minimum requirements on 2.7 if [[ $TRAVIS_PYTHON_VERSION == 2.7* ]]; then sed -i 's/>=/==/g' requirements.txt - # PIL instead of Pillow - sed -i 's/pillow.*/pil==1.1.7/g' requirements.txt - WHEELBINARIES=${WHEELBINARIES/pillow/pil} fi # create new empty venv diff --git a/tools/travis_script.sh b/tools/travis_script.sh index ad3afca3..e08fd689 100755 --- a/tools/travis_script.sh +++ b/tools/travis_script.sh @@ -69,7 +69,7 @@ touch $MPL_DIR/matplotlibrc echo 'backend : Template' > $MPL_DIR/matplotlibrc -for f in doc/examples/*.py; do +for f in doc/examples/*/*.py; do python "$f" if [ $? -ne 0 ]; then exit 1 @@ -81,7 +81,7 @@ section_end "Run.doc.examples" section "Run.doc.applications" -for f in doc/examples/applications/*.py; do +for f in doc/examples/xx_applications/*.py; do python "$f" if [ $? -ne 0 ]; then exit 1 @@ -96,7 +96,7 @@ else MPL_QT_API=PySide export QT_API=pyside fi -echo 'backend: Agg' > $MPL_DIR/matplotlibrc +echo 'backend: Qt4Agg' > $MPL_DIR/matplotlibrc echo 'backend.qt4 : '$MPL_QT_API >> $MPL_DIR/matplotlibrc section_end "Run.doc.applications" @@ -111,3 +111,7 @@ fi nosetests $TEST_ARGS section_end "Test.with.optional.dependencies" + +section "Prepare.release" +doc/release/contribs.py HEAD~10 +section_end "Prepare.release"