diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 689f552b..f8f78c52 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -132,7 +132,8 @@ Dense DAISY feature description, circle perimeter drawing. - François Boulogne - Drawing: Andres Method for circle perimeter, ellipse perimeter drawing, Bezier curve. + Drawing: Andres Method for circle perimeter, ellipse perimeter drawing, + Bezier curve, anti-aliasing. Circular and elliptical Hough Transforms Various fixes @@ -154,3 +155,6 @@ - Riaan van den Dool skimage.io plugin: GDAL + +- Fedor Morozov + Drawing: Wu's anti-aliased circle diff --git a/MANIFEST.in b/MANIFEST.in index 97fedda6..a2ef2eaf 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -3,7 +3,7 @@ include setup*.py include MANIFEST.in include *.txt include Makefile -recursive-include skimage *.pyx *.pxd *.pxi *.py *.c *.h *.ini *.md5 +recursive-include skimage *.pyx *.pxd *.pxi *.py *.c *.h *.ini *.md5 *.rst *.txt recursive-include skimage/data * include doc/Makefile @@ -12,3 +12,4 @@ recursive-include doc/tools *.txt recursive-include doc/source/_templates *.html recursive-include doc *.py prune doc/build +prune doc/gh-pages diff --git a/RELEASE.txt b/RELEASE.txt index d29143aa..5772cf77 100644 --- a/RELEASE.txt +++ b/RELEASE.txt @@ -5,21 +5,21 @@ How to make a new release of ``skimage`` - Update release notes. - - To show a list contributors, run ``doc/release/contributors.sh ``, - where ```` is the first commit since the previous release. + - To show a list of contributors and changes, run + ``doc/release/contribs.py ``. - Update the version number in ``setup.py`` and ``bento.info`` and commit - Update the docs: - - Edit ``doc/source/themes/agogo/static/docversions.js`` and commit + - Edit ``doc/source/_static/docversions.js`` and commit - Build a clean version of the docs. Run ``make`` in the root dir, then ``rm -rf build; make html`` in the docs. - Run ``make html`` again to copy the newly generated ``random.js`` into place. Double check ``random.js``, otherwise the skimage.org front page gets broken! - Build using ``make gh-pages``. - - Push upstream: ``git push`` in ``doc/gh-pages``. + - Push upstream: ``git push origin gh-pages`` in ``doc/gh-pages``. - Add the version number as a tag in git:: diff --git a/TODO.txt b/TODO.txt index 502961c8..3cbf4b53 100644 --- a/TODO.txt +++ b/TODO.txt @@ -4,12 +4,12 @@ Version 0.10 * Remove deprecated parameter `epsilon` of `skimage.viewer.LineProfile` * Remove backwards-compatability of `skimage.measure.regionprops` * Remove {`ratio`, `sigma`} deprecation warnings of `skimage.segmentation.slic` +* Change default mode of random_walker segmentation to 'cg_mg' > 'cg' > 'bf', + depending on which optional dependencies are available. +* Remove deprecated `out` parameter of `skimage.morphology.binary_*` +* Remove deprecated parameter `depth` in `skimage.segmentation.random_walker` +* Remove deprecated logger function in `skimage/__init__.py` +* Remove deprecated function `filter.median_filter` +* Remove deprecated `skimage.color.is_gray` and `skimage.color.is_rgb` + functions -Version 0.9 ------------ -* Remove deprecated functions - - `skimage.filter.denoise_tv_chambolle` - - `skimage.morphology.is_local_maximum` - - `skimage.transform.hough` - - `skimage.transform.probabilistic_hough` - - `skimage.transform.hough_peaks` diff --git a/bento.info b/bento.info index 808c0d09..b01a09d1 100644 --- a/bento.info +++ b/bento.info @@ -1,5 +1,5 @@ Name: scikit-image -Version: 0.8.1 +Version: 0.9.3 Summary: Image processing routines for SciPy Url: http://scikit-image.org DownloadUrl: http://github.com/scikit-image/scikit-image diff --git a/doc/examples/plot_censure_keypoints.py b/doc/examples/plot_censure_keypoints.py deleted file mode 100644 index fbba4e1b..00000000 --- a/doc/examples/plot_censure_keypoints.py +++ /dev/null @@ -1,54 +0,0 @@ -""" -========================= -CenSurE Feature Detection -========================= - -In this example we detect and plot the CenSurE (Center Surround Extrema) -features at various scales using Difference of Boxes, Octagon and Star shaped -bi-level filters. - -""" - -from skimage.feature import keypoints_censure -from skimage.data import lena -from skimage.color import rgb2gray -import matplotlib.pyplot as plt - -# Initializing the parameters for Censure keypoints -img = lena() -gray_img = rgb2gray(img) -min_scale = 2 -max_scale = 6 -non_max_threshold = 0.15 -line_threshold = 10 - - -_, ax = plt.subplots(nrows=(max_scale - min_scale - 1), ncols=3, - figsize=(6, 6)) -plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.94, - bottom=0.02, left=0.06, right=0.98) - -# Detecting Censure keypoints for the following filters -for col, mode in enumerate(['dob', 'octagon', 'star']): - - ax[0, col].set_title(mode.upper(), fontsize=12) - - keypoints, scales = keypoints_censure(gray_img, min_scale, max_scale, - mode, non_max_threshold, - line_threshold) - - # Plotting Censure features at all the scales - for row, scale in enumerate(range(min_scale + 1, max_scale)): - mask = scales == scale - x = keypoints[mask, 1] - y = keypoints[mask, 0] - s = 0.5 * 2 ** (scale + min_scale + 1) - ax[row, col].imshow(img) - ax[row, col].scatter(x, y, s, facecolors='none', edgecolors='b') - ax[row, col].set_xticks([]) - ax[row, col].set_yticks([]) - ax[row, col].axis((0, img.shape[1], img.shape[0], 0)) - if col == 0: - ax[row, col].set_ylabel('Scale %d' % scale, fontsize=12) - -plt.show() diff --git a/doc/examples/plot_circular_elliptical_hough_transform.py b/doc/examples/plot_circular_elliptical_hough_transform.py index bd036e03..7fb67046 100755 --- a/doc/examples/plot_circular_elliptical_hough_transform.py +++ b/doc/examples/plot_circular_elliptical_hough_transform.py @@ -74,7 +74,6 @@ for idx in np.argsort(accums)[::-1][:5]: image[cy, cx] = (220, 20, 20) ax.imshow(image, cmap=plt.cm.gray) -plt.show() """ @@ -96,13 +95,13 @@ an ellipse passes to them. A good match corresponds to high accumulator values. A full description of the algorithm can be found in reference [1]_. - References ---------- .. [1] Xie, Yonghong, and Qiang Ji. "A new efficient ellipse detection method." Pattern Recognition, 2002. Proceedings. 16th International Conference on. Vol. 2. IEEE, 2002 """ + import matplotlib.pyplot as plt from skimage import data, filter, color @@ -110,7 +109,7 @@ from skimage.transform import hough_ellipse from skimage.draw import ellipse_perimeter # Load picture, convert to grayscale and detect edges -image_rgb = data.load('coffee.png')[0:220, 100:450] +image_rgb = data.coffee()[0:220, 160:420] image_gray = color.rgb2gray(image_rgb) edges = filter.canny(image_gray, sigma=2.0, low_threshold=0.55, high_threshold=0.8) @@ -119,29 +118,31 @@ edges = filter.canny(image_gray, sigma=2.0, # The accuracy corresponds to the bin size of a major axis. # The value is chosen in order to get a single high accumulator. # The threshold eliminates low accumulators -accum = hough_ellipse(edges, accuracy=10, threshold=170, min_size=50) -accum.sort(key=lambda x:x[5]) +result = hough_ellipse(edges, accuracy=20, threshold=250, + min_size=100, max_size=120) +result.sort(order='accumulator') + # Estimated parameters for the ellipse -center_y = int(accum[-1][0]) -center_x = int(accum[-1][1]) -xradius = int(accum[-1][2]) -yradius = int(accum[-1][3]) -angle = np.pi - accum[-1][4] +best = result[-1] +yc = int(best[1]) +xc = int(best[2]) +a = int(best[3]) +b = int(best[4]) +orientation = best[5] # Draw the ellipse on the original image -cx, cy = ellipse_perimeter(center_y, center_x, - yradius, xradius, orientation=angle) -image_rgb[cy, cx] = (0, 0, 1) +cy, cx = ellipse_perimeter(yc, xc, a, b, orientation) +image_rgb[cy, cx] = (0, 0, 255) # Draw the edge (white) and the resulting ellipse (red) edges = color.gray2rgb(edges) edges[cy, cx] = (250, 0, 0) -fig = plt.subplots(figsize=(10, 6)) -plt.subplot(1, 2, 1) -plt.title('Original picture') -plt.imshow(image_rgb) -plt.subplot(1, 2, 2) -plt.title('Edge (white) and result (red)') -plt.imshow(edges) +fig2, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(10, 6)) + +ax1.set_title('Original picture') +ax1.imshow(image_rgb) + +ax2.set_title('Edge (white) and result (red)') +ax2.imshow(edges) plt.show() diff --git a/doc/examples/plot_join_segmentations.py b/doc/examples/plot_join_segmentations.py index 2cafab15..8ccc5038 100644 --- a/doc/examples/plot_join_segmentations.py +++ b/doc/examples/plot_join_segmentations.py @@ -18,27 +18,23 @@ from skimage.filter import sobel from skimage.segmentation import slic, join_segmentations from skimage.morphology import watershed from skimage.color import label2rgb -from skimage import data +from skimage import data, img_as_float - -coins = data.coins() +coins = img_as_float(data.coins()) # make segmentation using edge-detection and watershed edges = sobel(coins) markers = np.zeros_like(coins) foreground, background = 1, 2 -markers[coins < 30] = background -markers[coins > 150] = foreground +markers[coins < 30.0 / 255] = background +markers[coins > 150.0 / 255] = foreground ws = watershed(edges, markers) seg1 = nd.label(ws == foreground)[0] # make segmentation using SLIC superpixels - -# make the RGB equivalent of `coins` -coins_colour = np.tile(coins[..., np.newaxis], (1, 1, 3)) -seg2 = slic(coins_colour, n_segments=30, max_iter=160, sigma=1, ratio=9, - convert2lab=False) +seg2 = slic(coins, n_segments=117, max_iter=160, sigma=1, compactness=0.75, + multichannel=False) # combine the two segj = join_segmentations(seg1, seg2) diff --git a/doc/examples/plot_marching_cubes.py b/doc/examples/plot_marching_cubes.py index a57a40a2..5dad1680 100644 --- a/doc/examples/plot_marching_cubes.py +++ b/doc/examples/plot_marching_cubes.py @@ -17,12 +17,11 @@ a mesh for regions of bone or bone-like density. This implementation also works correctly on anisotropic datasets, where the voxel spacing is not equal for every spatial dimension, through use of the -`sampling` kwarg. +`spacing` kwarg. """ import numpy as np import matplotlib.pyplot as plt -from mpl_toolkits.mplot3d import Axes3D from mpl_toolkits.mplot3d.art3d import Poly3DCollection from skimage import measure diff --git a/doc/examples/plot_shapes.py b/doc/examples/plot_shapes.py index 507452ef..2ae842a5 100644 --- a/doc/examples/plot_shapes.py +++ b/doc/examples/plot_shapes.py @@ -1,29 +1,34 @@ """ -=========== -Fill shapes -=========== - -This example shows how to fill several different shapes: +====== +Shapes +====== +This example shows how to draw several different shapes: * line +* Bezier curve * polygon * circle * ellipse """ +import math import numpy as np import matplotlib.pyplot as plt -from skimage.draw import line, polygon, circle, circle_perimeter, \ - ellipse, ellipse_perimeter -import numpy as np -import math +from skimage.draw import (line, polygon, circle, + circle_perimeter, + ellipse, ellipse_perimeter, + bezier_curve) -img = np.zeros((500, 500, 3), dtype=np.uint8) + +fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(10, 6)) + + +img = np.zeros((500, 500, 3), dtype=np.double) # draw line rr, cc = line(120, 123, 20, 400) -img[rr,cc,0] = 255 +img[rr, cc, 0] = 255 # fill polygon poly = np.array(( @@ -33,28 +38,61 @@ poly = np.array(( (220, 590), (300, 300), )) -rr, cc = polygon(poly[:,0], poly[:,1], img.shape) -img[rr,cc,1] = 255 +rr, cc = polygon(poly[:, 0], poly[:, 1], img.shape) +img[rr, cc, 1] = 1 # fill circle rr, cc = circle(200, 200, 100, img.shape) -img[rr,cc,:] = (255, 255, 0) +img[rr, cc, :] = (1, 1, 0) # fill ellipse rr, cc = ellipse(300, 300, 100, 200, img.shape) -img[rr,cc,2] = 255 +img[rr, cc, 2] = 1 # circle rr, cc = circle_perimeter(120, 400, 15) -img[rr, cc, :] = (255, 0, 0) +img[rr, cc, :] = (1, 0, 0) + +# Bezier curve +rr, cc = bezier_curve(70, 100, 10, 10, 150, 100, 1) +img[rr, cc, :] = (1, 0, 0) # ellipses rr, cc = ellipse_perimeter(120, 400, 60, 20, orientation=math.pi / 4.) -img[rr, cc, :] = (255, 0, 255) +img[rr, cc, :] = (1, 0, 1) rr, cc = ellipse_perimeter(120, 400, 60, 20, orientation=-math.pi / 4.) -img[rr, cc, :] = (0, 0, 255) +img[rr, cc, :] = (0, 0, 1) rr, cc = ellipse_perimeter(120, 400, 60, 20, orientation=math.pi / 2.) -img[rr, cc, :] = (255, 255, 255) +img[rr, cc, :] = (1, 1, 1) + +ax1.imshow(img) +ax1.set_title('No anti-aliasing') +ax1.axis('off') + +""" + +Anti-aliased drawing for: +* line +* circle + +""" + +from skimage.draw import line_aa, circle_perimeter_aa + + +img = np.zeros((100, 100), dtype=np.double) + +# anti-aliased line +rr, cc, val = line_aa(12, 12, 20, 50) +img[rr, cc] = val + +# anti-aliased circle +rr, cc, val = circle_perimeter_aa(60, 40, 30) +img[rr, cc] = val + + +ax2.imshow(img, cmap=plt.cm.gray, interpolation='nearest') +ax2.set_title('Anti-aliasing') +ax2.axis('off') -plt.imshow(img) plt.show() diff --git a/doc/gh-pages.py b/doc/gh-pages.py index dd798c02..af992158 100644 --- a/doc/gh-pages.py +++ b/doc/gh-pages.py @@ -85,6 +85,10 @@ if __name__ == '__main__': for l in setup_lines: if l.startswith('VERSION'): tag = l.split("'")[1] + + # Rename to, e.g., 0.9.x + tag = '.'.join(tag.split('.')[:-1] + ['x']) + break if "dev" in tag: diff --git a/doc/release/contribs.py b/doc/release/contribs.py new file mode 100755 index 00000000..08a59fc9 --- /dev/null +++ b/doc/release/contribs.py @@ -0,0 +1,48 @@ +#!/usr/bin/env python +import subprocess +import sys +import string +import shlex + +if len(sys.argv) != 2: + print "Usage: ./contributors.py tag-of-previous-release" + sys.exit(-1) + +tag = sys.argv[1] + +def call(cmd): + return subprocess.check_output(shlex.split(cmd)).split('\n') + +tag_date = call("git show --format='%%ci' %s" % tag)[0] +print "Release %s was on %s" % (tag, tag_date) + +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 "\nIt contained the following %d merges:" % len(merges) +print +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" + +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.letters] + return author[-1] + +authors = sorted(set(authors), key=key) + +for a in authors: + print '-', a diff --git a/doc/release/contributors.sh b/doc/release/contributors.sh deleted file mode 100755 index 023928d3..00000000 --- a/doc/release/contributors.sh +++ /dev/null @@ -1,2 +0,0 @@ -git log $1..HEAD --format='- %aN' | sed 's/@/\-at\-/' | sed 's/<>//' | sort -u - diff --git a/doc/release/release_0.9.txt b/doc/release/release_0.9.txt new file mode 100644 index 00000000..a064ba2c --- /dev/null +++ b/doc/release/release_0.9.txt @@ -0,0 +1,130 @@ +Announcement: scikit-image 0.9.0 +================================ + +We're happy to announce the release of scikit-image v0.9.0! + +scikit-image is an image processing toolbox for SciPy that includes algorithms +for segmentation, geometric transformations, color space manipulation, +analysis, filtering, morphology, feature detection, and more. + +For more information, examples, and documentation, please visit our website: + + http://scikit-image.org + + +New Features +------------ + +`scikit-image` now runs without translation under both Python 2 and 3. + +In addition to several bug fixes, speed improvements and examples, the 204 pull +requests merged for this release include the following new features (PR number +in brackets): + +Segmentation: + +- 3D support in SLIC segmentation (#546) +- SLIC voxel spacing (#719) +- Generalized anisotropic spacing support for random_walker (#775) +- Yen threshold method (#686) + +Transforms and filters: + +- SART algorithm for tomography reconstruction (#584) +- Gabor filters (#371) +- Hough transform for ellipses (#597) +- Fast resampling of nD arrays (#511) +- Rotation axis center for Radon transforms with inverses. (#654) +- Reconstruction circle in inverse Radon transform (#567) +- Pixelwise image adjustment curves and methods (#505) + +Feature detection: + +- [experimental API] BRIEF feature descriptor (#591) +- [experimental API] Censure (STAR) Feature Detector (#668) +- Octagon structural element (#669) +- Add non rotation invariant uniform LBPs (#704) + +Color and noise: + +- Add deltaE color comparison and lab2lch conversion (#665) +- Isotropic denoising (#653) +- Generator to add various types of random noise to images (#625) +- Color deconvolution for immunohistochemical images (#441) +- Color label visualization (#485) + +Drawing and visualization: + +- Wu's anti-aliased circle, line, bezier curve (#709) +- Linked image viewers and docked plugins (#575) +- Rotated ellipse + bezier curve drawing (#510) +- PySide & PyQt4 compatibility in skimage-viewer (#551) + +Other: + +- Python 3 support without 2to3. (#620) +- 3D Marching Cubes (#469) +- Line, Circle, Ellipse total least squares fitting and RANSAC algorithm (#440) +- N-dimensional array padding (#577) +- Add a wrapper around `scipy.ndimage.gaussian_filter` with useful default behaviors. (#712) +- Predefined structuring elements for 3D morphology (#484) + + +API changes +----------- + +The following backward-incompatible API changes were made between 0.8 and 0.9: + +- No longer wrap ``imread`` output in an ``Image`` class +- Change default value of `sigma` parameter in ``skimage.segmentation.slic`` + to 0 +- ``hough_circle`` now returns a stack of arrays that are the same size as the + input image. Set the ``full_output`` flag to True for the old behavior. +- The following functions were deprecated over two releases: + `skimage.filter.denoise_tv_chambolle`, + `skimage.morphology.is_local_maximum`, `skimage.transform.hough`, + `skimage.transform.probabilistic_hough`,`skimage.transform.hough_peaks`. + Their functionality still exists, but under different names. + + +Contributors to this release +---------------------------- + +This release was made possible by the collaborative efforts of many +contributors, both new and old. They are listed in alphabetical order by +surname: + +- Ankit Agrawal +- K.-Michael Aye +- Chris Beaumont +- François Boulogne +- Luis Pedro Coelho +- Marianne Corvellec +- Olivier Debeir +- Ferdinand Deger +- Kemal Eren +- Jostein Bø Fløystad +- Christoph Gohlke +- Emmanuelle Gouillart +- Christian Horea +- Thouis (Ray) Jones +- Almar Klein +- Xavier Moles Lopez +- Alexis Mignon +- Juan Nunez-Iglesias +- Zachary Pincus +- Nicolas Pinto +- Davin Potts +- Malcolm Reynolds +- Umesh Sharma +- Johannes Schönberger +- Chintak Sheth +- Kirill Shklovsky +- Steven Silvester +- Matt Terry +- Riaan van den Dool +- Stéfan van der Walt +- Josh Warner +- Adam Wisniewski +- Yang Zetian +- Tony S Yu diff --git a/doc/source/_static/docversions.js b/doc/source/_static/docversions.js index ab333671..b4fd531f 100644 --- a/doc/source/_static/docversions.js +++ b/doc/source/_static/docversions.js @@ -1,4 +1,4 @@ -var versions = ['dev', '0.8.0', '0.7.0', '0.6', '0.5', '0.4', '0.3']; +var versions = ['dev', '0.9.x', '0.8.0', '0.7.0', '0.6', '0.5', '0.4', '0.3']; function insert_version_links() { for (i = 0; i < versions.length; i++){ diff --git a/doc/source/api_changes.txt b/doc/source/api_changes.txt index bf9ddb60..1f755635 100644 --- a/doc/source/api_changes.txt +++ b/doc/source/api_changes.txt @@ -3,6 +3,13 @@ Version 0.9 - No longer wrap ``imread`` output in an ``Image`` class - Change default value of `sigma` parameter in ``skimage.segmentation.slic`` to 0 +- ``hough_circle`` now returns a stack of arrays that are the same size as the + input image. Set the ``full_output`` flag to True for the old behavior. +- The following functions were deprecated over two releases: + `skimage.filter.denoise_tv_chambolle`, + `skimage.morphology.is_local_maximum`, `skimage.transform.hough`, + `skimage.transform.probabilistic_hough`,`skimage.transform.hough_peaks`. + Their functionality still exists, but under different names. Version 0.4 ----------- diff --git a/doc/tools/build_modref_templates.py b/doc/tools/build_modref_templates.py index c51c590a..b4db115d 100644 --- a/doc/tools/build_modref_templates.py +++ b/doc/tools/build_modref_templates.py @@ -35,7 +35,7 @@ if __name__ == '__main__': # are not (re)generated. This avoids automatic generation of documentation # for older or newer versions if such versions are installed on the system. - installed_version = V(module.version.version) + installed_version = V(module.__version__) setup_lines = open('../setup.py').readlines() version = 'vUndefined' diff --git a/doc/tools/plot_pr.py b/doc/tools/plot_pr.py index 08f4fc0e..5f9b4aa6 100644 --- a/doc/tools/plot_pr.py +++ b/doc/tools/plot_pr.py @@ -1,14 +1,15 @@ -import urllib import json -import copy +import urllib +import dateutil.parser from collections import OrderedDict +from datetime import datetime, timedelta +from dateutil.relativedelta import relativedelta +import numpy as np import matplotlib.pyplot as plt from matplotlib.ticker import FuncFormatter +from matplotlib.transforms import blended_transform_factory -import dateutil.parser -from dateutil.relativedelta import relativedelta -from datetime import datetime, timedelta cache = '_pr_cache.txt' @@ -22,16 +23,16 @@ cache = '_pr_cache.txt' releases = OrderedDict([ #('0.1', u'2009-10-07 13:52:19 +0200'), #('0.2', u'2009-11-12 14:48:45 +0200'), - ('0.3', u'2011-10-10 03:28:47 -0700'), + #('0.3', u'2011-10-10 03:28:47 -0700'), ('0.4', u'2011-12-03 14:31:32 -0800'), ('0.5', u'2012-02-26 21:00:51 -0800'), - ('0.6', u'2012-06-24 21:37:05 -0700')]) + ('0.6', u'2012-06-24 21:37:05 -0700'), + ('0.7', u'2012-09-29 18:08:49 -0700'), + ('0.8', u'2013-03-04 20:46:09 +0100')]) month_duration = 24 -for r in releases: - releases[r] = dateutil.parser.parse(releases[r]) def fetch_PRs(user='scikit-image', repo='scikit-image', state='open'): params = {'state': state, @@ -46,12 +47,12 @@ def fetch_PRs(user='scikit-image', repo='scikit-image', state='open'): 'repo': repo, 'params': urllib.urlencode(params)} - fetch_status = 'Fetching page %(page)d (state=%(state)s)' % params + \ - ' from %(user)s/%(repo)s...' % config + fetch_status = ('Fetching page %(page)d (state=%(state)s)' % params + + ' from %(user)s/%(repo)s...' % config) print(fetch_status) f = urllib.urlopen( - 'https://api.github.com/repos/%(user)s/%(repo)s/pulls?%(params)s' \ + 'https://api.github.com/repos/%(user)s/%(repo)s/pulls?%(params)s' % config ) @@ -67,6 +68,31 @@ def fetch_PRs(user='scikit-image', repo='scikit-image', state='open'): return data + +def seconds_from_epoch(dates): + seconds = [(dt - epoch).total_seconds() for dt in dates] + return seconds + + +def get_month_bins(dates): + now = datetime.now(tz=dates[0].tzinfo) + this_month = datetime(year=now.year, month=now.month, day=1, + tzinfo=dates[0].tzinfo) + + bins = [this_month - relativedelta(months=i) + for i in reversed(range(-1, month_duration))] + return seconds_from_epoch(bins) + + +def date_formatter(value, _): + dt = epoch + timedelta(seconds=value) + return dt.strftime('%Y/%m') + + +for r in releases: + releases[r] = dateutil.parser.parse(releases[r]) + + try: PRs = json.loads(open(cache, 'r').read()) print('Loaded PRs from cache...') @@ -87,47 +113,41 @@ dates = [dateutil.parser.parse(pr['created_at']) for pr in PRs] epoch = datetime(2009, 1, 1, tzinfo=dates[0].tzinfo) -def seconds_from_epoch(dates): - seconds = [(dt - epoch).total_seconds() for dt in dates] - return seconds - dates_f = seconds_from_epoch(dates) +bins = get_month_bins(dates) -def date_formatter(value, _): - dt = epoch + timedelta(seconds=value) - return dt.strftime('%Y/%m') +fig, ax = plt.subplots(figsize=(7, 5)) -plt.figure(figsize=(7, 5)) +n, bins, _ = ax.hist(dates_f, bins=bins, color='blue', alpha=0.6) -now = datetime.now(tz=dates[0].tzinfo) -this_month = datetime(year=now.year, month=now.month, day=1, - tzinfo=dates[0].tzinfo) - -bins = [this_month - relativedelta(months=i) \ - for i in reversed(range(-1, month_duration))] -bins = seconds_from_epoch(bins) -plt.hist(dates_f, bins=bins) - -ax = plt.gca() ax.xaxis.set_major_formatter(FuncFormatter(date_formatter)) -ax.set_xticks(bins[:-1]) +ax.set_xticks(bins[2:-1:3]) # Date label every 3 months. labels = ax.get_xticklabels() for l in labels: l.set_rotation(40) l.set_size(10) +mixed_transform = blended_transform_factory(ax.transData, ax.transAxes) for version, date in releases.items(): date = seconds_from_epoch([date])[0] - plt.axvline(date, color='r', label=version) + ax.axvline(date, color='black', linestyle=':', label=version) + ax.text(date, 1, version, color='r', va='bottom', ha='center', + transform=mixed_transform) -plt.title('Pull request activity').set_y(1.05) -plt.xlabel('Date') -plt.ylabel('PRs created') -plt.legend(loc=2, title='Release') -plt.subplots_adjust(top=0.875, bottom=0.225) +ax.set_title('Pull request activity').set_y(1.05) +ax.set_xlabel('Date') +ax.set_ylabel('PRs per month', color='blue') +fig.subplots_adjust(top=0.875, bottom=0.225) -plt.savefig('PRs.png') +cumulative = np.cumsum(n) +cumulative += len(dates) - cumulative[-1] + +ax2 = ax.twinx() +ax2.plot(bins[1:], cumulative, color='black', linewidth=2) +ax2.set_ylabel('Total PRs', color='black') + +fig.savefig('PRs.png') plt.show() diff --git a/setup.py b/setup.py index 1900443c..2420c33c 100755 --- a/setup.py +++ b/setup.py @@ -17,7 +17,7 @@ MAINTAINER_EMAIL = 'stefan@sun.ac.za' URL = 'http://scikit-image.org' LICENSE = 'Modified BSD' DOWNLOAD_URL = 'http://github.com/scikit-image/scikit-image' -VERSION = '0.8.1' +VERSION = '0.9.3' PYTHON_VERSION = (2, 5) DEPENDENCIES = { 'numpy': (1, 6), diff --git a/skimage/__init__.py b/skimage/__init__.py index 53d31eae..c73a79b7 100644 --- a/skimage/__init__.py +++ b/skimage/__init__.py @@ -91,6 +91,7 @@ test_verbose.__doc__ = test.__doc__ class _Log(Warning): pass + class _FakeLog(object): def __init__(self, name): """ diff --git a/skimage/color/colorlabel.py b/skimage/color/colorlabel.py index d6379037..8d7787aa 100644 --- a/skimage/color/colorlabel.py +++ b/skimage/color/colorlabel.py @@ -33,8 +33,32 @@ def _rgb_vector(color): """ if isinstance(color, six.string_types): color = color_dict[color] - # slice to handle RGBA colors - return np.array(color[:3]).reshape(1, 3) + # Slice to handle RGBA colors. + return np.array(color[:3]) + + +def _match_label_with_color(label, colors, bg_label, bg_color): + """Return `unique_labels` and `color_cycle` for label array and color list. + + Colors are cycled for normal labels, but the background color should only + be used for the background. + """ + # Temporarily set background color; it will be removed later. + if bg_color is None: + bg_color = (0, 0, 0) + bg_color = _rgb_vector([bg_color]) + + unique_labels = list(set(label.flat)) + # Ensure that the background label is in front to match call to `chain`. + if bg_label in unique_labels: + unique_labels.remove(bg_label) + unique_labels.insert(0, bg_label) + + # Modify labels and color cycle so background color is used only once. + color_cycle = itertools.cycle(colors) + color_cycle = itertools.chain(bg_color, color_cycle) + + return unique_labels, color_cycle def label2rgb(label, image=None, colors=None, alpha=0.3, @@ -66,7 +90,7 @@ def label2rgb(label, image=None, colors=None, alpha=0.3, colors = [_rgb_vector(c) for c in colors] if image is None: - colorized = np.zeros(label.shape + (3,), dtype=np.float64) + image = np.zeros(label.shape + (3,), dtype=np.float64) # Opacity doesn't make sense if no image exists. alpha = 1 else: @@ -77,20 +101,34 @@ def label2rgb(label, image=None, colors=None, alpha=0.3, warnings.warn("Negative intensities in `image` are not supported") image = img_as_float(rgb2gray(image)) - colorized = gray2rgb(image) * image_alpha + (1 - image_alpha) + image = gray2rgb(image) * image_alpha + (1 - image_alpha) - labels = list(set(label.flat)) - color_cycle = itertools.cycle(colors) + # Ensure that all labels are non-negative so we can index into + # `label_to_color` correctly. + offset = min(label.min(), bg_label) + if offset != 0: + label = label - offset # Make sure you don't modify the input array. + bg_label -= offset - if bg_label in labels: - labels.remove(bg_label) - if bg_color is not None: - labels.insert(0, bg_label) - bg_color = _rgb_vector(bg_color) - color_cycle = itertools.chain(bg_color, color_cycle) + new_type = np.min_scalar_type(label.max()) + if new_type == np.bool: + new_type = np.uint8 + label = label.astype(new_type) - for c, i in zip(color_cycle, labels): - mask = (label == i) - colorized[mask] = c * alpha + colorized[mask] * (1 - alpha) + unique_labels, color_cycle = _match_label_with_color(label, colors, + bg_label, bg_color) - return colorized + if len(unique_labels) == 0: + return image + + dense_labels = range(max(unique_labels) + 1) + label_to_color = np.array([c for i, c in zip(dense_labels, color_cycle)]) + + result = label_to_color[label] * alpha + image * (1 - alpha) + + # Remove background label if its color was not specified. + remove_background = bg_label in unique_labels and bg_color is None + if remove_background: + result[label == bg_label] = image[label == bg_label] + + return result diff --git a/skimage/color/tests/test_colorlabel.py b/skimage/color/tests/test_colorlabel.py index fa6ffcf3..dcfbe4ea 100644 --- a/skimage/color/tests/test_colorlabel.py +++ b/skimage/color/tests/test_colorlabel.py @@ -3,7 +3,8 @@ import itertools import numpy as np from numpy import testing from skimage.color.colorlabel import label2rgb -from numpy.testing import assert_array_almost_equal as assert_close +from numpy.testing import (assert_array_almost_equal as assert_close, + assert_array_equal) def test_shape_mismatch(): @@ -69,6 +70,26 @@ def test_bg_and_color_cycle(): assert_close(pixel, color) +def test_label_consistency(): + """Assert that the same labels map to the same colors.""" + label_1 = np.arange(5).reshape(1, -1) + label_2 = np.array([2, 4]) + colors = [(1, 0, 0), (0, 1, 0), (0, 0, 1), (1, 1, 0), (1, 0, 1)] + # Set alphas just in case the defaults change + rgb_1 = label2rgb(label_1, colors=colors) + rgb_2 = label2rgb(label_2, colors=colors) + for label_id in label_2.flat: + assert_close(rgb_1[label_1 == label_id], rgb_2[label_2 == label_id]) + +def test_leave_labels_alone(): + labels = np.array([-1, 0, 1]) + labels_saved = labels.copy() + + label2rgb(labels) + label2rgb(labels, bg_label=1) + assert_array_equal(labels, labels_saved) + + if __name__ == '__main__': testing.run_module_suite() diff --git a/skimage/data/__init__.py b/skimage/data/__init__.py index cea2a9f5..ecd261f7 100644 --- a/skimage/data/__init__.py +++ b/skimage/data/__init__.py @@ -200,4 +200,3 @@ def coffee(): """ return load("coffee.png") - diff --git a/skimage/draw/__init__.py b/skimage/draw/__init__.py index 4fb222d1..5f788ea7 100644 --- a/skimage/draw/__init__.py +++ b/skimage/draw/__init__.py @@ -1,9 +1,12 @@ from .draw import circle, ellipse, set_color -from ._draw import line, polygon, ellipse_perimeter, circle_perimeter, \ - bezier_segment from .draw3d import ellipsoid, ellipsoid_stats +from ._draw import (line, line_aa, polygon, ellipse_perimeter, + circle_perimeter, circle_perimeter_aa, + _bezier_segment, bezier_curve) __all__ = ['line', + 'line_aa', + 'bezier_curve', 'polygon', 'ellipse', 'ellipse_perimeter', @@ -11,4 +14,5 @@ __all__ = ['line', 'ellipsoid_stats', 'circle', 'circle_perimeter', + 'circle_perimeter_aa', 'set_color'] diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx index 132ba1d4..7300a2a1 100644 --- a/skimage/draw/_draw.pyx +++ b/skimage/draw/_draw.pyx @@ -6,7 +6,7 @@ import math import numpy as np cimport numpy as cnp -from libc.math cimport sqrt, sin, cos, floor +from libc.math cimport sqrt, sin, cos, floor, ceil from skimage._shared.geometry cimport point_in_polygon @@ -27,6 +27,10 @@ def line(Py_ssize_t y, Py_ssize_t x, Py_ssize_t y2, Py_ssize_t x2): May be used to directly index into an array, e.g. ``img[rr, cc] = 1``. + Notes + ----- + Anti-aliased line generator is available with `line_aa`. + Examples -------- >>> from skimage.draw import line @@ -67,8 +71,8 @@ def line(Py_ssize_t y, Py_ssize_t x, Py_ssize_t y2, Py_ssize_t x2): sx, sy = sy, sx d = (2 * dy) - dx - cdef Py_ssize_t[:] rr = np.zeros(int(dx) + 1, dtype=np.intp) - cdef Py_ssize_t[:] cc = np.zeros(int(dx) + 1, dtype=np.intp) + cdef Py_ssize_t[::1] rr = np.zeros(int(dx) + 1, dtype=np.intp) + cdef Py_ssize_t[::1] cc = np.zeros(int(dx) + 1, dtype=np.intp) for i in range(dx): if steep: @@ -89,6 +93,99 @@ def line(Py_ssize_t y, Py_ssize_t x, Py_ssize_t y2, Py_ssize_t x2): return np.asarray(rr), np.asarray(cc) +def line_aa(Py_ssize_t y1, Py_ssize_t x1, Py_ssize_t y2, Py_ssize_t x2): + """Generate anti-aliased line pixel coordinates. + + Parameters + ---------- + y1, x1 : int + Starting position (row, column). + y2, x2 : int + End position (row, column). + + Returns + ------- + rr, cc, val : (N,) ndarray (int, int, float) + Indices of pixels (`rr`, `cc`) and intensity values (`val`). + ``img[rr, cc] = val``. + + References + ---------- + .. [1] A Rasterizing Algorithm for Drawing Curves, A. Zingl, 2012 + http://members.chello.at/easyfilter/Bresenham.pdf + + Examples + -------- + >>> from skimage.draw import line_aa + >>> img = np.zeros((10, 10), dtype=np.uint8) + >>> rr, cc, val = line_aa(1, 1, 8, 8) + >>> img[rr, cc] = val * 255 + >>> img + array([[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 255, 56, 0, 0, 0, 0, 0, 0, 0], + [ 0, 56, 255, 56, 0, 0, 0, 0, 0, 0], + [ 0, 0, 56, 255, 56, 0, 0, 0, 0, 0], + [ 0, 0, 0, 56, 255, 56, 0, 0, 0, 0], + [ 0, 0, 0, 0, 56, 255, 56, 0, 0, 0], + [ 0, 0, 0, 0, 0, 56, 255, 56, 0, 0], + [ 0, 0, 0, 0, 0, 0, 56, 255, 56, 0], + [ 0, 0, 0, 0, 0, 0, 0, 56, 255, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) + """ + cdef list rr = list() + cdef list cc = list() + cdef list val = list() + + cdef int dx = abs(x1 - x2) + cdef int dy = abs(y1 - y2) + cdef int err = dx - dy + cdef int x, y, e, ed, sign_x, sign_y + + if x1 < x2: + sign_x = 1 + else: + sign_x = -1 + + if y1 < y2: + sign_y = 1 + else: + sign_y = -1 + + if dx + dy == 0: + ed = 1 + else: + ed = (sqrt(dx*dx + dy*dy)) + + x, y = x1, y1 + while True: + cc.append(x) + rr.append(y) + val.append(1. * abs(err - dx + dy) / (ed)) + e = err + if 2 * e >= -dx: + if x == x2: + break + if e + dy < ed: + cc.append(x) + rr.append(y + sign_y) + val.append(1. * abs(e + dy) / (ed)) + err -= dy + x += sign_x + if 2 * e <= dy: + if y == y2: + break + if dx - e < ed: + cc.append(x) + rr.append(y) + val.append(abs(dx - e) / (ed)) + err += dx + y += sign_y + + return (np.array(rr, dtype=np.intp), + np.array(cc, dtype=np.intp), + 1. - np.array(val, dtype=np.float)) + + def polygon(y, x, shape=None): """Generate coordinates of pixels within polygon. @@ -134,9 +231,9 @@ def polygon(y, x, shape=None): cdef Py_ssize_t nr_verts = x.shape[0] cdef Py_ssize_t minr = int(max(0, y.min())) - cdef Py_ssize_t maxr = int(math.ceil(y.max())) + cdef Py_ssize_t maxr = int(ceil(y.max())) cdef Py_ssize_t minc = int(max(0, x.min())) - cdef Py_ssize_t maxc = int(math.ceil(x.max())) + cdef Py_ssize_t maxc = int(ceil(x.max())) # make sure output coordinates do not exceed image size if shape is not None: @@ -182,6 +279,7 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius, Returns ------- rr, cc : (N,) ndarray of int + Bresenham and Andres' method: Indices of pixels that belong to the circle perimeter. May be used to directly index into an array, e.g. ``img[rr, cc] = 1``. @@ -192,13 +290,14 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius, circles create a disc whereas Bresenham can make holes. There is also less distortions when Andres circles are rotated. Bresenham method is also known as midpoint circle algorithm. + Anti-aliased circle generator is available with `circle_perimeter_aa`. References ---------- .. [1] J.E. Bresenham, "Algorithm for computer control of a digital - plotter", 4 (1965) 25-30. - .. [2] E. Andres, "Discrete circles, rings and spheres", - 18 (1994) 695-706. + plotter", IBM Systems journal, 4 (1965) 25-30. + .. [2] E. Andres, "Discrete circles, rings and spheres", Computers & + Graphics, 18 (1994) 695-706. Examples -------- @@ -226,6 +325,10 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius, cdef Py_ssize_t x = 0 cdef Py_ssize_t y = radius cdef Py_ssize_t d = 0 + + cdef double dceil = 0 + cdef double dceil_prev = 0 + cdef char cmethod if method == 'bresenham': d = 3 - 2 * radius @@ -258,8 +361,84 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius, d = d + 2 * (y - x - 1) y = y - 1 x = x + 1 + return (np.array(rr, dtype=np.intp) + cy, + np.array(cc, dtype=np.intp) + cx) - return np.array(rr, dtype=np.intp) + cy, np.array(cc, dtype=np.intp) + cx + +def circle_perimeter_aa(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius): + """Generate anti-aliased circle perimeter coordinates. + + Parameters + ---------- + cy, cx : int + Centre coordinate of circle. + radius: int + Radius of circle. + + Returns + ------- + rr, cc, val : (N,) ndarray (int, int, float) + Indices of pixels (`rr`, `cc`) and intensity values (`val`). + ``img[rr, cc] = val``. + + Notes + ----- + Wu's method draws anti-aliased circle. This implementation doesn't use + lookup table optimization. + + References + ---------- + .. [1] X. Wu, "An efficient antialiasing technique", In ACM SIGGRAPH + Computer Graphics, 25 (1991) 143-152. + + Examples + -------- + >>> from skimage.draw import circle_perimeter_aa + >>> img = np.zeros((10, 10), dtype=np.uint8) + >>> rr, cc, val = circle_perimeter_aa(4, 4, 3) + >>> img[rr, cc] = val * 255 + >>> img + array([[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 60, 211, 255, 211, 60, 0, 0, 0], + [ 0, 60, 194, 43, 0, 43, 194, 60, 0, 0], + [ 0, 211, 43, 0, 0, 0, 43, 211, 0, 0], + [ 0, 255, 0, 0, 0, 0, 0, 255, 0, 0], + [ 0, 211, 43, 0, 0, 0, 43, 211, 0, 0], + [ 0, 60, 194, 43, 0, 43, 194, 60, 0, 0], + [ 0, 0, 60, 211, 255, 211, 60, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) + """ + + cdef Py_ssize_t x = 0 + cdef Py_ssize_t y = radius + cdef Py_ssize_t d = 0 + + cdef double dceil = 0 + cdef double dceil_prev = 0 + + cdef list rr = [y, x, y, x, -y, -x, -y, -x] + cdef list cc = [x, y, -x, -y, x, y, -x, -y] + cdef list val = [1] * 8 + + while y > x + 1: + x += 1 + dceil = sqrt(radius**2 - x**2) + dceil = ceil(dceil) - dceil + if dceil < dceil_prev: + y -= 1 + rr.extend([y, y - 1, x, x, y, y - 1, x, x]) + cc.extend([x, x, y, y - 1, -x, -x, -y, 1 - y]) + + rr.extend([-y, 1 - y, -x, -x, -y, 1 - y, -x, -x]) + cc.extend([x, x, y, y - 1, -x, -x, -y, 1 - y]) + + val.extend([1 - dceil, dceil] * 8) + dceil_prev = dceil + + return (np.array(rr, dtype=np.intp) + cy, + np.array(cc, dtype=np.intp) + cx, + np.array(val, dtype=np.float)) def ellipse_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t yradius, @@ -270,9 +449,9 @@ def ellipse_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t yradius, ---------- cy, cx : int Centre coordinate of ellipse. - yradius, xradius: int + yradius, xradius : int Minor and major semi-axes. ``(x/xradius)**2 + (y/yradius)**2 = 1``. - orientation: double, optional (default 0) + orientation : double, optional (default 0) Major axis orientation in clockwise direction as radians. Returns @@ -382,38 +561,38 @@ def ellipse_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t yradius, iyd = int(floor(ya * w + 0.5)) # Draw the 4 quadrants - rr, cc = bezier_segment(iy0 + iyd, ix0, iy0, ix0, iy0, ix0 + ixd, 1-w) + rr, cc = _bezier_segment(iy0 + iyd, ix0, iy0, ix0, iy0, ix0 + ixd, 1-w) py.extend(rr) px.extend(cc) - rr, cc = bezier_segment(iy0 + iyd, ix0, iy1, ix0, iy1, ix1 - ixd, w) + rr, cc = _bezier_segment(iy0 + iyd, ix0, iy1, ix0, iy1, ix1 - ixd, w) py.extend(rr) px.extend(cc) - rr, cc = bezier_segment(iy1 - iyd, ix1, iy1, ix1, iy1, ix1 - ixd, 1-w) + rr, cc = _bezier_segment(iy1 - iyd, ix1, iy1, ix1, iy1, ix1 - ixd, 1-w) py.extend(rr) px.extend(cc) - rr, cc = bezier_segment(iy1 - iyd, ix1, iy0, ix1, iy0, ix0 + ixd, w) + rr, cc = _bezier_segment(iy1 - iyd, ix1, iy0, ix1, iy0, ix0 + ixd, w) py.extend(rr) px.extend(cc) return np.array(py, dtype=np.intp), np.array(px, dtype=np.intp) -def bezier_segment(Py_ssize_t y0, Py_ssize_t x0, - Py_ssize_t y1, Py_ssize_t x1, - Py_ssize_t y2, Py_ssize_t x2, - double weight): +def _bezier_segment(Py_ssize_t y0, Py_ssize_t x0, + Py_ssize_t y1, Py_ssize_t x1, + Py_ssize_t y2, Py_ssize_t x2, + double weight): """Generate Bezier segment coordinates. Parameters ---------- y0, x0 : int - Coordinates of the first point + Coordinates of the first control point. y1, x1 : int - Coordinates of the middle point + Coordinates of the middle control point. y2, x2 : int - Coordinates of the last point + Coordinates of the last control point. weight : double - Middle point weight, it describes the line tension. + Middle control point weight, it describes the line tension. Returns ------- @@ -425,7 +604,7 @@ def bezier_segment(Py_ssize_t y0, Py_ssize_t x0, Notes ----- The algorithm is the rational quadratic algorithm presented in - reference [1]. + reference [1]_. References ---------- @@ -492,8 +671,8 @@ def bezier_segment(Py_ssize_t y0, Py_ssize_t x0, sy = floor((y0 + 2 * weight * y1 + y2) * xy * 0.5 + 0.5) dx = floor((weight * x1 + x0) * xy + 0.5) dy = floor((y1 * weight + y0) * xy + 0.5) - return bezier_segment(y0, x0, (dy), (dx), - (sy), (sx), cur) + return _bezier_segment(y0, x0, (dy), (dx), + (sy), (sx), cur) err = dx + dy - xy while dy <= xy and dx >= xy: @@ -526,29 +705,132 @@ def bezier_segment(Py_ssize_t y0, Py_ssize_t x0, return np.array(py, dtype=np.intp), np.array(px, dtype=np.intp) -def set_color(img, coords, color): - """Set pixel color in the image at the given coordinates. - - Coordinates that exceed the shape of the image will be ignored. +def bezier_curve(Py_ssize_t y0, Py_ssize_t x0, + Py_ssize_t y1, Py_ssize_t x1, + Py_ssize_t y2, Py_ssize_t x2, + double weight): + """Generate Bezier curve coordinates. Parameters ---------- - img : (M, N, D) ndarray - Image - coords : ((P,) ndarray, (P,) ndarray) - Coordinates of pixels to be colored. - color : (D,) ndarray - Color to be assigned to coordinates in the image. + y0, x0 : int + Coordinates of the first control point. + y1, x1 : int + Coordinates of the middle control point. + y2, x2 : int + Coordinates of the last control point. + weight : double + Middle control point weight, it describes the line tension. Returns ------- - img : (M, N, D) ndarray - The updated image. + rr, cc : (N,) ndarray of int + Indices of pixels that belong to the Bezier curve. + May be used to directly index into an array, e.g. + ``img[rr, cc] = 1``. + Notes + ----- + The algorithm is the rational quadratic algorithm presented in + reference [1]_. + + References + ---------- + .. [1] A Rasterizing Algorithm for Drawing Curves, A. Zingl, 2012 + http://members.chello.at/easyfilter/Bresenham.pdf + + Examples + -------- + >>> import numpy as np + >>> from skimage.draw import bezier_curve + >>> img = np.zeros((10, 10), dtype=np.uint8) + >>> rr, cc = bezier_curve(1, 5, 5, -2, 8, 8, 2) + >>> img[rr, cc] = 1 + >>> img + array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) """ + # Pixels + cdef list px = list() + cdef list py = list() - rr, cc = coords - rr_inside = np.logical_and(rr >= 0, rr < img.shape[0]) - cc_inside = np.logical_and(cc >= 0, cc < img.shape[1]) - inside = np.logical_and(rr_inside, cc_inside) - img[rr[inside], cc[inside]] = color + cdef int x, y + cdef double xx, yy, ww, t, q + x = x0 - 2 * x1 + x2 + y = y0 - 2 * y1 + y2 + + xx = x0 - x1 + yy = y0 - y1 + + if xx * (x2 - x1) > 0: + if yy * (y2 - y1): + if abs(xx * y) > abs(yy * x): + x0 = x2 + x2 = (xx + x1) + y0 = y2 + y2 = (yy + y1) + if (x0 == x2) or (weight == 1.): + t = (x0 - x1) / x + else: + q = sqrt(4. * weight * weight * (x0 - x1) * (x2 - x1) + (x2 - x0) * floor(x2 - x0)) + if (x1 < x0): + q = -q + t = (2. * weight * (x0 - x1) - x0 + x2 + q) / (2. * (1. - weight) * (x2 - x0)) + + q = 1. / (2. * t * (1. - t) * (weight - 1.) + 1.0) + xx = (t * t * (x0 - 2. * weight * x1 + x2) + 2. * t * (weight * x1 - x0) + x0) * q + yy = (t * t * (y0 - 2. * weight * y1 + y2) + 2. * t * (weight * y1 - y0) + y0) * q + ww = t * (weight - 1.) + 1. + ww *= ww * q + weight = ((1. - t) * (weight - 1.) + 1.) * sqrt(q) + x = (xx + 0.5) + y = (yy + 0.5) + yy = (xx - x0) * (y1 - y0) / (x1 - x0) + y0 + + rr, cc = _bezier_segment(y0, x0, (yy + 0.5), x, y, x, ww) + px.extend(rr) + py.extend(cc) + + yy = (xx - x2) * (y1 - y2) / (x1 - x2) + y2 + y1 = (yy + 0.5) + x0 = x1 = x + y0 = y + if (y0 - y1) * floor(y2 - y1) > 0: + if (y0 == y2) or (weight == 1): + t = (y0 - y1) / (y0 - 2. * y1 + y2) + else: + q = sqrt(4. * weight * weight * (y0 - y1) * (y2 - y1) + (y2 - y0) * floor(y2 - y0)) + if y1 < y0: + q = -q + t = (2. * weight * (y0 - y1) - y0 + y2 + q) / (2. * (1. - weight) * (y2 - y0)) + q = 1. / (2. * t * (1. - t) * (weight - 1.) + 1.) + xx = (t * t * (x0 - 2. * weight * x1 + x2) + 2. * t * (weight * x1 - x0) + x0) * q + yy = (t * t * (y0 - 2. * weight * y1 + y2) + 2. * t * (weight * y1 - y0) + y0) * q + ww = t * (weight - 1.) + 1. + ww *= ww * q + weight = ((1. - t) * (weight - 1.) + 1.) * sqrt(q) + x = (xx + 0.5) + y = (yy + 0.5) + xx = (x1 - x0) * (yy - y0) / (y1 - y0) + x0 + + rr, cc = _bezier_segment(y0, x0, y, (xx + 0.5), y, x, ww) + px.extend(rr) + py.extend(cc) + + xx = (x1 - x2) * (yy - y2) / (y1 - y2) + x2 + x1 = (xx + 0.5) + x0 = x + y0 = y1 = y + + rr, cc = _bezier_segment(y0, x0, y1, x1, y2, x2, weight * weight) + px.extend(rr) + py.extend(cc) + return np.array(px, dtype=np.intp), np.array(py, dtype=np.intp) diff --git a/skimage/draw/draw.py b/skimage/draw/draw.py index d01bc2b0..cbf3ced2 100644 --- a/skimage/draw/draw.py +++ b/skimage/draw/draw.py @@ -52,7 +52,7 @@ def ellipse(cy, cx, yradius, xradius, shape=None): dc = 1 / float(xradius) r, c = np.ogrid[-1:1:dr, -1:1:dc] - rr, cc = np.nonzero(r ** 2 + c ** 2 < 1) + rr, cc = np.nonzero(r ** 2 + c ** 2 < 1) rr.flags.writeable = True cc.flags.writeable = True diff --git a/skimage/draw/draw3d.py b/skimage/draw/draw3d.py index 0b6fcb2d..db2f0d39 100644 --- a/skimage/draw/draw3d.py +++ b/skimage/draw/draw3d.py @@ -3,10 +3,10 @@ import numpy as np from scipy.special import (ellipkinc as ellip_F, ellipeinc as ellip_E) -def ellipsoid(a, b, c, sampling=(1., 1., 1.), levelset=False): +def ellipsoid(a, b, c, spacing=(1., 1., 1.), levelset=False): """ Generates ellipsoid with semimajor axes aligned with grid dimensions - on grid with specified `sampling`. + on grid with specified `spacing`. Parameters ---------- @@ -16,8 +16,8 @@ def ellipsoid(a, b, c, sampling=(1., 1., 1.), levelset=False): Length of semimajor axis aligned with y-axis. c : float Length of semimajor axis aligned with z-axis. - sampling : tuple of floats, length 3 - Sampling in (x, y, z) spatial dimensions. + spacing : tuple of floats, length 3 + Spacing in (x, y, z) spatial dimensions. levelset : bool If True, returns the level set for this ellipsoid (signed level set about zero, with positive denoting interior) as np.float64. @@ -26,7 +26,7 @@ def ellipsoid(a, b, c, sampling=(1., 1., 1.), levelset=False): Returns ------- ellip : (N, M, P) array - Ellipsoid centered in a correctly sized array for given `sampling`. + Ellipsoid centered in a correctly sized array for given `spacing`. Boolean dtype unless `levelset=True`, in which case a float array is returned with the level set above 0.0 representing the ellipsoid. @@ -34,7 +34,7 @@ def ellipsoid(a, b, c, sampling=(1., 1., 1.), levelset=False): if (a <= 0) or (b <= 0) or (c <= 0): raise ValueError('Parameters a, b, and c must all be > 0') - offset = np.r_[1, 1, 1] * np.r_[sampling] + offset = np.r_[1, 1, 1] * np.r_[spacing] # Calculate limits, and ensure output volume is odd & symmetric low = np.ceil((- np.r_[a, b, c] - offset)) @@ -43,14 +43,14 @@ def ellipsoid(a, b, c, sampling=(1., 1., 1.), levelset=False): for dim in range(3): if (high[dim] - low[dim]) % 2 == 0: low[dim] -= 1 - num = np.arange(low[dim], high[dim], sampling[dim]) + num = np.arange(low[dim], high[dim], spacing[dim]) if 0 not in num: low[dim] -= np.max(num[num < 0]) # Generate (anisotropic) spatial grid - x, y, z = np.mgrid[low[0]:high[0]:sampling[0], - low[1]:high[1]:sampling[1], - low[2]:high[2]:sampling[2]] + x, y, z = np.mgrid[low[0]:high[0]:spacing[0], + low[1]:high[1]:spacing[1], + low[2]:high[2]:spacing[2]] if not levelset: arr = ((x / float(a)) ** 2 + @@ -64,10 +64,10 @@ def ellipsoid(a, b, c, sampling=(1., 1., 1.), levelset=False): return arr -def ellipsoid_stats(a, b, c, sampling=(1., 1., 1.)): +def ellipsoid_stats(a, b, c): """ Calculates analytical surface area and volume for ellipsoid with - semimajor axes aligned with grid dimensions of specified `sampling`. + semimajor axes aligned with grid dimensions of specified `spacing`. Parameters ---------- @@ -77,8 +77,6 @@ def ellipsoid_stats(a, b, c, sampling=(1., 1., 1.)): Length of semimajor axis aligned with y-axis. c : float Length of semimajor axis aligned with z-axis. - sampling : tuple of floats, length 3 - Sampling in (x, y, z) spatial dimensions. Returns ------- diff --git a/skimage/draw/tests/test_draw.py b/skimage/draw/tests/test_draw.py index f15fba01..2d739f0f 100644 --- a/skimage/draw/tests/test_draw.py +++ b/skimage/draw/tests/test_draw.py @@ -1,8 +1,22 @@ -from numpy.testing import assert_array_equal +from numpy.testing import assert_array_equal, assert_equal import numpy as np -from skimage.draw import line, polygon, circle, circle_perimeter, \ - ellipse, ellipse_perimeter, bezier_segment +from skimage.draw import (set_color, line, line_aa, polygon, + circle, circle_perimeter, circle_perimeter_aa, + ellipse, ellipse_perimeter, + _bezier_segment, bezier_curve) + + +def test_set_color(): + img = np.zeros((10, 10)) + + rr, cc = line(0, 0, 0, 30) + set_color(img, (rr, cc), 1) + + img_ = np.zeros((10, 10)) + img_[0, :] = 1 + + assert_array_equal(img, img_) def test_line_horizontal(): @@ -52,6 +66,43 @@ def test_line_diag(): assert_array_equal(img, img_) +def test_line_aa_horizontal(): + img = np.zeros((10, 10)) + + rr, cc, val = line_aa(0, 0, 0, 9) + img[rr, cc] = val + + img_ = np.zeros((10, 10)) + img_[0, :] = 1 + + assert_array_equal(img, img_) + + +def test_line_aa_vertical(): + img = np.zeros((10, 10)) + + rr, cc, val = line_aa(0, 0, 9, 0) + img[rr, cc] = val + + img_ = np.zeros((10, 10)) + img_[:, 0] = 1 + + assert_array_equal(img, img_) + + +def test_line_aa_diagonal(): + img = np.zeros((10, 10)) + + rr, cc, val = line_aa(0, 0, 9, 6) + img[rr, cc] = 1 + + # Check that each pixel belonging to line, + # also belongs to line_aa + r, c = line(0, 0, 9, 6) + for x, y in zip(r, c): + assert_equal(img[r, c], 1) + + def test_polygon_rectangle(): img = np.zeros((10, 10), 'uint8') poly = np.array(((1, 1), (4, 1), (4, 4), (1, 4), (1, 1))) @@ -215,6 +266,38 @@ def test_circle_perimeter_andres(): assert_array_equal(img, img_) +def test_circle_perimeter_aa(): + img = np.zeros((15, 15), 'uint8') + rr, cc, val = circle_perimeter_aa(7, 7, 0) + img[rr, cc] = 1 + assert(np.sum(img) == 1) + assert(img[7][7] == 1) + + img = np.zeros((17, 17), 'uint8') + rr, cc, val = circle_perimeter_aa(8, 8, 7) + img[rr, cc] = val * 255 + img_ = np.array( + [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 82, 180, 236, 255, 236, 180, 82, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 189, 172, 74, 18, 0, 18, 74, 172, 189, 0, 0, 0, 0], + [ 0, 0, 0, 229, 25, 0, 0, 0, 0, 0, 0, 0, 25, 229, 0, 0, 0], + [ 0, 0, 189, 25, 0, 0, 0, 0, 0, 0, 0, 0, 0, 25, 189, 0, 0], + [ 0, 82, 172, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 172, 82, 0], + [ 0, 180, 74, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 74, 180, 0], + [ 0, 236, 18, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 18, 236, 0], + [ 0, 255, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 255, 0], + [ 0, 236, 18, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 18, 236, 0], + [ 0, 180, 74, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 74, 180, 0], + [ 0, 82, 172, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 172, 82, 0], + [ 0, 0, 189, 25, 0, 0, 0, 0, 0, 0, 0, 0, 0, 25, 189, 0, 0], + [ 0, 0, 0, 229, 25, 0, 0, 0, 0, 0, 0, 0, 25, 229, 0, 0, 0], + [ 0, 0, 0, 0, 189, 172, 74, 18, 0, 18, 74, 172, 189, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 82, 180, 236, 255, 236, 180, 82, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]] + ) + assert_array_equal(img, img_) + + def test_ellipse(): img = np.zeros((15, 15), 'uint8') @@ -360,18 +443,21 @@ def test_bezier_segment_straight(): y1 = 50 x2 = 150 y2 = 150 - rr, cc = bezier_segment(x0, y0, x1, y1, x2, y2, 0) - image [rr, cc] = 1 + rr, cc = _bezier_segment(x0, y0, x1, y1, x2, y2, 0) + image[rr, cc] = 1 image2 = np.zeros((200, 200), dtype=int) rr, cc = line(x0, y0, x2, y2) - image2 [rr, cc] = 1 + image2[rr, cc] = 1 assert_array_equal(image, image2) def test_bezier_segment_curved(): img = np.zeros((25, 25), 'uint8') - rr, cc = bezier_segment(20, 20, 20, 2, 2, 2, 1) + x1, y1 = 20, 20 + x2, y2 = 20, 2 + x3, y3 = 2, 2 + rr, cc = _bezier_segment(x1, y1, x2, y2, x3, y3, 1) img[rr, cc] = 1 img_ = np.array( [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], @@ -400,9 +486,101 @@ def test_bezier_segment_curved(): [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]] ) + assert_equal(img[x1, y1], 1) + assert_equal(img[x3, y3], 1) assert_array_equal(img, img_) +def test_bezier_curve_straight(): + image = np.zeros((200, 200), dtype=int) + x0 = 50 + y0 = 50 + x1 = 150 + y1 = 50 + x2 = 150 + y2 = 150 + rr, cc = bezier_curve(x0, y0, x1, y1, x2, y2, 0) + image [rr, cc] = 1 + + image2 = np.zeros((200, 200), dtype=int) + rr, cc = line(x0, y0, x2, y2) + image2 [rr, cc] = 1 + assert_array_equal(image, image2) + + +def test_bezier_curved_weight_eq_1(): + img = np.zeros((23, 8), 'uint8') + x1, y1 = (1, 1) + x2, y2 = (11, 11) + x3, y3 = (21, 1) + rr, cc = bezier_curve(x1, y1, x2, y2, x3, y3, 1) + img[rr, cc] = 1 + assert_equal(img[x1, y1], 1) + assert_equal(img[x3, y3], 1) + img_ = np.array( + [[0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0]] + ) + assert_equal(img, img_) + + +def test_bezier_curved_weight_neq_1(): + img = np.zeros((23, 10), 'uint8') + x1, y1 = (1, 1) + x2, y2 = (11, 11) + x3, y3 = (21, 1) + rr, cc = bezier_curve(x1, y1, x2, y2, x3, y3, 2) + img[rr, cc] = 1 + assert_equal(img[x1, y1], 1) + assert_equal(img[x3, y3], 1) + img_ = np.array( + [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]] + ) + assert_equal(img, img_) + if __name__ == "__main__": from numpy.testing import run_module_suite run_module_suite() diff --git a/skimage/draw/tests/test_draw3d.py b/skimage/draw/tests/test_draw3d.py index 59a7f6b3..2e1198eb 100644 --- a/skimage/draw/tests/test_draw3d.py +++ b/skimage/draw/tests/test_draw3d.py @@ -6,7 +6,7 @@ from skimage.draw import ellipsoid, ellipsoid_stats def test_ellipsoid_bool(): test = ellipsoid(2, 2, 2)[1:-1, 1:-1, 1:-1] - test_anisotropic = ellipsoid(2, 2, 4, sampling=(1., 1., 2.)) + test_anisotropic = ellipsoid(2, 2, 4, spacing=(1., 1., 2.)) test_anisotropic = test_anisotropic[1:-1, 1:-1, 1:-1] expected = np.array([[[0, 0, 0, 0, 0], @@ -45,7 +45,7 @@ def test_ellipsoid_bool(): def test_ellipsoid_levelset(): test = ellipsoid(2, 2, 2, levelset=True)[1:-1, 1:-1, 1:-1] - test_anisotropic = ellipsoid(2, 2, 4, sampling=(1., 1., 2.), + test_anisotropic = ellipsoid(2, 2, 4, spacing=(1., 1., 2.), levelset=True) test_anisotropic = test_anisotropic[1:-1, 1:-1, 1:-1] diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py index 9c50ab7d..fd5d53dd 100644 --- a/skimage/exposure/exposure.py +++ b/skimage/exposure/exposure.py @@ -287,7 +287,7 @@ def adjust_log(image, gain=1, inv=False): inv : float If True, it performs inverse logarithmic correction, else correction will be logarithmic. Defaults to False. - + Returns ------- out : ndarray diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index bde81e59..4a6518d6 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -7,9 +7,7 @@ from .corner import (corner_kitchen_rosenfeld, corner_harris, corner_peaks) from .corner_cy import corner_moravec from .template import match_template -from ._brief import brief, match_keypoints_brief -from .util import pairwise_hamming_distance -from .censure import keypoints_censure + __all__ = ['daisy', 'hog', @@ -24,8 +22,4 @@ __all__ = ['daisy', 'corner_subpix', 'corner_peaks', 'corner_moravec', - 'match_template', - 'brief', - 'pairwise_hamming_distance', - 'match_keypoints_brief', - 'keypoints_censure'] + 'match_template'] diff --git a/skimage/feature/_brief.py b/skimage/feature/_brief.py index 27a8d085..ecc2ec11 100644 --- a/skimage/feature/_brief.py +++ b/skimage/feature/_brief.py @@ -9,7 +9,9 @@ from ._brief_cy import _brief_loop def brief(image, keypoints, descriptor_size=256, mode='normal', patch_size=49, sample_seed=1, variance=2): - """Extract BRIEF Descriptor about given keypoints for a given image. + """**Experimental function**. + + Extract BRIEF Descriptor about given keypoints for a given image. Parameters ---------- @@ -178,7 +180,9 @@ def brief(image, keypoints, descriptor_size=256, mode='normal', patch_size=49, def match_keypoints_brief(keypoints1, descriptors1, keypoints2, descriptors2, threshold=0.15): - """Match keypoints described using BRIEF descriptors in one image to + """**Experimental function**. + + Match keypoints described using BRIEF descriptors in one image to those in second image. Parameters diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx index 6caa7ea3..ec83fa65 100644 --- a/skimage/feature/_texture.pyx +++ b/skimage/feature/_texture.pyx @@ -106,17 +106,17 @@ def _local_binary_pattern(double[:, ::1] image, """ # texture weights - cdef int[:] weights = 2 ** np.arange(P, dtype=np.int32) + cdef int[::1] weights = 2 ** np.arange(P, dtype=np.int32) # local position of texture elements rr = - R * np.sin(2 * np.pi * np.arange(P, dtype=np.double) / P) cc = R * np.cos(2 * np.pi * np.arange(P, dtype=np.double) / P) - cdef double[:] rp = np.round(rr, 5) - cdef double[:] cp = np.round(cc, 5) + cdef double[::1] rp = np.round(rr, 5) + cdef double[::1] cp = np.round(cc, 5) # pre-allocate arrays for computation - cdef double[:] texture = np.zeros(P, dtype=np.double) - cdef char[:] signed_texture = np.zeros(P, dtype=np.int8) - cdef int[:] rotation_chain = np.zeros(P, dtype=np.int32) + cdef double[::1] texture = np.zeros(P, dtype=np.double) + cdef char[::1] signed_texture = np.zeros(P, dtype=np.int8) + cdef int[::1] rotation_chain = np.zeros(P, dtype=np.int32) output_shape = (image.shape[0], image.shape[1]) cdef double[:, ::1] output = np.zeros(output_shape, dtype=np.double) @@ -162,21 +162,21 @@ def _local_binary_pattern(double[:, ::1] image, # n_ones=2: 0011, 1001, 1100, 0110 # n_ones=3: 0111, 1011, 1101, 1110 # n_ones=4: 1111 - # + # # For a pattern of size P there are 2 constant patterns # corresponding to n_ones=0 and n_ones=P. For each other # value of `n_ones` , i.e n_ones=[1..P-1], there are P # possible patterns which are related to each other through # circular permutations. The total number of uniform - # patterns is thus (2 + P * (P - 1)). + # patterns is thus (2 + P * (P - 1)). # Given any pattern (uniform or not) we must be able to - # associate a unique code: + # associate a unique code: # 1. Constant patterns patterns (with n_ones=0 and # n_ones=P) and non uniform patterns are given fixed # code values. # 2. Other uniform patterns are indexed considering the # value of n_ones, and an index called 'rot_index' - # reprenting the number of circular right shifts + # reprenting the number of circular right shifts # required to obtain the pattern starting from a # reference position (corresponding to all zeros stacked # on the right). This number of rotations (or circular @@ -215,7 +215,7 @@ def _local_binary_pattern(double[:, ::1] image, lbp += signed_texture[i] else: lbp = P + 1 - + if method == 'V': var = np.var(texture) if var != 0: diff --git a/skimage/feature/censure.py b/skimage/feature/censure.py index 00be16a8..4bb7fdda 100644 --- a/skimage/feature/censure.py +++ b/skimage/feature/censure.py @@ -111,7 +111,8 @@ def _suppress_lines(feature_mask, image, sigma, line_threshold): def keypoints_censure(image, min_scale=1, max_scale=7, mode='DoB', non_max_threshold=0.15, line_threshold=10): - """ + """**Experimental function**. + Extracts CenSurE keypoints along with the corresponding scale using either Difference of Boxes, Octagon or STAR bi-level filter. diff --git a/skimage/feature/tests/test_brief.py b/skimage/feature/tests/_test_brief.py similarity index 100% rename from skimage/feature/tests/test_brief.py rename to skimage/feature/tests/_test_brief.py diff --git a/skimage/feature/tests/test_censure.py b/skimage/feature/tests/_test_censure.py similarity index 100% rename from skimage/feature/tests/test_censure.py rename to skimage/feature/tests/_test_censure.py diff --git a/skimage/feature/tests/test_peak.py b/skimage/feature/tests/test_peak.py index 39dc8d8f..1a3e91f2 100644 --- a/skimage/feature/tests/test_peak.py +++ b/skimage/feature/tests/test_peak.py @@ -135,6 +135,137 @@ def test_ndarray_exclude_border(): assert (result == expected).all() +def test_empty(): + image = np.zeros((10, 20)) + labels = np.zeros((10, 20), int) + result = peak.peak_local_max(image, labels=labels, + footprint=np.ones((3, 3), bool), + min_distance=1, threshold_rel=0, + indices=False, exclude_border=False) + assert np.all(~ result) + + +def test_one_point(): + image = np.zeros((10, 20)) + labels = np.zeros((10, 20), int) + image[5, 5] = 1 + labels[5, 5] = 1 + result = peak.peak_local_max(image, labels=labels, + footprint=np.ones((3, 3), bool), + min_distance=1, threshold_rel=0, + indices=False, exclude_border=False) + assert np.all(result == (labels == 1)) + + +def test_adjacent_and_same(): + image = np.zeros((10, 20)) + labels = np.zeros((10, 20), int) + image[5, 5:6] = 1 + labels[5, 5:6] = 1 + result = peak.peak_local_max(image, labels=labels, + footprint=np.ones((3, 3), bool), + min_distance=1, threshold_rel=0, + indices=False, exclude_border=False) + assert np.all(result == (labels == 1)) + + +def test_adjacent_and_different(): + image = np.zeros((10, 20)) + labels = np.zeros((10, 20), int) + image[5, 5] = 1 + image[5, 6] = .5 + labels[5, 5:6] = 1 + expected = (image == 1) + result = peak.peak_local_max(image, labels=labels, + footprint=np.ones((3, 3), bool), + min_distance=1, threshold_rel=0, + indices=False, exclude_border=False) + assert np.all(result == expected) + result = peak.peak_local_max(image, labels=labels, + min_distance=1, threshold_rel=0, + indices=False, exclude_border=False) + assert np.all(result == expected) + + +def test_not_adjacent_and_different(): + image = np.zeros((10, 20)) + labels = np.zeros((10, 20), int) + image[5, 5] = 1 + image[5, 8] = .5 + labels[image > 0] = 1 + expected = (labels == 1) + result = peak.peak_local_max(image, labels=labels, + footprint=np.ones((3, 3), bool), + min_distance=1, threshold_rel=0, + indices=False, exclude_border=False) + assert np.all(result == expected) + + +def test_two_objects(): + image = np.zeros((10, 20)) + labels = np.zeros((10, 20), int) + image[5, 5] = 1 + image[5, 15] = .5 + labels[5, 5] = 1 + labels[5, 15] = 2 + expected = (labels > 0) + result = peak.peak_local_max(image, labels=labels, + footprint=np.ones((3, 3), bool), + min_distance=1, threshold_rel=0, + indices=False, exclude_border=False) + assert np.all(result == expected) + + +def test_adjacent_different_objects(): + image = np.zeros((10, 20)) + labels = np.zeros((10, 20), int) + image[5, 5] = 1 + image[5, 6] = .5 + labels[5, 5] = 1 + labels[5, 6] = 2 + expected = (labels > 0) + result = peak.peak_local_max(image, labels=labels, + footprint=np.ones((3, 3), bool), + min_distance=1, threshold_rel=0, + indices=False, exclude_border=False) + assert np.all(result == expected) + + +def test_four_quadrants(): + np.random.seed(21) + image = np.random.uniform(size=(40, 60)) + i, j = np.mgrid[0:40, 0:60] + labels = 1 + (i >= 20) + (j >= 30) * 2 + i, j = np.mgrid[-3:4, -3:4] + footprint = (i * i + j * j <= 9) + expected = np.zeros(image.shape, float) + for imin, imax in ((0, 20), (20, 40)): + for jmin, jmax in ((0, 30), (30, 60)): + expected[imin:imax, jmin:jmax] = scipy.ndimage.maximum_filter( + image[imin:imax, jmin:jmax], footprint=footprint) + expected = (expected == image) + result = peak.peak_local_max(image, labels=labels, footprint=footprint, + min_distance=1, threshold_rel=0, + indices=False, exclude_border=False) + assert np.all(result == expected) + + +def test_disk(): + '''regression test of img-1194, footprint = [1] + Test peak.peak_local_max when every point is a local maximum + ''' + np.random.seed(31) + image = np.random.uniform(size=(10, 20)) + footprint = np.array([[1]]) + result = peak.peak_local_max(image, labels=np.ones((10, 20)), + footprint=footprint, + min_distance=1, threshold_rel=0, + indices=False, exclude_border=False) + assert np.all(result) + result = peak.peak_local_max(image, footprint=footprint) + assert np.all(result) + + if __name__ == '__main__': from numpy import testing testing.run_module_suite() diff --git a/skimage/feature/util.py b/skimage/feature/util.py index eb3817e8..a5267d44 100644 --- a/skimage/feature/util.py +++ b/skimage/feature/util.py @@ -14,7 +14,9 @@ def _mask_border_keypoints(image, keypoints, dist): def pairwise_hamming_distance(array1, array2): - """Calculate hamming dissimilarity measure between two sets of + """**Experimental function**. + + Calculate hamming dissimilarity measure between two sets of vectors. Parameters diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index 67088b20..cef9b2e2 100644 --- a/skimage/filter/__init__.py +++ b/skimage/filter/__init__.py @@ -3,9 +3,9 @@ from .ctmf import median_filter from ._gaussian import gaussian_filter from ._canny import canny from .edges import (sobel, hsobel, vsobel, scharr, hscharr, vscharr, prewitt, - hprewitt, vprewitt, roberts , roberts_positive_diagonal, + hprewitt, vprewitt, roberts, roberts_positive_diagonal, roberts_negative_diagonal) -from ._denoise import denoise_tv_chambolle, tv_denoise +from ._denoise import denoise_tv_chambolle from ._denoise_cy import denoise_bilateral, denoise_tv_bregman from ._rank_order import rank_order from ._gabor import gabor_kernel, gabor_filter @@ -32,7 +32,6 @@ __all__ = ['inverse', 'roberts_positive_diagonal', 'roberts_negative_diagonal', 'denoise_tv_chambolle', - 'tv_denoise', 'denoise_bilateral', 'denoise_tv_bregman', 'rank_order', diff --git a/skimage/filter/_denoise.py b/skimage/filter/_denoise.py index ff89b78b..043399e0 100644 --- a/skimage/filter/_denoise.py +++ b/skimage/filter/_denoise.py @@ -1,6 +1,5 @@ import numpy as np from skimage import img_as_float -from skimage._shared.utils import deprecated def _denoise_tv_chambolle_3d(im, weight=100, eps=2.e-4, n_iter_max=200): @@ -250,14 +249,10 @@ def denoise_tv_chambolle(im, weight=50, eps=2.e-4, n_iter_max=200, out = np.zeros_like(im) for c in range(im.shape[2]): out[..., c] = _denoise_tv_chambolle_2d(im[..., c], weight, eps, - n_iter_max) + n_iter_max) else: out = _denoise_tv_chambolle_3d(im, weight, eps, n_iter_max) else: raise ValueError('only 2-d and 3-d images may be denoised with this ' 'function') return out - - -tv_denoise = deprecated('skimage.filter.denoise_tv_chambolle')\ - (denoise_tv_chambolle) diff --git a/skimage/filter/_rank_order.py b/skimage/filter/_rank_order.py index f878702f..cdd992ff 100644 --- a/skimage/filter/_rank_order.py +++ b/skimage/filter/_rank_order.py @@ -8,7 +8,7 @@ Copyright (c) 2009-2011 Broad Institute All rights reserved. Original author: Lee Kamentstky """ -import numpy +import numpy as np def rank_order(image): @@ -47,14 +47,14 @@ def rank_order(image): (array([0, 1, 2, 1], dtype=uint32), array([-1. , 2.5, 3.1])) """ flat_image = image.ravel() - sort_order = flat_image.argsort().astype(numpy.uint32) + sort_order = flat_image.argsort().astype(np.uint32) flat_image = flat_image[sort_order] - sort_rank = numpy.zeros_like(sort_order) + sort_rank = np.zeros_like(sort_order) is_different = flat_image[:-1] != flat_image[1:] - numpy.cumsum(is_different, out=sort_rank[1:]) - original_values = numpy.zeros((sort_rank[-1] + 1,), image.dtype) + np.cumsum(is_different, out=sort_rank[1:]) + original_values = np.zeros((sort_rank[-1] + 1,), image.dtype) original_values[0] = flat_image[0] original_values[1:] = flat_image[1:][is_different] - int_image = numpy.zeros_like(sort_order) + int_image = np.zeros_like(sort_order) int_image[sort_order] = sort_rank return (int_image.reshape(image.shape), original_values) diff --git a/skimage/filter/ctmf.py b/skimage/filter/ctmf.py index 267c5b6c..e39aa8bc 100644 --- a/skimage/filter/ctmf.py +++ b/skimage/filter/ctmf.py @@ -1,4 +1,4 @@ -'''ctmf.py - constant time per pixel median filtering with an octagonal shape +"""ctmf.py - constant time per pixel median filtering with an octagonal shape Reference: S. Perreault and P. Hebert, "Median Filtering in Constant Time", IEEE Transactions on Image Processing, September 2007. @@ -9,39 +9,46 @@ Copyright (c) 2003-2009 Massachusetts Institute of Technology Copyright (c) 2009-2011 Broad Institute All rights reserved. Original author: Lee Kamentsky -''' +""" +import warnings import numpy as np from . import _ctmf from ._rank_order import rank_order +from .._shared.utils import deprecated +@deprecated('filter.rank.median') def median_filter(image, radius=2, mask=None, percent=50): - '''Masked median filter with octagon shape. + """Masked median filter with octagon shape. Parameters ---------- - image : (M,N) ndarray, dtype uint8 + image : (M, N) ndarray Input image. - radius : {int, 2}, optional - The radius of a circle inscribed into the filtering - octagon. Must be at least 2. Default radius is 2. - mask : (M,N) ndarray, dtype uint8, optional - A value of 1 indicates a significant pixel, 0 - that a pixel is masked. By default, all pixels - are considered. - percent : {int, 50}, optional + radius : int + Radius (in pixels) of a circle inscribed into the filtering + octagon. Must be at least 2. Default radius is 2. + mask : (M, N) ndarray + Mask with 1's for significant pixels, 0's for masked pixels. + By default, all pixels are considered significant. + percent : int The unmasked pixels within the octagon are sorted, and the - value at the `percent`-th index chosen. For example, the - default value of 50 chooses the median pixel. + value at `percent` percent of the index range is chosen. + Default value of 50 gives the median pixel. Returns ------- - out : (M,N) ndarray, dtype uint8 - Filtered array. In areas where the median filter does - not overlap the mask, the filtered result is underfined, but + out : (M, N) ndarray + Filtered array. In areas where the median filter does + not overlap the mask, the filtered result is undefined, but in practice, it will be the lowest value in the valid area. + Notes + ----- + Because of the histogram implementation, the number of unique values + for the output is limited to 256. + Examples -------- >>> a = np.ones((5, 5)) @@ -49,53 +56,54 @@ def median_filter(image, radius=2, mask=None, percent=50): >>> b = median_filter(a) >>> b[2, 2] # the median filter is good at removing outliers 1.0 - ''' + """ if image.ndim != 2: - raise TypeError("The input 'image' must be a two dimensional array.") + raise TypeError("Input 'image' must be a two-dimensional array.") if radius < 2: - raise ValueError("The input 'radius' must be >= 2.") + raise ValueError("Input 'radius' must be >= 2.") if mask is None: mask = np.ones(image.shape, dtype=np.bool) mask = np.ascontiguousarray(mask, dtype=np.bool) if np.all(~ mask): + warnings.warn('Mask is all over image! Returning copy of input image.') return image.copy() - # - # Normalize the ranked image to 0-255 - # + if (not np.issubdtype(image.dtype, np.int) or np.min(image) < 0 or np.max(image) > 255): - ranked_image, translation = rank_order(image[mask]) - max_ranked_image = np.max(ranked_image) - if max_ranked_image == 0: - return image - if max_ranked_image > 255: - ranked_image = ranked_image * 255 // max_ranked_image + ranked_values, translation = rank_order(image[mask]) + max_ranked_values = np.max(ranked_values) + if max_ranked_values == 0: + warnings.warn('Particular case? Returning copy of input image.') + return image.copy() + if max_ranked_values > 255: + ranked_values = ranked_values * 255 // max_ranked_values was_ranked = True else: - ranked_image = image[mask] + ranked_values = image[mask] was_ranked = False - input = np.zeros(image.shape, np.uint8) - input[mask] = ranked_image + ranked_image = np.zeros(image.shape, np.uint8) + ranked_image[mask] = ranked_values mask.dtype = np.uint8 output = np.zeros(image.shape, np.uint8) - _ctmf.median_filter(input, mask, output, radius, percent) + _ctmf.median_filter(ranked_image, mask, output, radius, percent) if was_ranked: # # The translation gives the original value at each ranking. # We rescale the output to the original ranking and then # use the translation to look up the original value in the image. # - if max_ranked_image > 255: + if max_ranked_values > 255: result = translation[output.astype(np.uint32) * - max_ranked_image // 255] + max_ranked_values // 255] else: result = translation[output] else: result = output return result + diff --git a/skimage/filter/edges.py b/skimage/filter/edges.py index 7a70f00b..764c7d34 100644 --- a/skimage/filter/edges.py +++ b/skimage/filter/edges.py @@ -31,8 +31,8 @@ HPREWITT_WEIGHTS = np.array([[ 1, 1, 1], [-1,-1,-1]]) / 3.0 VPREWITT_WEIGHTS = HPREWITT_WEIGHTS.T -ROBERTS_PD_WEIGHTS = np.array([[ 1, 0], - [ 0, -1]], dtype=np.double) +ROBERTS_PD_WEIGHTS = np.array([[1, 0], + [0, -1]], dtype=np.double) ROBERTS_ND_WEIGHTS = np.array([[0, 1], [-1, 0]], dtype=np.double) @@ -346,7 +346,7 @@ def roberts(image, mask=None): """Find the edge magnitude using Roberts' cross operator. Parameters - ---------- + ---------- image : 2-D array Image to process. mask : 2-D array, optional diff --git a/skimage/io/_plugins/gtk_plugin.py b/skimage/io/_plugins/gtk_plugin.py index 961aaa0a..6130289c 100644 --- a/skimage/io/_plugins/gtk_plugin.py +++ b/skimage/io/_plugins/gtk_plugin.py @@ -5,7 +5,7 @@ try: # or else the gui import might trample another # gui's pyos_inputhook. window_manager.acquire('gtk') -except GuiLockError, gle: +except GuiLockError as gle: print(gle) else: try: diff --git a/skimage/io/_plugins/pil_plugin.py b/skimage/io/_plugins/pil_plugin.py index 3ce79c47..2ddbfe26 100644 --- a/skimage/io/_plugins/pil_plugin.py +++ b/skimage/io/_plugins/pil_plugin.py @@ -11,6 +11,8 @@ except ImportError: from skimage.util import img_as_ubyte +from skimage._shared import six + def imread(fname, dtype=None): """Load an image from file. @@ -104,7 +106,7 @@ def imsave(fname, arr, format_str=None): arr = arr.astype(np.uint8) # default to PNG if file-like object - if not isinstance(fname, basestring) and format_str is None: + if not isinstance(fname, six.string_types) and format_str is None: format_str = "PNG" img = Image.fromstring(mode, (arr.shape[1], arr.shape[0]), arr.tostring()) diff --git a/skimage/io/tests/test_pil.py b/skimage/io/tests/test_pil.py index 61ec5ce3..aa582ebc 100644 --- a/skimage/io/tests/test_pil.py +++ b/skimage/io/tests/test_pil.py @@ -8,7 +8,7 @@ from tempfile import NamedTemporaryFile from skimage import data_dir from skimage.io import (imread, imsave, use_plugin, reset_plugins, Image as ioImage) -from skimage._shared.six.moves import StringIO +from skimage._shared.six import BytesIO try: @@ -132,7 +132,7 @@ class TestSave: def test_imsave_filelike(): shape = (2, 2) image = np.zeros(shape) - s = StringIO() + s = BytesIO() # save to file-like object imsave(s, image) diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py index 772c52f1..c3ec9070 100644 --- a/skimage/measure/_marching_cubes.py +++ b/skimage/measure/_marching_cubes.py @@ -2,7 +2,7 @@ import numpy as np from . import _marching_cubes_cy -def marching_cubes(volume, level, sampling=(1., 1., 1.)): +def marching_cubes(volume, level, spacing=(1., 1., 1.)): """ Marching cubes algorithm to find iso-valued surfaces in 3d volumetric data @@ -12,7 +12,7 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)): Input data volume to find isosurfaces. Will be cast to `np.float64`. level : float Contour value to search for isosurfaces in `volume`. - sampling : length-3 tuple of floats + spacing : length-3 tuple of floats Voxel spacing in spatial dimensions corresponding to numpy array indexing dimensions (M, N, P) as in `volume`. @@ -34,7 +34,7 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)): http://www.essi.fr/~lingrand/MarchingCubes/algo.html There are several known ambiguous cases in the marching cubes algorithm. - Using point labeling as in [1]_, Figure 4, as shown: + Using point labeling as in [1]_, Figure 4, as shown:: v8 ------ v7 / | / | y @@ -72,15 +72,15 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)): the 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: + is recommended. To contour a volume named `myvolume` about the level 0.0:: - >>> from mayavi import mlab - >>> verts, tris = marching_cubes(myvolume, 0.0, (1., 1., 2.)) - >>> mlab.triangular_mesh([vert[0] for vert in verts], - [vert[1] for vert in verts], - [vert[2] for vert in verts], - tris) - >>> mlab.show() + >>> from mayavi import mlab + >>> verts, tris = marching_cubes(myvolume, 0.0, (1., 1., 2.)) + >>> mlab.triangular_mesh([vert[0] for vert in verts], + ... [vert[1] for vert in verts], + ... [vert[2] for vert in verts], + ... tris) + >>> mlab.show() References ---------- @@ -107,7 +107,7 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)): # have repeated vertices - and equivalent vertices are redundantly # placed in every triangle they connect with. raw_tris = _marching_cubes_cy.iterate_and_store_3d(volume, float(level), - sampling) + spacing) # Find and collect unique vertices, storing triangle verts as indices. # Returns a true mesh with no degenerate faces. diff --git a/skimage/measure/_marching_cubes_cy.pyx b/skimage/measure/_marching_cubes_cy.pyx index a0f63d1c..085108ab 100644 --- a/skimage/measure/_marching_cubes_cy.pyx +++ b/skimage/measure/_marching_cubes_cy.pyx @@ -56,32 +56,32 @@ def unpack_unique_verts(list trilist): def iterate_and_store_3d(double[:, :, ::1] arr, double level, - tuple sampling=(1., 1., 1.)): + tuple spacing=(1., 1., 1.)): """Iterate across the given array in a marching-cubes fashion, looking for volumes with edges that cross 'level'. If such a volume is found, appropriate triangulations are added to a growing list of faces to be returned by this function. - If `sampling` is not provided, vertices are returned in the indexing + If `spacing` is not provided, vertices are returned in the indexing coordinate system (assuming all 3 spatial dimensions sampled equally). - If `sampling` is provided, vertices will be returned in volume coordinates + If `spacing` is provided, vertices will be returned in volume coordinates relative to the origin, regularly spaced as specified in each dimension. """ if arr.shape[0] < 2 or arr.shape[1] < 2 or arr.shape[2] < 2: raise ValueError("Input array must be at least 2x2x2.") - if len(sampling) != 3: - raise ValueError("`sampling` must be (double, double, double)") + if len(spacing) != 3: + raise ValueError("`spacing` must be (double, double, double)") cdef list face_list = [] cdef list norm_list = [] cdef Py_ssize_t n - cdef bint odd_sampling, plus_z + cdef bint odd_spacing, plus_z plus_z = False - if [float(i) for i in sampling] == [1.0, 1.0, 1.0]: - odd_sampling = False + if [float(i) for i in spacing] == [1.0, 1.0, 1.0]: + odd_spacing = False else: - odd_sampling = True + odd_spacing = True # The plan is to iterate a 2x2x2 cube across the input array. This means # the upper-left corner of the cube needs to iterate across a sub-array @@ -107,11 +107,11 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level, coords[1] = 0 coords[2] = 0 - # Extract doubles from `sampling` for speed - cdef double[3] sampling2 - sampling2[0] = sampling[0] - sampling2[1] = sampling[1] - sampling2[2] = sampling[2] + # Extract doubles from `spacing` for speed + cdef double[3] spacing2 + spacing2[0] = spacing[0] + spacing2[1] = spacing[1] + spacing2[2] = spacing[2] # Calculate the number of iterations we'll need cdef Py_ssize_t num_cube_steps = ((arr.shape[0] - 1) * @@ -138,15 +138,15 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level, x0, y0, z0 = coords[0], coords[1], coords[2] x1, y1, z1 = x0 + 1, y0 + 1, z0 + 1 - if odd_sampling: + if odd_spacing: # These doubles are the modified world coordinates; they are only - # calculated if non-default `sampling` provided. - r0 = coords[0] * sampling2[0] - c0 = coords[1] * sampling2[1] - d0 = coords[2] * sampling2[2] - r1 = r0 + sampling2[0] - c1 = c0 + sampling2[1] - d1 = d0 + sampling2[2] + # calculated if non-default `spacing` provided. + r0 = coords[0] * spacing2[0] + c0 = coords[1] * spacing2[1] + d0 = coords[2] * spacing2[2] + r1 = r0 + spacing2[0] + c1 = c0 + spacing2[1] + d1 = d0 + spacing2[2] else: r0, c0, d0, r1, c1, d1 = x0, y0, z0, x1, y1, z1 @@ -193,11 +193,11 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level, e4 = e8 else: # Calculate edges normally - if odd_sampling: - e1 = r0 + _get_fraction(v1, v2, level) * sampling2[0], c0, d0 - e2 = r1, c0 + _get_fraction(v2, v3, level) * sampling2[1], d0 - e3 = r0 + _get_fraction(v4, v3, level) * sampling2[0], c1, d0 - e4 = r0, c0 + _get_fraction(v1, v4, level) * sampling2[1], d0 + if odd_spacing: + e1 = r0 + _get_fraction(v1, v2, level) * spacing2[0], c0, d0 + e2 = r1, c0 + _get_fraction(v2, v3, level) * spacing2[1], d0 + e3 = r0 + _get_fraction(v4, v3, level) * spacing2[0], c1, d0 + e4 = r0, c0 + _get_fraction(v1, v4, level) * spacing2[1], d0 else: e1 = r0 + _get_fraction(v1, v2, level), c0, d0 e2 = r1, c0 + _get_fraction(v2, v3, level), d0 @@ -208,15 +208,15 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level, # large, growing lookup table for all adjacent values; could save # ~30% in terms of runtime at the expense of memory usage and # much greater complexity. - if odd_sampling: - e5 = r0 + _get_fraction(v5, v6, level) * sampling2[0], c0, d1 - e6 = r1, c0 + _get_fraction(v6, v7, level) * sampling2[1], d1 - e7 = r0 + _get_fraction(v8, v7, level) * sampling2[0], c1, d1 - e8 = r0, c0 + _get_fraction(v5, v8, level) * sampling2[1], d1 - e9 = r0, c0, d0 + _get_fraction(v1, v5, level) * sampling2[2] - e10 = r1, c0, d0 + _get_fraction(v2, v6, level) * sampling2[2] - e11 = r0, c1, d0 + _get_fraction(v4, v8, level) * sampling2[2] - e12 = r1, c1, d0 + _get_fraction(v3, v7, level) * sampling2[2] + if odd_spacing: + e5 = r0 + _get_fraction(v5, v6, level) * spacing2[0], c0, d1 + e6 = r1, c0 + _get_fraction(v6, v7, level) * spacing2[1], d1 + e7 = r0 + _get_fraction(v8, v7, level) * spacing2[0], c1, d1 + e8 = r0, c0 + _get_fraction(v5, v8, level) * spacing2[1], d1 + e9 = r0, c0, d0 + _get_fraction(v1, v5, level) * spacing2[2] + e10 = r1, c0, d0 + _get_fraction(v2, v6, level) * spacing2[2] + e11 = r0, c1, d0 + _get_fraction(v4, v8, level) * spacing2[2] + e12 = r1, c1, d0 + _get_fraction(v3, v7, level) * spacing2[2] else: e5 = r0 + _get_fraction(v5, v6, level), c0, d1 e6 = r1, c0 + _get_fraction(v6, v7, level), d1 diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index 9078df6f..798e04cc 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -4,11 +4,13 @@ from math import sqrt, atan2, pi as PI import numpy as np from scipy import ndimage -from skimage.morphology import convex_hull_image +from collections import MutableMapping + +from skimage.morphology import convex_hull_image, label from skimage.measure import _moments -__all__ = ['regionprops'] +__all__ = ['regionprops', 'perimeter'] STREL_4 = np.array([[0, 1, 0], @@ -47,7 +49,7 @@ PROPS = { # 'PixelList', 'Solidity': 'solidity', # 'SubarrayIdx' - 'WeightedCentralMoments': 'weighted_central_moments', + 'WeightedCentralMoments': 'weighted_moments_central', 'WeightedCentroid': 'weighted_centroid', 'WeightedHuMoments': 'weighted_moments_hu', 'WeightedMoments': 'weighted_moments', @@ -103,15 +105,16 @@ class _cached_property(object): return value -class _RegionProperties(object): +class _RegionProperties(MutableMapping): def __init__(self, slice, label, label_image, intensity_image, - cache_active): + cache_active, properties=None): self.label = label self._slice = slice self._label_image = label_image self._intensity_image = intensity_image self._cache_active = cache_active + self._properties = properties @_cached_property def area(self): @@ -155,8 +158,8 @@ class _RegionProperties(object): @_cached_property def euler_number(self): euler_array = self.filled_image != self.image - _, num = ndimage.label(euler_array, STREL_8) - return -num + _, num = label(euler_array, neighbors=8, return_num=True) + return -num + 1 @_cached_property def extent(self): @@ -288,7 +291,7 @@ class _RegionProperties(object): return _moments.moments_central(self._intensity_image_double, 0, 0, 3) @_cached_property - def weighted_central_moments(self): + def weighted_moments_central(self): row, col = self.weighted_local_centroid return _moments.moments_central(self._intensity_image_double, row, col, 3) @@ -299,7 +302,21 @@ class _RegionProperties(object): @_cached_property def weighted_moments_normalized(self): - return _moments.moments_normalized(self.weighted_central_moments, 3) + return _moments.moments_normalized(self.weighted_moments_central, 3) + + + # Preserve dictionary interface + def __delitem__(self, key): + pass + + def __len__(self): + return len(self._properties or PROPS.values()) + + def __setitem__(self, key, value): + raise RuntimeError("Cannot assign region properties.") + + def __iter__(self): + return iter(self._properties or PROPS.values()) def __getitem__(self, key): value = getattr(self, key, None) @@ -430,7 +447,7 @@ def regionprops(label_image, properties=None, wm_ji = sum{ array(x, y) * x^j * y^i } where the sum is over the `x`, `y` coordinates of the region. - **weighted_central_moments** : (3, 3) ndarray + **weighted_moments_central** : (3, 3) ndarray Central moments (translation invariant) of intensity image up to 3rd order:: @@ -489,7 +506,7 @@ def regionprops(label_image, properties=None, label = i + 1 props = _RegionProperties(sl, label, label_image, - intensity_image, cache) + intensity_image, cache, properties=properties) regions.append(props) return regions diff --git a/skimage/measure/tests/test_marching_cubes.py b/skimage/measure/tests/test_marching_cubes.py index 7a1fd40a..b3c2ddc1 100644 --- a/skimage/measure/tests/test_marching_cubes.py +++ b/skimage/measure/tests/test_marching_cubes.py @@ -16,12 +16,12 @@ def test_marching_cubes_isotropic(): def test_marching_cubes_anisotropic(): - sampling = (1., 10 / 6., 16 / 6.) - ellipsoid_anisotropic = ellipsoid(6, 10, 16, sampling=sampling, + spacing = (1., 10 / 6., 16 / 6.) + ellipsoid_anisotropic = ellipsoid(6, 10, 16, spacing=spacing, levelset=True) - _, surf = ellipsoid_stats(6, 10, 16, sampling=sampling) + _, surf = ellipsoid_stats(6, 10, 16) verts, faces = marching_cubes(ellipsoid_anisotropic, 0., - sampling=sampling) + spacing=spacing) surf_calc = mesh_surface_area(verts, faces) # Test within 1.5% tolerance for anisotropic. Will always underestimate. @@ -32,7 +32,7 @@ def test_invalid_input(): assert_raises(ValueError, marching_cubes, np.zeros((2, 2, 1)), 0) assert_raises(ValueError, marching_cubes, np.zeros((2, 2, 1)), 1) assert_raises(ValueError, marching_cubes, np.ones((3, 3, 3)), 1, - sampling=(1, 2)) + spacing=(1, 2)) assert_raises(ValueError, marching_cubes, np.zeros((20, 20)), 0) diff --git a/skimage/measure/tests/test_regionprops.py b/skimage/measure/tests/test_regionprops.py index 647c09a7..c0c6443d 100644 --- a/skimage/measure/tests/test_regionprops.py +++ b/skimage/measure/tests/test_regionprops.py @@ -1,5 +1,5 @@ from numpy.testing import assert_array_equal, assert_almost_equal, \ - assert_array_almost_equal, assert_raises + assert_array_almost_equal, assert_raises, assert_equal import numpy as np import math @@ -31,7 +31,8 @@ def test_all_props(): def test_dtype(): regionprops(np.zeros((10, 10), dtype=np.int)) regionprops(np.zeros((10, 10), dtype=np.uint)) - assert_raises(TypeError, regionprops, np.zeros((10, 10), dtype=np.double)) + assert_raises((TypeError, RuntimeError), regionprops, + np.zeros((10, 10), dtype=np.double)) def test_ndim(): @@ -269,7 +270,7 @@ def test_solidity(): assert_almost_equal(solidity, 0.580645161290323) -def test_weighted_moments(): +def test_weighted_moments_central(): wmu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE )[0].weighted_moments_central ref = np.array( @@ -335,6 +336,17 @@ def test_weighted_moments_normalized(): assert_array_almost_equal(wnu, ref) +def test_old_dict_interface(): + feats = regionprops(SAMPLE, + ['Area', 'Eccentricity', 'EulerNumber', + 'Extent', 'MinIntensity', 'MeanIntensity', + 'MaxIntensity', 'Solidity'], + intensity_image=INTENSITY_SAMPLE) + + np.array([list(props.values()) for props in feats], np.float) + assert_equal(len(feats[0]), 8) + + 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 f01fc3ed..4788c308 100644 --- a/skimage/morphology/__init__.py +++ b/skimage/morphology/__init__.py @@ -7,7 +7,7 @@ from .grey import (erosion, dilation, opening, closing, white_tophat, from .selem import (square, rectangle, diamond, disk, cube, octahedron, ball, octagon, star) from .ccomp import label -from .watershed import watershed, is_local_maximum +from .watershed import watershed from ._skeletonize import skeletonize, medial_axis from .convex_hull import convex_hull_image, convex_hull_object from .greyreconstruct import reconstruction @@ -40,7 +40,6 @@ __all__ = ['binary_erosion', 'octagon', 'label', 'watershed', - 'is_local_maximum', 'skeletonize', 'medial_axis', 'convex_hull_image', diff --git a/skimage/morphology/_skeletonize.py b/skimage/morphology/_skeletonize.py index 8872aada..11181a7b 100644 --- a/skimage/morphology/_skeletonize.py +++ b/skimage/morphology/_skeletonize.py @@ -282,8 +282,8 @@ def medial_axis(image, mask=None, return_distance=False): i, j = np.mgrid[0:image.shape[0], 0:image.shape[1]] result = masked_image.copy() distance = distance[result] - i = np.ascontiguousarray(i[result], np.intp) - j = np.ascontiguousarray(j[result], np.intp) + i = np.ascontiguousarray(i[result], dtype=np.intp) + j = np.ascontiguousarray(j[result], dtype=np.intp) result = np.ascontiguousarray(result, np.uint8) # Determine the order in which pixels are processed. @@ -296,9 +296,9 @@ def medial_axis(image, mask=None, return_distance=False): order = np.lexsort((tiebreaker, corner_score[masked_image], distance)) - order = np.ascontiguousarray(order, np.int32) + order = np.ascontiguousarray(order, dtype=np.int32) - table = np.ascontiguousarray(table, np.uint8) + table = np.ascontiguousarray(table, dtype=np.uint8) # Remove pixels not belonging to the medial axis _skeletonize_loop(result, i, j, order, table) diff --git a/skimage/morphology/_skeletonize_cy.pyx b/skimage/morphology/_skeletonize_cy.pyx index 10e70d98..c4471a11 100644 --- a/skimage/morphology/_skeletonize_cy.pyx +++ b/skimage/morphology/_skeletonize_cy.pyx @@ -20,8 +20,8 @@ cimport numpy as cnp def _skeletonize_loop(cnp.uint8_t[:, ::1] result, - Py_ssize_t[:] i, Py_ssize_t[:] j, - cnp.int32_t[:] order, cnp.uint8_t[:] table): + Py_ssize_t[::1] i, Py_ssize_t[::1] j, + cnp.int32_t[::1] order, cnp.uint8_t[::1] table): """ Inner loop of skeletonize function diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx index beaf874b..a0857707 100644 --- a/skimage/morphology/_watershed.pyx +++ b/skimage/morphology/_watershed.pyx @@ -23,12 +23,12 @@ include "heap_watershed.pxi" @cython.boundscheck(False) -def watershed(DTYPE_INT32_t[:] image, +def watershed(DTYPE_INT32_t[::1] image, DTYPE_INT32_t[:, ::1] pq, Py_ssize_t age, DTYPE_INT32_t[:, ::1] structure, - DTYPE_BOOL_t[:] mask, - DTYPE_INT32_t[:] output): + DTYPE_BOOL_t[::1] mask, + DTYPE_INT32_t[::1] output): """Do heavy lifting of watershed algorithm Parameters diff --git a/skimage/morphology/binary.py b/skimage/morphology/binary.py index e2e0f20b..24c5c0ea 100644 --- a/skimage/morphology/binary.py +++ b/skimage/morphology/binary.py @@ -1,3 +1,4 @@ +import warnings import numpy as np from scipy import ndimage @@ -8,32 +9,40 @@ def binary_erosion(image, selem, out=None): This function returns the same result as greyscale erosion but performs faster for binary images. - Morphological erosion sets a pixel at (i,j) to the minimum over all pixels - in the neighborhood centered at (i,j). Erosion shrinks bright regions and - enlarges dark regions. + Morphological erosion sets a pixel at ``(i,j)`` to the minimum over all + pixels in the neighborhood centered at ``(i,j)``. Erosion shrinks bright + regions and enlarges dark regions. Parameters ---------- image : ndarray - Image array. + Binary input image. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - out : ndarray + out : ndarray of bool The array to store the result of the morphology. If None is passed, a new array will be allocated. Returns ------- - eroded : bool array - The result of the morphological erosion. + eroded : ndarray of bool or uint + The result of the morphological erosion with values in ``[0, 1]``. """ + selem = (selem != 0) + selem_sum = np.sum(selem) - conv = ndimage.convolve(image > 0, selem, output=out, - mode='constant', cval=1) - if conv is not None: + if selem_sum <= 255: + conv = np.empty_like(image, dtype=np.uint8) + else: + conv = np.empty_like(image, dtype=np.uint) + + binary = (image > 0).view(np.uint8) + ndimage.convolve(binary, selem, mode='constant', cval=1, output=conv) + + if out is None: out = conv - return np.equal(out, np.sum(selem), out=out) + return np.equal(conv, selem_sum, out=out) def binary_dilation(image, selem, out=None): @@ -42,33 +51,40 @@ def binary_dilation(image, selem, out=None): This function returns the same result as greyscale dilation but performs faster for binary images. - Morphological dilation sets a pixel at (i,j) to the maximum over all pixels - in the neighborhood centered at (i,j). Dilation enlarges bright regions - and shrinks dark regions. + Morphological dilation sets a pixel at ``(i,j)`` to the maximum over all + pixels in the neighborhood centered at ``(i,j)``. Dilation enlarges bright + regions and shrinks dark regions. Parameters ---------- image : ndarray - Image array. + Binary input image. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - out : ndarray + out : ndarray of bool The array to store the result of the morphology. If None, is passed, a new array will be allocated. Returns ------- - dilated : bool array - The result of the morphological dilation. + dilated : ndarray of bool or uint + The result of the morphological dilation with values in ``[0, 1]``. """ + selem = (selem != 0) - conv = ndimage.convolve(image > 0, selem, output=out, - mode='constant', cval=0) - if conv is not None: + if np.sum(selem) <= 255: + conv = np.empty_like(image, dtype=np.uint8) + else: + conv = np.empty_like(image, dtype=np.uint) + + binary = (image > 0).view(np.uint8) + ndimage.convolve(binary, selem, mode='constant', cval=0, output=conv) + + if out is None: out = conv - return np.not_equal(out, 0, out=out) + return np.not_equal(conv, 0, out=out) def binary_opening(image, selem, out=None): @@ -85,20 +101,19 @@ def binary_opening(image, selem, out=None): Parameters ---------- image : ndarray - Image array. + Binary input image. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - out : ndarray + out : ndarray of bool The array to store the result of the morphology. If None is passed, a new array will be allocated. Returns ------- - opening : bool array + opening : ndarray of bool The result of the morphological opening. """ - eroded = binary_erosion(image, selem) out = binary_dilation(eroded, selem, out=out) return out @@ -118,16 +133,16 @@ def binary_closing(image, selem, out=None): Parameters ---------- image : ndarray - Image array. + Binary input image. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - out : ndarray + out : ndarray of bool The array to store the result of the morphology. If None, is passed, a new array will be allocated. Returns ------- - closing : bool array + closing : ndarray of bool The result of the morphological closing. """ diff --git a/skimage/morphology/convex_hull.py b/skimage/morphology/convex_hull.py index 3d592dca..adc63607 100644 --- a/skimage/morphology/convex_hull.py +++ b/skimage/morphology/convex_hull.py @@ -44,7 +44,7 @@ def convex_hull_image(image): (-0.5, 0.5, 0, 0))): coords_corners[i * N:(i + 1) * N] = coords + [x_offset, y_offset] - # repeated coordinates can *sometimes* cause problems in + # repeated coordinates can *sometimes* cause problems in # scipy.spatial.Delaunay, so we remove them. coords = unique_rows(coords_corners) diff --git a/skimage/morphology/tests/test_binary.py b/skimage/morphology/tests/test_binary.py new file mode 100644 index 00000000..deab3d82 --- /dev/null +++ b/skimage/morphology/tests/test_binary.py @@ -0,0 +1,68 @@ +import numpy as np +from numpy import testing + +from skimage import data, color +from skimage.util import img_as_bool +from skimage.morphology import binary, grey, selem + + +lena = color.rgb2gray(data.lena()) +bw_lena = lena > 100 + + +def test_non_square_image(): + strel = selem.square(3) + binary_res = binary.binary_erosion(bw_lena[:100, :200], strel) + grey_res = img_as_bool(grey.erosion(bw_lena[:100, :200], strel)) + testing.assert_array_equal(binary_res, grey_res) + + +def test_binary_erosion(): + strel = selem.square(3) + binary_res = binary.binary_erosion(bw_lena, strel) + grey_res = img_as_bool(grey.erosion(bw_lena, strel)) + testing.assert_array_equal(binary_res, grey_res) + + +def test_binary_dilation(): + strel = selem.square(3) + binary_res = binary.binary_dilation(bw_lena, strel) + grey_res = img_as_bool(grey.dilation(bw_lena, strel)) + testing.assert_array_equal(binary_res, grey_res) + + +def test_binary_closing(): + strel = selem.square(3) + binary_res = binary.binary_closing(bw_lena, strel) + grey_res = img_as_bool(grey.closing(bw_lena, strel)) + testing.assert_array_equal(binary_res, grey_res) + + +def test_binary_opening(): + strel = selem.square(3) + binary_res = binary.binary_opening(bw_lena, strel) + grey_res = img_as_bool(grey.opening(bw_lena, strel)) + testing.assert_array_equal(binary_res, grey_res) + + +def test_selem_overflow(): + strel = np.ones((17, 17), dtype=np.uint8) + img = np.zeros((20, 20)) + img[2:19, 2:19] = 1 + binary_res = binary.binary_erosion(img, strel) + grey_res = img_as_bool(grey.erosion(img, strel)) + testing.assert_array_equal(binary_res, grey_res) + + +def test_out_argument(): + for func in (binary.binary_erosion, binary.binary_dilation): + strel = np.ones((3, 3), dtype=np.uint8) + img = np.ones((10, 10)) + out = np.zeros_like(img) + out_saved = out.copy() + func(img, strel, out=out) + testing.assert_(np.any(out != out_saved)) + testing.assert_array_equal(out, func(img, strel)) + +if __name__ == '__main__': + testing.run_module_suite() diff --git a/skimage/morphology/tests/test_grey.py b/skimage/morphology/tests/test_grey.py index 244ec566..e2a3928d 100644 --- a/skimage/morphology/tests/test_grey.py +++ b/skimage/morphology/tests/test_grey.py @@ -6,7 +6,7 @@ from numpy import testing import skimage from skimage import data_dir from skimage.util import img_as_bool -from skimage.morphology import binary, grey, selem +from skimage.morphology import grey, selem lena = np.load(os.path.join(data_dir, 'lena_GRAY_U8.npy')) @@ -155,40 +155,5 @@ class TestDTypes(): self._test_image(image) -def test_non_square_image(): - strel = selem.square(3) - binary_res = binary.binary_erosion(bw_lena[:100, :200], strel) - grey_res = img_as_bool(grey.erosion(bw_lena[:100, :200], strel)) - testing.assert_array_equal(binary_res, grey_res) - - -def test_binary_erosion(): - strel = selem.square(3) - binary_res = binary.binary_erosion(bw_lena, strel) - grey_res = img_as_bool(grey.erosion(bw_lena, strel)) - testing.assert_array_equal(binary_res, grey_res) - - -def test_binary_dilation(): - strel = selem.square(3) - binary_res = binary.binary_dilation(bw_lena, strel) - grey_res = img_as_bool(grey.dilation(bw_lena, strel)) - testing.assert_array_equal(binary_res, grey_res) - - -def test_binary_closing(): - strel = selem.square(3) - binary_res = binary.binary_closing(bw_lena, strel) - grey_res = img_as_bool(grey.closing(bw_lena, strel)) - testing.assert_array_equal(binary_res, grey_res) - - -def test_binary_opening(): - strel = selem.square(3) - binary_res = binary.binary_opening(bw_lena, strel) - grey_res = img_as_bool(grey.opening(bw_lena, strel)) - testing.assert_array_equal(binary_res, grey_res) - - if __name__ == '__main__': testing.run_module_suite() diff --git a/skimage/morphology/tests/test_watershed.py b/skimage/morphology/tests/test_watershed.py index 82531713..5dc0f07c 100644 --- a/skimage/morphology/tests/test_watershed.py +++ b/skimage/morphology/tests/test_watershed.py @@ -48,8 +48,7 @@ import unittest import numpy as np import scipy.ndimage -from skimage.morphology.watershed import watershed, \ - _slow_watershed, is_local_maximum +from skimage.morphology.watershed import watershed, _slow_watershed eps = 1e-12 @@ -387,101 +386,5 @@ class TestWatershed(unittest.TestCase): self.eight) -class TestIsLocalMaximum(unittest.TestCase): - def test_00_00_empty(self): - image = np.zeros((10, 20)) - labels = np.zeros((10, 20), int) - result = is_local_maximum(image, labels, np.ones((3, 3), bool)) - self.assertTrue(np.all(~ result)) - - def test_01_01_one_point(self): - image = np.zeros((10, 20)) - labels = np.zeros((10, 20), int) - image[5, 5] = 1 - labels[5, 5] = 1 - result = is_local_maximum(image, labels, np.ones((3, 3), bool)) - self.assertTrue(np.all(result == (labels == 1))) - - def test_01_02_adjacent_and_same(self): - image = np.zeros((10, 20)) - labels = np.zeros((10, 20), int) - image[5, 5:6] = 1 - labels[5, 5:6] = 1 - result = is_local_maximum(image, labels, np.ones((3, 3), bool)) - self.assertTrue(np.all(result == (labels == 1))) - - def test_01_03_adjacent_and_different(self): - image = np.zeros((10, 20)) - labels = np.zeros((10, 20), int) - image[5, 5] = 1 - image[5, 6] = .5 - labels[5, 5:6] = 1 - expected = (image == 1) - result = is_local_maximum(image, labels, np.ones((3, 3), bool)) - self.assertTrue(np.all(result == expected)) - result = is_local_maximum(image, labels) - self.assertTrue(np.all(result == expected)) - - def test_01_04_not_adjacent_and_different(self): - image = np.zeros((10, 20)) - labels = np.zeros((10, 20), int) - image[5, 5] = 1 - image[5, 8] = .5 - labels[image > 0] = 1 - expected = (labels == 1) - result = is_local_maximum(image, labels, np.ones((3, 3), bool)) - self.assertTrue(np.all(result == expected)) - - def test_01_05_two_objects(self): - image = np.zeros((10, 20)) - labels = np.zeros((10, 20), int) - image[5, 5] = 1 - image[5, 15] = .5 - labels[5, 5] = 1 - labels[5, 15] = 2 - expected = (labels > 0) - result = is_local_maximum(image, labels, np.ones((3, 3), bool)) - self.assertTrue(np.all(result == expected)) - - def test_01_06_adjacent_different_objects(self): - image = np.zeros((10, 20)) - labels = np.zeros((10, 20), int) - image[5, 5] = 1 - image[5, 6] = .5 - labels[5, 5] = 1 - labels[5, 6] = 2 - expected = (labels > 0) - result = is_local_maximum(image, labels, np.ones((3, 3), bool)) - self.assertTrue(np.all(result == expected)) - - def test_02_01_four_quadrants(self): - np.random.seed(21) - image = np.random.uniform(size=(40, 60)) - i, j = np.mgrid[0:40, 0:60] - labels = 1 + (i >= 20) + (j >= 30) * 2 - i, j = np.mgrid[-3:4, -3:4] - footprint = (i * i + j * j <= 9) - expected = np.zeros(image.shape, float) - for imin, imax in ((0, 20), (20, 40)): - for jmin, jmax in ((0, 30), (30, 60)): - expected[imin:imax, jmin:jmax] = scipy.ndimage.maximum_filter( - image[imin:imax, jmin:jmax], footprint=footprint) - expected = (expected == image) - result = is_local_maximum(image, labels, footprint) - self.assertTrue(np.all(result == expected)) - - def test_03_01_disk_1(self): - '''regression test of img-1194, footprint = [1] - - Test is_local_maximum when every point is a local maximum - ''' - np.random.seed(31) - image = np.random.uniform(size=(10, 20)) - footprint = np.array([[1]]) - result = is_local_maximum(image, np.ones((10, 20)), footprint) - self.assertTrue(np.all(result)) - result = is_local_maximum(image, footprint=footprint) - self.assertTrue(np.all(result)) - if __name__ == "__main__": np.testing.run_module_suite() diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py index 097e3129..2bdb57b2 100644 --- a/skimage/morphology/watershed.py +++ b/skimage/morphology/watershed.py @@ -116,7 +116,9 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None): >>> # to the background >>> from scipy import ndimage >>> distance = ndimage.distance_transform_edt(image) - >>> local_maxi = is_local_maximum(distance, image, np.ones((3, 3))) + >>> from skimage.feature import peak_local_max + >>> local_maxi = peak_local_max(distance, labels=image, + ... footprint=np.ones((3, 3))) >>> markers = ndimage.label(local_maxi)[0] >>> labels = watershed(-distance, markers, mask=image) @@ -200,7 +202,7 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None): stride = np.dot(image_stride, np.array(offs)) offs.insert(0, stride) c.append(offs) - c = np.array(c, np.int32) + c = np.array(c, dtype=np.int32) pq, age = __heapify_markers(c_markers, c_image) pq = np.ascontiguousarray(pq, dtype=np.int32) @@ -224,79 +226,6 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None): return c_output -@deprecated('feature.peak_local_max') -def is_local_maximum(image, labels=None, footprint=None): - """ - Return a boolean array of points that are local maxima - - Parameters - ---------- - image: ndarray (2-D, 3-D, ...) - intensity image - labels: ndarray, optional - find maxima only within labels. Zero is reserved for background. - footprint: ndarray of bools, optional - binary mask indicating the neighborhood to be examined - `footprint` must be a matrix with odd dimensions, the center is taken - to be the point in question. - - Returns - ------- - result: ndarray of bools - mask that is True for pixels that are local maxima of `image` - - See also - -------- - skimage.feature.peak_local_max: Unified peak finding backend. - The more capable backend for finding local maxima. - - Notes - ----- - This function is now a wrapper for skimage.feature.peak_local_max() and is - retained only for convenience and backward compatibility. - - Examples - -------- - >>> image = np.zeros((4, 4)) - >>> image[1, 2] = 2 - >>> image[3, 3] = 1 - >>> image - array([[ 0., 0., 0., 0.], - [ 0., 0., 2., 0.], - [ 0., 0., 0., 0.], - [ 0., 0., 0., 1.]]) - >>> is_local_maximum(image) - array([[ True, False, False, False], - [ True, False, True, False], - [ True, False, False, False], - [ True, True, False, True]], dtype=bool) - >>> image = np.arange(16).reshape((4, 4)) - >>> labels = np.array([[1, 2], [3, 4]]) - >>> labels = np.repeat(np.repeat(labels, 2, axis=0), 2, axis=1) - >>> labels - array([[1, 1, 2, 2], - [1, 1, 2, 2], - [3, 3, 4, 4], - [3, 3, 4, 4]]) - >>> image - array([[ 0, 1, 2, 3], - [ 4, 5, 6, 7], - [ 8, 9, 10, 11], - [12, 13, 14, 15]]) - >>> is_local_maximum(image, labels=labels) - array([[False, False, False, False], - [False, True, False, True], - [False, False, False, False], - [False, True, False, True]], dtype=bool) - - """ - # call import here to prevent circular imports - from ..feature import peak_local_max - return peak_local_max(image, labels=labels, min_distance=1, - threshold_rel=0, footprint=footprint, - indices=False, exclude_border=False) - - # ---------------------- deprecated ------------------------------ # Deprecate slower pure-Python code, that we keep only for # pedagogical purposes diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py index 107111bd..aea6c70f 100644 --- a/skimage/segmentation/__init__.py +++ b/skimage/segmentation/__init__.py @@ -4,7 +4,7 @@ from .slic_superpixels import slic from ._quickshift import quickshift from .boundaries import find_boundaries, visualize_boundaries, mark_boundaries from ._clear_border import clear_border -from ._join import join_segmentations, relabel_from_one +from ._join import join_segmentations, relabel_from_one, relabel_sequential __all__ = ['random_walker', @@ -16,4 +16,5 @@ __all__ = ['random_walker', 'mark_boundaries', 'clear_border', 'join_segmentations', - 'relabel_from_one'] + 'relabel_from_one', + 'relabel_sequential'] diff --git a/skimage/segmentation/_join.py b/skimage/segmentation/_join.py index 1eda9176..6095382d 100644 --- a/skimage/segmentation/_join.py +++ b/skimage/segmentation/_join.py @@ -1,12 +1,13 @@ import numpy as np +from skimage._shared.utils import deprecated def join_segmentations(s1, s2): """Return the join of the two input segmentations. - The join J of S1 and S2 is defined as the segmentation in which two voxels - are in the same segment in J if and only if they are in the same segment - in *both* S1 and S2. + The join J of S1 and S2 is defined as the segmentation in which two + voxels are in the same segment if and only if they are in the same + segment in *both* S1 and S2. Parameters ---------- @@ -35,16 +36,25 @@ def join_segmentations(s1, s2): if s1.shape != s2.shape: raise ValueError("Cannot join segmentations of different shape. " + "s1.shape: %s, s2.shape: %s" % (s1.shape, s2.shape)) - s1 = relabel_from_one(s1)[0] - s2 = relabel_from_one(s2)[0] + s1 = relabel_sequential(s1)[0] + s2 = relabel_sequential(s2)[0] j = (s2.max() + 1) * s1 + s2 - j = relabel_from_one(j)[0] + j = relabel_sequential(j)[0] return j +@deprecated('relabel_sequential') def relabel_from_one(label_field): """Convert labels in an arbitrary label field to {1, ... number_of_labels}. + This function is deprecated, see ``relabel_sequential`` for more. + """ + return relabel_sequential(label_field, offset=1) + + +def relabel_sequential(label_field, offset=1): + """Relabel arbitrary labels to {`offset`, ... `offset` + number_of_labels}. + This function also returns the forward map (mapping the original labels to the reduced labels) and the inverse map (mapping the reduced labels back to the original ones). @@ -52,6 +62,10 @@ def relabel_from_one(label_field): Parameters ---------- label_field : numpy array of int, arbitrary shape + An array of labels. + offset : int, optional + The return labels will start at `offset`, which should be + strictly positive. Returns ------- @@ -62,13 +76,15 @@ def relabel_from_one(label_field): The map from the original label space to the returned label space. Can be used to re-apply the same mapping. See examples for usage. - inverse_map : numpy array of int, shape ``(len(np.unique(label_field)),)`` + inverse_map : 1D numpy array of int, of length offset + number of labels The map from the new label space to the original space. This can be used to reconstruct the original label field from the relabeled one. Notes ----- + The label 0 is assumed to denote the background and is never remapped. + The forward map can be extremely big for some inputs, since its length is given by the maximum of the label field. However, in most situations, ``label_field.max()`` is much smaller than @@ -77,9 +93,9 @@ def relabel_from_one(label_field): Examples -------- - >>> from skimage.segmentation import relabel_from_one - >>> label_field = array([1, 1, 5, 5, 8, 99, 42]) - >>> relab, fw, inv = relabel_from_one(label_field) + >>> from skimage.segmentation import relabel_sequential + >>> label_field = np.array([1, 1, 5, 5, 8, 99, 42]) + >>> relab, fw, inv = relabel_sequential(label_field) >>> relab array([1, 1, 2, 2, 3, 5, 4]) >>> fw @@ -94,16 +110,20 @@ def relabel_from_one(label_field): True >>> (inv[relab] == label_field).all() True + >>> relab, fw, inv = relabel_sequential(label_field, offset=5) + >>> relab + array([5, 5, 6, 6, 7, 9, 8]) """ labels = np.unique(label_field) labels0 = labels[labels != 0] m = labels.max() - if m == len(labels0): # nothing to do, already 1...n labels + if m == len(labels0): # nothing to do, already 1...n labels return label_field, labels, labels forward_map = np.zeros(m+1, int) - forward_map[labels0] = np.arange(1, len(labels0) + 1) + forward_map[labels0] = np.arange(offset, offset + len(labels0) + 1) if not (labels == 0).any(): labels = np.concatenate(([0], labels)) - inverse_map = labels - return forward_map[label_field], forward_map, inverse_map - + inverse_map = np.zeros(offset - 1 + len(labels), dtype=np.intp) + inverse_map[(offset - 1):] = labels + relabeled = forward_map[label_field] + return relabeled, forward_map, inverse_map diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index 9ef3a8fc..17c4c98b 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -12,11 +12,26 @@ import warnings import numpy as np from scipy import sparse, ndimage + +# executive summary for next code block: try to import umfpack from +# scipy, but make sure not to raise a fuss if it fails since it's only +# needed to speed up a few cases. +# See discussions at: +# https://groups.google.com/d/msg/scikit-image/FrM5IGP6wh4/1hp-FtVZmfcJ +# http://stackoverflow.com/questions/13977970/ignore-exceptions-printed-to-stderr-in-del/13977992?noredirect=1#comment28386412_13977992 try: from scipy.sparse.linalg.dsolve import umfpack + old_del = umfpack.UmfpackContext.__del__ + def new_del(self): + try: + old_del(self) + except AttributeError: + pass + umfpack.UmfpackContext.__del__ = new_del UmfpackContext = umfpack.UmfpackContext() except: UmfpackContext = None + try: from pyamg import ruge_stuben_solver amg_loaded = True @@ -62,14 +77,14 @@ def _make_graph_edges_3d(n_x, n_y, n_z): return edges -def _compute_weights_3d(data, beta=130, eps=1.e-6, depth=1., +def _compute_weights_3d(data, spacing, beta=130, eps=1.e-6, multichannel=False): # Weight calculation is main difference in multispectral version # Original gradient**2 replaced with sum of gradients ** 2 gradients = 0 for channel in range(0, data.shape[-1]): gradients += _compute_gradients_3d(data[..., channel], - depth=depth) ** 2 + spacing) ** 2 # All channels considered together in this standard deviation beta /= 10 * data.std() if multichannel: @@ -82,10 +97,10 @@ def _compute_weights_3d(data, beta=130, eps=1.e-6, depth=1., return weights -def _compute_gradients_3d(data, depth=1.): - gr_deep = np.abs(data[:, :, :-1] - data[:, :, 1:]).ravel() / depth - gr_right = np.abs(data[:, :-1] - data[:, 1:]).ravel() - gr_down = np.abs(data[:-1] - data[1:]).ravel() +def _compute_gradients_3d(data, spacing): + gr_deep = np.abs(data[:, :, :-1] - data[:, :, 1:]).ravel() / spacing[2] + gr_right = np.abs(data[:, :-1] - data[:, 1:]).ravel() / spacing[1] + gr_down = np.abs(data[:-1] - data[1:]).ravel() / spacing[0] return np.r_[gr_deep, gr_right, gr_down] @@ -101,9 +116,10 @@ def _make_laplacian_sparse(edges, weights): lap = sparse.coo_matrix((data, (i_indices, j_indices)), shape=(pixel_nb, pixel_nb)) connect = - np.ravel(lap.sum(axis=1)) - lap = sparse.coo_matrix((np.hstack((data, connect)), - (np.hstack((i_indices, diag)), np.hstack((j_indices, diag)))), - shape=(pixel_nb, pixel_nb)) + lap = sparse.coo_matrix( + (np.hstack((data, connect)), (np.hstack((i_indices, diag)), + np.hstack((j_indices, diag)))), + shape=(pixel_nb, pixel_nb)) return lap.tocsr() @@ -153,14 +169,15 @@ def _mask_edges_weights(edges, weights, mask): # Reassign edges labels to 0, 1, ... edges_number - 1 order = np.searchsorted(np.unique(edges.ravel()), np.arange(max_node_index + 1)) - edges = order[edges] + edges = order[edges.astype(np.int64)] return edges, weights -def _build_laplacian(data, mask=None, beta=50, depth=1., multichannel=False): - l_x, l_y, l_z = data.shape[:3] +def _build_laplacian(data, spacing, mask=None, beta=50, + multichannel=False): + l_x, l_y, l_z = tuple(data.shape[i] for i in range(3)) edges = _make_graph_edges_3d(l_x, l_y, l_z) - weights = _compute_weights_3d(data, beta=beta, eps=1.e-10, depth=depth, + weights = _compute_weights_3d(data, spacing, beta=beta, eps=1.e-10, multichannel=multichannel) if mask is not None: edges, weights = _mask_edges_weights(edges, weights, mask) @@ -173,7 +190,8 @@ def _build_laplacian(data, mask=None, beta=50, depth=1., multichannel=False): def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, - multichannel=False, return_full_prob=False, depth=1.): + multichannel=False, return_full_prob=False, depth=1., + spacing=None): """Random walker algorithm for segmentation from markers. Random walker algorithm is implemented for gray-level or multichannel @@ -231,12 +249,16 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, return_full_prob : bool, default False If True, the probability that a pixel belongs to each of the labels will be returned, instead of only the most likely label. - depth : float, default 1. + depth : float, default 1. [DEPRECATED] Correction for non-isotropic voxel depths in 3D volumes. Default (1.) implies isotropy. This factor is derived as follows: depth = (out-of-plane voxel spacing) / (in-plane voxel spacing), where in-plane voxel spacing represents the first two spatial dimensions and out-of-plane voxel spacing represents the third spatial dimension. + `depth` is deprecated as of 0.9, in favor of `spacing`. + spacing : iterable of floats + Spacing between voxels in each spatial dimension. If `None`, then + the spacing between pixels/voxels in each dimension is assumed 1. Returns ------- @@ -259,12 +281,9 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, Multichannel inputs are scaled with all channel data combined. Ensure all channels are separately normalized prior to running this algorithm. - The `depth` argument is specifically for certain types of 3-dimensional - volumes which, due to how they were acquired, have different spacing - along in-plane and out-of-plane dimensions. This is commonly encountered - in medical imaging. The `depth` argument corrects gradients calculated - along the third spatial dimension for the otherwise inherent assumption - that all points are equally spaced. + The `spacing` argument is specifically for anisotropic datasets, where + data points are spaced differently in one or more spatial dimensions. + Anisotropic data is commonly encountered in medical imaging. The algorithm was first proposed in *Random walks for image segmentation*, Leo Grady, IEEE Trans Pattern Anal Mach Intell. @@ -324,12 +343,31 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, """ - if UmfpackContext is None: + if mode is None: + mode = 'bf' + warnings.warn("Default mode will change in the next release from 'bf' " + "to 'cg_mg' if pyamg is installed, else to 'cg' if " + "SciPy was built with UMFPACK, or to 'bf' otherwise.") + + if UmfpackContext is None and mode == 'cg': warnings.warn('SciPy was built without UMFPACK. Consider rebuilding ' 'SciPy with UMFPACK, this will greatly speed up the ' 'random walker functions. You may also install pyamg ' 'and run the random walker function in cg_mg mode ' '(see the docstrings)') + if depth != 1.: + warnings.warn('`depth` kwarg is deprecated, and will be removed in the' + ' next major release. Use `spacing` instead.') + + # Spacing kwarg checks + if spacing is None: + spacing = (1., 1.) + (depth, ) + elif len(spacing) == 2: + spacing = tuple(spacing) + (depth, ) + elif len(spacing) == 3: + pass + else: + raise ValueError('Input argument `spacing` incorrect, see docstring.') # Parse input data if not multichannel: @@ -363,10 +401,10 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, del filled labels = np.atleast_3d(labels) if np.any(labels < 0): - lap_sparse = _build_laplacian(data, mask=labels >= 0, beta=beta, - depth=depth, multichannel=multichannel) + lap_sparse = _build_laplacian(data, spacing, mask=labels >= 0, + beta=beta, multichannel=multichannel) else: - lap_sparse = _build_laplacian(data, beta=beta, depth=depth, + lap_sparse = _build_laplacian(data, spacing, beta=beta, multichannel=multichannel) lap_sparse, B = _buildAB(lap_sparse, labels) # We solve the linear system diff --git a/skimage/segmentation/tests/test_join.py b/skimage/segmentation/tests/test_join.py index f03244e9..548fcc8d 100644 --- a/skimage/segmentation/tests/test_join.py +++ b/skimage/segmentation/tests/test_join.py @@ -1,6 +1,6 @@ import numpy as np from numpy.testing import assert_array_equal, assert_raises -from skimage.segmentation import join_segmentations, relabel_from_one +from skimage.segmentation import join_segmentations, relabel_sequential def test_join_segmentations(): s1 = np.array([[0, 0, 1, 1], @@ -24,9 +24,10 @@ def test_join_segmentations(): s3 = np.array([[0, 0, 1, 1], [0, 2, 2, 1]]) assert_raises(ValueError, join_segmentations, s1, s3) -def test_relabel_from_one(): + +def test_relabel_sequential_offset1(): ar = np.array([1, 1, 5, 5, 8, 99, 42]) - ar_relab, fw, inv = relabel_from_one(ar) + ar_relab, fw, inv = relabel_sequential(ar) ar_relab_ref = np.array([1, 1, 2, 2, 3, 5, 4]) assert_array_equal(ar_relab, ar_relab_ref) fw_ref = np.zeros(100, int) @@ -36,5 +37,29 @@ def test_relabel_from_one(): assert_array_equal(inv, inv_ref) +def test_relabel_sequential_offset5(): + ar = np.array([1, 1, 5, 5, 8, 99, 42]) + ar_relab, fw, inv = relabel_sequential(ar, offset=5) + ar_relab_ref = np.array([5, 5, 6, 6, 7, 9, 8]) + assert_array_equal(ar_relab, ar_relab_ref) + fw_ref = np.zeros(100, int) + fw_ref[1] = 5; fw_ref[5] = 6; fw_ref[8] = 7; fw_ref[42] = 8; fw_ref[99] = 9 + assert_array_equal(fw, fw_ref) + inv_ref = np.array([0, 0, 0, 0, 0, 1, 5, 8, 42, 99]) + assert_array_equal(inv, inv_ref) + + +def test_relabel_sequential_offset5_with0(): + ar = np.array([1, 1, 5, 5, 8, 99, 42, 0]) + ar_relab, fw, inv = relabel_sequential(ar, offset=5) + ar_relab_ref = np.array([5, 5, 6, 6, 7, 9, 8, 0]) + assert_array_equal(ar_relab, ar_relab_ref) + fw_ref = np.zeros(100, int) + fw_ref[1] = 5; fw_ref[5] = 6; fw_ref[8] = 7; fw_ref[42] = 8; fw_ref[99] = 9 + assert_array_equal(fw, fw_ref) + inv_ref = np.array([0, 0, 0, 0, 0, 1, 5, 8, 42, 99]) + assert_array_equal(inv, inv_ref) + + if __name__ == "__main__": np.testing.run_module_suite() diff --git a/skimage/segmentation/tests/test_random_walker.py b/skimage/segmentation/tests/test_random_walker.py index 1cc0a1ee..46a82d8e 100644 --- a/skimage/segmentation/tests/test_random_walker.py +++ b/skimage/segmentation/tests/test_random_walker.py @@ -1,5 +1,6 @@ import numpy as np from skimage.segmentation import random_walker +from skimage.transform import resize def make_2d_syntheticdata(lx, ly=None): @@ -7,16 +8,16 @@ def make_2d_syntheticdata(lx, ly=None): ly = lx np.random.seed(1234) data = np.zeros((lx, ly)) + 0.1 * np.random.randn(lx, ly) - small_l = int(lx / 5) - data[lx / 2 - small_l:lx / 2 + small_l, - ly / 2 - small_l:ly / 2 + small_l] = 1 - data[lx / 2 - small_l + 1:lx / 2 + small_l - 1, - ly / 2 - small_l + 1:ly / 2 + small_l - 1] = \ - 0.1 * np.random.randn(2 * small_l - 2, 2 * small_l - 2) - data[lx / 2 - small_l, ly / 2 - small_l / 8:ly / 2 + small_l / 8] = 0 + small_l = int(lx // 5) + data[lx // 2 - small_l:lx // 2 + small_l, + ly // 2 - small_l:ly // 2 + small_l] = 1 + data[lx // 2 - small_l + 1:lx // 2 + small_l - 1, + ly // 2 - small_l + 1:ly // 2 + small_l - 1] = ( + 0.1 * np.random.randn(2 * small_l - 2, 2 * small_l - 2)) + data[lx // 2 - small_l, ly // 2 - small_l // 8:ly // 2 + small_l // 8] = 0 seeds = np.zeros_like(data) - seeds[lx / 5, ly / 5] = 1 - seeds[lx / 2 + small_l / 4, ly / 2 - small_l / 4] = 2 + seeds[lx // 5, ly // 5] = 1 + seeds[lx // 2 + small_l // 4, ly // 2 - small_l // 4] = 2 return data, seeds @@ -27,21 +28,23 @@ def make_3d_syntheticdata(lx, ly=None, lz=None): lz = lx np.random.seed(1234) data = np.zeros((lx, ly, lz)) + 0.1 * np.random.randn(lx, ly, lz) - small_l = int(lx / 5) - data[lx / 2 - small_l:lx / 2 + small_l, - ly / 2 - small_l:ly / 2 + small_l, - lz / 2 - small_l:lz / 2 + small_l] = 1 - data[lx / 2 - small_l + 1:lx / 2 + small_l - 1, - ly / 2 - small_l + 1:ly / 2 + small_l - 1, - lz / 2 - small_l + 1:lz / 2 + small_l - 1] = 0 + small_l = int(lx // 5) + data[lx // 2 - small_l:lx // 2 + small_l, + ly // 2 - small_l:ly // 2 + small_l, + lz // 2 - small_l:lz // 2 + small_l] = 1 + data[lx // 2 - small_l + 1:lx // 2 + small_l - 1, + ly // 2 - small_l + 1:ly // 2 + small_l - 1, + lz // 2 - small_l + 1:lz // 2 + small_l - 1] = 0 # make a hole - hole_size = np.max([1, small_l / 8]) - data[lx / 2 - small_l, - ly / 2 - hole_size:ly / 2 + hole_size, - lz / 2 - hole_size:lz / 2 + hole_size] = 0 + hole_size = np.max([1, small_l // 8]) + data[lx // 2 - small_l, + ly // 2 - hole_size:ly // 2 + hole_size, + lz // 2 - hole_size:lz // 2 + hole_size] = 0 seeds = np.zeros_like(data) - seeds[lx / 5, ly / 5, lz / 5] = 1 - seeds[lx / 2 + small_l / 4, ly / 2 - small_l / 4, lz / 2 - small_l / 4] = 2 + seeds[lx // 5, ly // 5, lz // 5] = 1 + seeds[lx // 2 + small_l // 4, + ly // 2 - small_l // 4, + lz // 2 - small_l // 4] = 2 return data, seeds @@ -101,7 +104,7 @@ def test_types(): lx = 70 ly = 100 data, labels = make_2d_syntheticdata(lx, ly) - data = 255 * (data - data.min()) / (data.max() - data.min()) + data = 255 * (data - data.min()) // (data.max() - data.min()) data = data.astype(np.uint8) labels_cg_mg = random_walker(data, labels, beta=90, mode='cg_mg') assert (labels_cg_mg[25:45, 40:60] == 2).all() @@ -181,6 +184,77 @@ def test_multispectral_3d(): return data, multi_labels, single_labels, labels +def test_depth(): + n = 30 + lx, ly, lz = n, n, n + data, _ = make_3d_syntheticdata(lx, ly, lz) + + # Rescale `data` along Z axis + data_aniso = np.zeros((n, n, n // 2)) + for i, yz in enumerate(data): + data_aniso[i, :, :] = resize(yz, (n, n // 2)) + + # Generate new labels + small_l = int(lx // 5) + labels_aniso = np.zeros_like(data_aniso) + labels_aniso[lx // 5, ly // 5, lz // 5] = 1 + labels_aniso[lx // 2 + small_l // 4, + ly // 2 - small_l // 4, + lz // 4 - small_l // 8] = 2 + + # Test with `depth` kwarg + labels_aniso = random_walker(data_aniso, labels_aniso, mode='cg', + depth=0.5) + + assert (labels_aniso[13:17, 13:17, 7:9] == 2).all() + + +def test_spacing(): + n = 30 + lx, ly, lz = n, n, n + data, _ = make_3d_syntheticdata(lx, ly, lz) + + # Rescale `data` along Y axis + # `resize` is not yet 3D capable, so this must be done by looping in 2D. + data_aniso = np.zeros((n, n * 2, n)) + for i, yz in enumerate(data): + data_aniso[i, :, :] = resize(yz, (n * 2, n)) + + # Generate new labels + small_l = int(lx // 5) + labels_aniso = np.zeros_like(data_aniso) + labels_aniso[lx // 5, ly // 5, lz // 5] = 1 + labels_aniso[lx // 2 + small_l // 4, + ly - small_l // 2, + lz // 2 - small_l // 4] = 2 + + # Test with `spacing` kwarg + # First, anisotropic along Y + labels_aniso = random_walker(data_aniso, labels_aniso, mode='cg', + spacing=(1., 2., 1.)) + assert (labels_aniso[13:17, 26:34, 13:17] == 2).all() + + # Rescale `data` along X axis + # `resize` is not yet 3D capable, so this must be done by looping in 2D. + data_aniso = np.zeros((n, n * 2, n)) + for i in range(data.shape[1]): + data_aniso[i, :, :] = resize(data[:, 1, :], (n * 2, n)) + + # Generate new labels + small_l = int(lx // 5) + labels_aniso2 = np.zeros_like(data_aniso) + labels_aniso2[lx // 5, ly // 5, lz // 5] = 1 + labels_aniso2[lx - small_l // 2, + ly // 2 + small_l // 4, + lz // 2 - small_l // 4] = 2 + + # Anisotropic along X + labels_aniso2 = random_walker(data_aniso, + labels_aniso2, + mode='cg', spacing=(2., 1., 1.)) + assert (labels_aniso2[26:34, 13:17, 13:17] == 2).all() + + if __name__ == '__main__': from numpy import testing testing.run_module_suite() diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index 8fa2cfdb..4f00a076 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -1,7 +1,6 @@ from ._hough_transform import (hough_circle, hough_ellipse, hough_line, probabilistic_hough_line) -from .hough_transform import (hough, probabilistic_hough, hough_peaks, - hough_line_peaks) +from .hough_transform import hough_line_peaks from .radon_transform import radon, iradon, iradon_sart from .finite_radon_transform import frt2, ifrt2 from .integral import integral_image, integrate @@ -18,7 +17,6 @@ __all__ = ['hough_circle', 'hough_ellipse', 'hough_line', 'probabilistic_hough_line', - 'hough', 'probabilistic_hough', 'hough_peaks', 'hough_line_peaks', diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 25a13765..c1f49fa5 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -341,7 +341,7 @@ class AffineTransform(ProjectiveTransform): return self._matrix[0:2, 2] -class PiecewiseAffineTransform(ProjectiveTransform): +class PiecewiseAffineTransform(GeometricTransform): """2D piecewise affine transformation. @@ -1031,21 +1031,24 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, out = None - # use fast Cython version for specific interpolation orders + # use fast Cython version for specific interpolation orders and input if order in range(4) and not map_args: + matrix = None + # inverse_map is a transformation matrix as numpy array if isinstance(inverse_map, np.ndarray) and inverse_map.shape == (3, 3): matrix = inverse_map - elif inverse_map in HOMOGRAPHY_TRANSFORMS: + # inverse_map is a homography + elif isinstance(inverse_map, HOMOGRAPHY_TRANSFORMS): matrix = inverse_map._matrix + # inverse_map is the inverse of a homography elif (hasattr(inverse_map, '__name__') and inverse_map.__name__ == 'inverse' - and get_bound_method_class(inverse_map) - in HOMOGRAPHY_TRANSFORMS): - + and isinstance(get_bound_method_class(inverse_map), + HOMOGRAPHY_TRANSFORMS)): matrix = np.linalg.inv(six.get_method_self(inverse_map)._matrix) if matrix is not None: @@ -1067,6 +1070,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, rows, cols = output_shape[:2] + # inverse_map is a transformation matrix as numpy array if isinstance(inverse_map, np.ndarray) and inverse_map.shape == (3, 3): inverse_map = ProjectiveTransform(matrix=inverse_map) @@ -1075,19 +1079,19 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, coords = warp_coords(coord_map, (rows, cols, bands)) - # Prefilter not necessary for order 0, 1 interpolation + # Pre-filtering not necessary for order 0, 1 interpolation prefilter = order > 1 out = ndimage.map_coordinates(image, coords, prefilter=prefilter, mode=mode, order=order, cval=cval) - # The spline filters sometimes return results outside [0, 1], - # so clip to ensure valid data - clipped = np.clip(out, 0, 1) + # The spline filters sometimes return results outside [0, 1], + # so clip to ensure valid data + clipped = np.clip(out, 0, 1) - if mode == 'constant' and not (0 <= cval <= 1): - clipped[out == cval] = cval + if mode == 'constant' and not (0 <= cval <= 1): + clipped[out == cval] = cval - out = clipped + out = clipped if out.ndim == 3 and orig_ndim == 2: # remove singleton dimension introduced by atleast_3d diff --git a/skimage/transform/_hough_transform.pyx b/skimage/transform/_hough_transform.pyx index a3aa2649..29344fa8 100644 --- a/skimage/transform/_hough_transform.pyx +++ b/skimage/transform/_hough_transform.pyx @@ -7,7 +7,7 @@ import numpy as np cimport numpy as cnp cimport cython -from libc.math cimport abs, fabs, sqrt, ceil +from libc.math cimport abs, fabs, sqrt, ceil, atan2, M_PI from libc.stdlib cimport rand from skimage.draw import circle_perimeter @@ -122,17 +122,18 @@ def hough_ellipse(cnp.ndarray img, int threshold=4, double accuracy=1, Returns ------- - res : list of tuples [(x0, y0, a, b, angle, accumulator)] - Where (x0, y0) is the center, (a, b) major and minor axis. - The angle value follows `draw.ellipse_perimeter()` convention. + result : ndarray with fields [(accumulator, y0, x0, a, b, orientation)] + Where ``(yc, xc)`` is the center, ``(a, b)`` the major and minor + axes, respectively. The `orientation` value follows + `skimage.draw.ellipse_perimeter` convention. Examples -------- - >>> img = np.zeros((25, 25), dtype=int) - >>> rr, cc = draw.ellipse_perimeter(10, 10, 6, 8) + >>> img = np.zeros((25, 25), dtype=np.uint8) + >>> rr, cc = ellipse_perimeter(10, 10, 6, 8) >>> img[cc, rr] = 1 >>> result = hough_ellipse(img, threshold=8) - [(10.0, 10.0, 8.0, 6.0, 0.0, 10)] + [(10, 10.0, 8.0, 6.0, 0.0, 10.0)] Notes ----- @@ -149,47 +150,47 @@ def hough_ellipse(cnp.ndarray img, int threshold=4, double accuracy=1, if img.ndim != 2: raise ValueError('The input image must be 2D.') - cdef long[:, :] pixels = np.transpose(np.nonzero(img)) - cdef Py_ssize_t num_pixels = pixels.shape[0] + cdef Py_ssize_t[:, ::1] pixels = np.row_stack(np.nonzero(img)) + cdef Py_ssize_t num_pixels = pixels.shape[1] cdef list acc = list() cdef list results = list() - cdef bin_size = accuracy**2 + cdef double bin_size = accuracy ** 2 cdef int max_b_squared if max_size is None: if img.shape[0] < img.shape[1]: - max_b_squared = np.round(0.5 * img.shape[0])**2 + max_b_squared = np.round(0.5 * img.shape[0]) ** 2 else: - max_b_squared = np.round(0.5 * img.shape[1])**2 + max_b_squared = np.round(0.5 * img.shape[1]) ** 2 else: max_b_squared = max_size**2 cdef Py_ssize_t p1, p2, p3, p1x, p1y, p2x, p2y, p3x, p3y - cdef double x0, y0, a, b, d, k - cdef double cos_tau_squared, b_squared, f_squared, angle + cdef double xc, yc, a, b, d, k + cdef double cos_tau_squared, b_squared, f_squared, orientation for p1 in range(num_pixels): - p1x = pixels[p1, 1] - p1y = pixels[p1, 0] + p1x = pixels[1, p1] + p1y = pixels[0, p1] for p2 in range(p1): - p2x = pixels[p2, 1] - p2y = pixels[p2, 0] + p2x = pixels[1, p2] + p2y = pixels[0, p2] - # Candidate: center (x0, y0) and main axis a + # Candidate: center (xc, yc) and main axis a a = 0.5 * sqrt((p1x - p2x)**2 + (p1y - p2y)**2) if a > 0.5 * min_size: - x0 = 0.5 * (p1x + p2x) - y0 = 0.5 * (p1y + p2y) + xc = 0.5 * (p1x + p2x) + yc = 0.5 * (p1y + p2y) for p3 in range(num_pixels): - p3x = pixels[p3, 1] - p3y = pixels[p3, 0] + p3x = pixels[1, p3] + p3y = pixels[0, p3] - d = sqrt((p3x - x0)**2 + (p3y - y0)**2) + d = sqrt((p3x - xc)**2 + (p3y - yc)**2) if d > min_size: f_squared = (p3x - p1x)**2 + (p3y - p1y)**2 - cos_tau_squared = ((a**2 + d**2 - f_squared) \ + cos_tau_squared = ((a**2 + d**2 - f_squared) / (2 * a * d))**2 # Consider b2 > 0 and avoid division by zero k = a**2 - d**2 * cos_tau_squared @@ -205,21 +206,29 @@ def hough_ellipse(cnp.ndarray img, int threshold=4, double accuracy=1, hist, bin_edges = np.histogram(acc, bins=bins) hist_max = np.max(hist) if hist_max > threshold: - angle = np.arctan2(p1x - p2x, p1y - p2y) - # pi - angle to keep ellipse_perimeter() convention - if angle != 0: - angle = np.pi - angle + orientation = atan2(p1x - p2x, p1y - p2y) b = sqrt(bin_edges[hist.argmax()]) - results.append((x0, - y0, - a, - b, - angle, - hist_max, # Accumulator - )) + # to keep ellipse_perimeter() convention + if orientation != 0: + orientation = M_PI - orientation + # When orientation is not in [-pi:pi] + # it would mean in ellipse_perimeter() + # that a < b. But we keep a > b. + if orientation > M_PI: + orientation = orientation - M_PI / 2. + a, b = b, a + results.append((hist_max, # Accumulator + yc, xc, + a, b, + orientation)) acc = [] - return results + return np.array(results, dtype=[('accumulator', np.intp), + ('yc', np.double), + ('xc', np.double), + ('a', np.double), + ('b', np.double), + ('orientation', np.double)]) def hough_line(cnp.ndarray img, diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx index b3136c22..433c586d 100644 --- a/skimage/transform/_warps_cy.pyx +++ b/skimage/transform/_warps_cy.pyx @@ -89,11 +89,11 @@ def _warp_fast(cnp.ndarray image, cnp.ndarray H, output_shape=None, cdef Py_ssize_t out_r, out_c if output_shape is None: - out_r = img.shape[0] - out_c = img.shape[1] + out_r = int(img.shape[0]) + out_c = int(img.shape[1]) else: - out_r = output_shape[0] - out_c = output_shape[1] + out_r = int(output_shape[0]) + out_c = int(output_shape[1]) cdef double[:, ::1] out = np.zeros((out_r, out_c), dtype=np.double) diff --git a/skimage/transform/hough_transform.py b/skimage/transform/hough_transform.py index 15968fd8..0a7d35a2 100644 --- a/skimage/transform/hough_transform.py +++ b/skimage/transform/hough_transform.py @@ -3,30 +3,6 @@ from scipy import ndimage from skimage import measure, morphology -from ._hough_transform import hough_line, probabilistic_hough_line -from skimage._shared.utils import deprecated - - -@deprecated('hough_line') -def hough(img, theta=None): - return hough_line(img, theta) - - -@deprecated('probabilistic_hough_line') -def probabilistic_hough(img, threshold=10, line_length=50, line_gap=10, - theta=None): - return probabilistic_hough_line(img, threshold=threshold, - line_length=line_length, line_gap=line_gap, - theta=theta) - - -@deprecated('hough_line_peaks') -def hough_peaks(hspace, angles, dists, min_distance=10, min_angle=10, - threshold=None, num_peaks=np.inf): - return hough_line_peaks(hspace, angles, dists, min_distance, min_angle, - threshold, num_peaks) - - def hough_line_peaks(hspace, angles, dists, min_distance=9, min_angle=10, threshold=None, num_peaks=np.inf): """Return peaks in hough transform. diff --git a/skimage/transform/tests/test_hough_transform.py b/skimage/transform/tests/test_hough_transform.py index 148f4d43..fb19d8c1 100644 --- a/skimage/transform/tests/test_hough_transform.py +++ b/skimage/transform/tests/test_hough_transform.py @@ -1,7 +1,5 @@ import numpy as np -from numpy.testing import (assert_almost_equal, - assert_equal, - ) +from numpy.testing import assert_almost_equal, assert_equal import skimage.transform as tf from skimage.draw import line, circle_perimeter, ellipse_perimeter @@ -81,8 +79,10 @@ def test_hough_line_peaks_dist(): img[:, 30] = True img[:, 40] = True hspace, angles, dists = tf.hough_line(img) - assert len(tf.hough_line_peaks(hspace, angles, dists, min_distance=5)[0]) == 2 - assert len(tf.hough_line_peaks(hspace, angles, dists, min_distance=15)[0]) == 1 + assert len(tf.hough_line_peaks(hspace, angles, dists, + min_distance=5)[0]) == 2 + assert len(tf.hough_line_peaks(hspace, angles, dists, + min_distance=15)[0]) == 1 def test_hough_line_peaks_angle(): @@ -91,18 +91,24 @@ def test_hough_line_peaks_angle(): img[0, :] = True hspace, angles, dists = tf.hough_line(img) - assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=45)[0]) == 2 - assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=90)[0]) == 1 + assert len(tf.hough_line_peaks(hspace, angles, dists, + min_angle=45)[0]) == 2 + assert len(tf.hough_line_peaks(hspace, angles, dists, + min_angle=90)[0]) == 1 theta = np.linspace(0, np.pi, 100) hspace, angles, dists = tf.hough_line(img, theta) - assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=45)[0]) == 2 - assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=90)[0]) == 1 + assert len(tf.hough_line_peaks(hspace, angles, dists, + min_angle=45)[0]) == 2 + assert len(tf.hough_line_peaks(hspace, angles, dists, + min_angle=90)[0]) == 1 theta = np.linspace(np.pi / 3, 4. / 3 * np.pi, 100) hspace, angles, dists = tf.hough_line(img, theta) - assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=45)[0]) == 2 - assert len(tf.hough_line_peaks(hspace, angles, dists, min_angle=90)[0]) == 1 + assert len(tf.hough_line_peaks(hspace, angles, dists, + min_angle=45)[0]) == 2 + assert len(tf.hough_line_peaks(hspace, angles, dists, + min_angle=90)[0]) == 1 def test_hough_line_peaks_num(): @@ -122,7 +128,7 @@ def test_hough_circle(): y, x = circle_perimeter(y_0, x_0, radius) img[x, y] = 1 - out = tf.hough_circle(img, np.array([radius])) + out = tf.hough_circle(img, np.array([radius], dtype=np.intp)) x, y = np.where(out[0] == out[0].max()) assert_equal(x[0], x_0) @@ -138,7 +144,8 @@ def test_hough_circle_extended(): y, x = circle_perimeter(y_0, x_0, radius) img[x[np.where(x > 0)], y[np.where(x > 0)]] = 1 - out = tf.hough_circle(img, np.array([radius]), full_output=True) + out = tf.hough_circle(img, np.array([radius], dtype=np.intp), + full_output=True) x, y = np.where(out[0] == out[0].max()) # Offset for x_0, y_0 @@ -148,36 +155,204 @@ def test_hough_circle_extended(): def test_hough_ellipse_zero_angle(): img = np.zeros((25, 25), dtype=int) - a = 6 - b = 8 + rx = 6 + ry = 8 x0 = 12 - y0 = 12 + y0 = 15 angle = 0 - rr, cc = ellipse_perimeter(x0, x0, b, a) + rr, cc = ellipse_perimeter(y0, x0, ry, rx) img[rr, cc] = 1 result = tf.hough_ellipse(img, threshold=9) - assert_equal(result[0][0], x0) - assert_equal(result[0][1], y0) - assert_almost_equal(result[0][2], b, decimal=1) - assert_almost_equal(result[0][3], a, decimal=1) - assert_equal(result[0][4], angle) + best = result[-1] + assert_equal(best[1], y0) + assert_equal(best[2], x0) + assert_almost_equal(best[3], ry, decimal=1) + assert_almost_equal(best[4], rx, decimal=1) + assert_equal(best[5], angle) + # Check if I re-draw the ellipse, points are the same! + # ie check API compatibility between hough_ellipse and ellipse_perimeter + rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), + orientation=best[5]) + assert_equal(rr, rr2) + assert_equal(cc, cc2) -def test_hough_ellipse_non_zero_angle(): - img = np.zeros((20, 20), dtype=int) - a = 6 - b = 9 +def test_hough_ellipse_non_zero_posangle1(): + # ry > rx, angle in [0:pi/2] + img = np.zeros((30, 24), dtype=int) + rx = 6 + ry = 12 x0 = 10 - y0 = 10 + y0 = 15 angle = np.pi / 1.35 - rr, cc = ellipse_perimeter(x0, x0, b, a, orientation=angle) + rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle) img[rr, cc] = 1 result = tf.hough_ellipse(img, threshold=15, accuracy=3) - assert_almost_equal(result[0][0] / 100., x0 / 100., decimal=1) - assert_almost_equal(result[0][1] / 100., y0 / 100., decimal=1) - assert_almost_equal(result[0][2] / 100., b / 100., decimal=1) - assert_almost_equal(result[0][3] / 100., a / 100., decimal=1) - assert_almost_equal(result[0][4], angle, decimal=1) + result.sort(order='accumulator') + best = result[-1] + assert_almost_equal(best[1] / 100., y0 / 100., decimal=1) + assert_almost_equal(best[2] / 100., x0 / 100., decimal=1) + assert_almost_equal(best[3] / 10., ry / 10., decimal=1) + assert_almost_equal(best[4] / 100., rx / 100., decimal=1) + assert_almost_equal(best[5], angle, decimal=1) + # Check if I re-draw the ellipse, points are the same! + # ie check API compatibility between hough_ellipse and ellipse_perimeter + rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), + orientation=best[5]) + assert_equal(rr, rr2) + assert_equal(cc, cc2) + + +def test_hough_ellipse_non_zero_posangle2(): + # ry < rx, angle in [0:pi/2] + img = np.zeros((30, 24), dtype=int) + rx = 12 + ry = 6 + x0 = 10 + y0 = 15 + angle = np.pi / 1.35 + rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle) + img[rr, cc] = 1 + result = tf.hough_ellipse(img, threshold=15, accuracy=3) + result.sort(order='accumulator') + best = result[-1] + assert_almost_equal(best[1] / 100., y0 / 100., decimal=1) + assert_almost_equal(best[2] / 100., x0 / 100., decimal=1) + assert_almost_equal(best[3] / 10., ry / 10., decimal=1) + assert_almost_equal(best[4] / 100., rx / 100., decimal=1) + assert_almost_equal(best[5], angle, decimal=1) + # Check if I re-draw the ellipse, points are the same! + # ie check API compatibility between hough_ellipse and ellipse_perimeter + rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), + orientation=best[5]) + assert_equal(rr, rr2) + assert_equal(cc, cc2) + + +def test_hough_ellipse_non_zero_posangle3(): + # ry < rx, angle in [pi/2:pi] + img = np.zeros((30, 24), dtype=int) + rx = 12 + ry = 6 + x0 = 10 + y0 = 15 + angle = np.pi / 1.35 + np.pi / 2. + rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle) + img[rr, cc] = 1 + result = tf.hough_ellipse(img, threshold=15, accuracy=3) + result.sort(order='accumulator') + best = result[-1] + # Check if I re-draw the ellipse, points are the same! + # ie check API compatibility between hough_ellipse and ellipse_perimeter + rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), + orientation=best[5]) + assert_equal(rr, rr2) + assert_equal(cc, cc2) + + +def test_hough_ellipse_non_zero_posangle4(): + # ry < rx, angle in [pi:3pi/4] + img = np.zeros((30, 24), dtype=int) + rx = 12 + ry = 6 + x0 = 10 + y0 = 15 + angle = np.pi / 1.35 + np.pi + rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle) + img[rr, cc] = 1 + result = tf.hough_ellipse(img, threshold=15, accuracy=3) + result.sort(order='accumulator') + best = result[-1] + # Check if I re-draw the ellipse, points are the same! + # ie check API compatibility between hough_ellipse and ellipse_perimeter + rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), + orientation=best[5]) + assert_equal(rr, rr2) + assert_equal(cc, cc2) + + +def test_hough_ellipse_non_zero_negangle1(): + # ry > rx, angle in [0:-pi/2] + img = np.zeros((30, 24), dtype=int) + rx = 6 + ry = 12 + x0 = 10 + y0 = 15 + angle = - np.pi / 1.35 + rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle) + img[rr, cc] = 1 + result = tf.hough_ellipse(img, threshold=15, accuracy=3) + result.sort(order='accumulator') + best = result[-1] + # Check if I re-draw the ellipse, points are the same! + # ie check API compatibility between hough_ellipse and ellipse_perimeter + rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), + orientation=best[5]) + assert_equal(rr, rr2) + assert_equal(cc, cc2) + + +def test_hough_ellipse_non_zero_negangle2(): + # ry < rx, angle in [0:-pi/2] + img = np.zeros((30, 24), dtype=int) + rx = 12 + ry = 6 + x0 = 10 + y0 = 15 + angle = - np.pi / 1.35 + rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle) + img[rr, cc] = 1 + result = tf.hough_ellipse(img, threshold=15, accuracy=3) + result.sort(order='accumulator') + best = result[-1] + # Check if I re-draw the ellipse, points are the same! + # ie check API compatibility between hough_ellipse and ellipse_perimeter + rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), + orientation=best[5]) + assert_equal(rr, rr2) + assert_equal(cc, cc2) + + +def test_hough_ellipse_non_zero_negangle3(): + # ry < rx, angle in [-pi/2:-pi] + img = np.zeros((30, 24), dtype=int) + rx = 12 + ry = 6 + x0 = 10 + y0 = 15 + angle = - np.pi / 1.35 - np.pi / 2. + rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle) + img[rr, cc] = 1 + result = tf.hough_ellipse(img, threshold=15, accuracy=3) + result.sort(order='accumulator') + best = result[-1] + # Check if I re-draw the ellipse, points are the same! + # ie check API compatibility between hough_ellipse and ellipse_perimeter + rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), + orientation=best[5]) + assert_equal(rr, rr2) + assert_equal(cc, cc2) + + +def test_hough_ellipse_non_zero_negangle4(): + # ry < rx, angle in [-pi:-3pi/4] + img = np.zeros((30, 24), dtype=int) + rx = 12 + ry = 6 + x0 = 10 + y0 = 15 + angle = - np.pi / 1.35 - np.pi + rr, cc = ellipse_perimeter(y0, x0, ry, rx, orientation=angle) + img[rr, cc] = 1 + result = tf.hough_ellipse(img, threshold=15, accuracy=3) + result.sort(order='accumulator') + best = result[-1] + # Check if I re-draw the ellipse, points are the same! + # ie check API compatibility between hough_ellipse and ellipse_perimeter + rr2, cc2 = ellipse_perimeter(y0, x0, int(best[3]), int(best[4]), + orientation=best[5]) + assert_equal(rr, rr2) + assert_equal(cc, cc2) if __name__ == "__main__": diff --git a/skimage/util/noise.py b/skimage/util/noise.py index d0a3c5f6..9283f537 100644 --- a/skimage/util/noise.py +++ b/skimage/util/noise.py @@ -5,7 +5,7 @@ from .dtype import img_as_float __all__ = ['random_noise'] -def random_noise(image, mode='gaussian', seed=None, **kwargs): +def random_noise(image, mode='gaussian', seed=None, clip=True, **kwargs): """ Function to add random noise of various types to a floating-point image. @@ -17,6 +17,8 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs): One of the following strings, selecting the type of noise to add: 'gaussian' Gaussian-distributed additive noise. + 'localvar' Gaussian-distributed additive noise, with specified + local variance at each point of `image` 'poisson' Poisson-distributed noise generated from the data. 'salt' Replaces random pixels with 1. 'pepper' Replaces random pixels with 0. @@ -26,12 +28,20 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs): seed : int If provided, this will set the random seed before generating noise, for valid pseudo-random comparisons. + clip : bool + If True (default), the output will be clipped after noise applied + for modes `'speckle'`, `'poisson'`, and `'gaussian'`. This is + needed to maintain the proper image data range. If False, clipping + is not applied, and the output may extend beyond the range [-1, 1]. mean : float Mean of random distribution. Used in 'gaussian' and 'speckle'. Default : 0. var : float Variance of random distribution. Used in 'gaussian' and 'speckle'. Note: variance = (standard deviation) ** 2. Default : 0.01 + local_vars : ndarray + Array of positive floats, same shape as `image`, defining the local + variance at every image point. Used in 'localvar'. amount : float Proportion of image pixels to replace with noise on range [0, 1]. Used in 'salt', 'pepper', and 'salt & pepper'. Default : 0.05 @@ -42,17 +52,51 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs): Returns ------- out : ndarray - Output floating-point image data on range [0, 1]. + Output floating-point image data on range [0, 1] or [-1, 1] if the + input `image` was unsigned or signed, respectively. + + Notes + ----- + Speckle, Poisson, Localvar, and Gaussian noise may generate noise outside + the valid image range. The default is to clip (not alias) these values, + but they may be preserved by setting `clip=False`. Note that in this case + the output may contain values outside the ranges [0, 1] or [-1, 1]. + Use this option with care. + + Because of the prevalence of exclusively positive floating-point images in + intermediate calculations, it is not possible to intuit if an input is + signed based on dtype alone. Instead, negative values are explicity + searched for. Only if found does this function assume signed input. + Unexpected results only occur in rare, poorly exposes cases (e.g. if all + values are above 50 percent gray in a signed `image`). In this event, + manually scaling the input to the positive domain will solve the problem. + + The Poisson distribution is only defined for positive integers. To apply + this noise type, the number of unique values in the image is found and + the next round power of two is used to scale up the floating-point result, + after which it is scaled back down to the floating-point image range. + + To generate Poisson noise against a signed image, the signed image is + temporarily converted to an unsigned image in the floating point domain, + Poisson noise is generated, then it is returned to the original range. """ mode = mode.lower() + + # Detect if a signed image was input + if image.min() < 0: + low_clip = -1. + else: + low_clip = 0. + image = img_as_float(image) if seed is not None: np.random.seed(seed=seed) allowedtypes = { 'gaussian': 'gaussian_values', - 'poisson': '', + 'localvar': 'localvar_values', + 'poisson': 'poisson_values', 'salt': 'sp_values', 'pepper': 'sp_values', 's&p': 's&p_values', @@ -62,12 +106,15 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs): 'mean': 0., 'var': 0.01, 'amount': 0.05, - 'salt_vs_pepper': 0.5} + 'salt_vs_pepper': 0.5, + 'local_vars': np.zeros_like(image) + 0.01} allowedkwargs = { 'gaussian_values': ['mean', 'var'], + 'localvar_values': ['local_vars'], 'sp_values': ['amount'], - 's&p_values': ['amount', 'salt_vs_pepper']} + 's&p_values': ['amount', 'salt_vs_pepper'], + 'poisson_values': []} for key in kwargs: if key not in allowedkwargs[allowedtypes[mode]]: @@ -81,16 +128,32 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs): if mode == 'gaussian': noise = np.random.normal(kwargs['mean'], kwargs['var'] ** 0.5, image.shape) - out = np.clip(image + noise, 0., 1.) + out = image + noise + + elif mode == 'localvar': + # Ensure local variance input is correct + if (kwargs['local_vars'] <= 0).any(): + raise ValueError('All values of `local_vars` must be > 0.') + + # Safe shortcut usage broadcasts kwargs['local_vars'] as a ufunc + out = image + np.random.normal(0, kwargs['local_vars'] ** 0.5) elif mode == 'poisson': + # Determine unique values in image & calculate the next power of two + vals = len(np.unique(image)) + vals = 2 ** np.ceil(np.log2(vals)) + + # Ensure image is exclusively positive + if low_clip == -1.: + old_max = image.max() + image = (image + 1.) / (old_max + 1.) + # Generating noise for each unique value in image. - out = np.zeros_like(image) - for val in np.unique(image): - # Generate mask for a unique value, replace w/values drawn from - # Poisson distribution about the unique value - mask = image == val - out[mask] = np.poisson(val, mask.sum()) + out = np.random.poisson(image * vals) / float(vals) + + # Return image to original range if input was signed + if low_clip == -1.: + out = out * (old_max + 1.) - 1. elif mode == 'salt': # Re-call function with mode='s&p' and p=1 (all salt noise) @@ -119,11 +182,15 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs): kwargs['amount'] * image.size * (1. - kwargs['salt_vs_pepper'])) coords = [np.random.randint(0, i - 1, int(num_pepper)) for i in image.shape] - out[coords] = 0 + out[coords] = low_clip elif mode == 'speckle': noise = np.random.normal(kwargs['mean'], kwargs['var'] ** 0.5, image.shape) - out = np.clip(image + image * noise, 0., 1.) + out = image + image * noise + + # Clip back to original range, if necessary + if clip: + out = np.clip(out, low_clip, 1.0) return out diff --git a/skimage/util/shape.py b/skimage/util/shape.py index 5fe27a36..f91286c3 100644 --- a/skimage/util/shape.py +++ b/skimage/util/shape.py @@ -233,9 +233,7 @@ def view_as_windows(arr_in, window_shape, step=1): tuple(window_shape) arr_strides = np.array(arr_in.strides) - new_strides = np.concatenate( - (arr_strides * step, arr_strides) - ) + new_strides = np.concatenate((arr_strides * step, arr_strides)) arr_out = as_strided(arr_in, shape=new_shape, strides=new_strides) diff --git a/skimage/util/tests/test_random_noise.py b/skimage/util/tests/test_random_noise.py index 87005d60..39477cb0 100644 --- a/skimage/util/tests/test_random_noise.py +++ b/skimage/util/tests/test_random_noise.py @@ -23,12 +23,14 @@ def test_salt(): # Ensure approximately correct amount of noise was added proportion = float(saltmask.sum()) / (cam.shape[0] * cam.shape[1]) - assert 0.11 < proportion <= 0.18 + assert 0.11 < proportion <= 0.15 def test_pepper(): seed = 42 cam = img_as_float(camera()) + data_signed = cam * 2. - 1. # Same image, on range [-1, 1] + cam_noisy = random_noise(cam, seed=seed, mode='pepper', amount=0.15) peppermask = cam != cam_noisy @@ -37,7 +39,16 @@ def test_pepper(): # Ensure approximately correct amount of noise was added proportion = float(peppermask.sum()) / (cam.shape[0] * cam.shape[1]) - assert 0.11 < proportion <= 0.18 + assert 0.11 < proportion <= 0.15 + + # Check to make sure pepper gets added properly to signed images + orig_zeros = (data_signed == -1).sum() + cam_noisy_signed = random_noise(data_signed, seed=seed, mode='pepper', + amount=.15) + + proportion = (float((cam_noisy_signed == -1).sum() - orig_zeros) / + (cam.shape[0] * cam.shape[1])) + assert 0.11 < proportion <= 0.15 def test_salt_and_pepper(): @@ -72,10 +83,35 @@ def test_gaussian(): assert 0.012 < data_gaussian.var() < 0.018 +def test_localvar(): + seed = 42 + data = np.zeros((128, 128)) + 0.5 + local_vars = np.zeros((128, 128)) + 0.001 + local_vars[:64, 64:] = 0.1 + local_vars[64:, :64] = 0.25 + local_vars[64:, 64:] = 0.45 + + data_gaussian = random_noise(data, mode='localvar', seed=seed, + local_vars=local_vars, clip=False) + assert 0. < data_gaussian[:64, :64].var() < 0.002 + assert 0.095 < data_gaussian[:64, 64:].var() < 0.105 + assert 0.245 < data_gaussian[64:, :64].var() < 0.255 + assert 0.445 < data_gaussian[64:, 64:].var() < 0.455 + + # Ensure local variance bounds checking works properly + bad_local_vars = np.zeros_like(data) + assert_raises(ValueError, random_noise, data, mode='localvar', seed=seed, + local_vars=bad_local_vars) + bad_local_vars += 0.1 + bad_local_vars[0, 0] = -1 + assert_raises(ValueError, random_noise, data, mode='localvar', seed=seed, + local_vars=bad_local_vars) + + def test_speckle(): seed = 42 data = np.zeros((128, 128)) + 0.1 - np.random.seed(seed=42) + np.random.seed(seed=seed) noise = np.random.normal(0.1, 0.02 ** 0.5, (128, 128)) expected = np.clip(data + data * noise, 0, 1) @@ -84,6 +120,78 @@ def test_speckle(): assert_allclose(expected, data_speckle) +def test_poisson(): + seed = 42 + data = camera() # 512x512 grayscale uint8 + cam_noisy = random_noise(data, mode='poisson', seed=seed) + cam_noisy2 = random_noise(data, mode='poisson', seed=seed, clip=False) + + np.random.seed(seed=seed) + expected = np.random.poisson(img_as_float(data) * 256) / 256. + assert_allclose(cam_noisy, np.clip(expected, 0., 1.)) + assert_allclose(cam_noisy2, expected) + + +def test_clip_poisson(): + seed = 42 + data = camera() # 512x512 grayscale uint8 + data_signed = img_as_float(data) * 2. - 1. # Same image, on range [-1, 1] + + # Signed and unsigned, clipped + cam_poisson = random_noise(data, mode='poisson', seed=seed, clip=True) + cam_poisson2 = random_noise(data_signed, mode='poisson', seed=seed, + clip=True) + assert (cam_poisson.max() == 1.) and (cam_poisson.min() == 0.) + assert (cam_poisson2.max() == 1.) and (cam_poisson2.min() == -1.) + + # Signed and unsigned, unclipped + cam_poisson = random_noise(data, mode='poisson', seed=seed, clip=False) + cam_poisson2 = random_noise(data_signed, mode='poisson', seed=seed, + clip=False) + assert (cam_poisson.max() > 1.15) and (cam_poisson.min() == 0.) + assert (cam_poisson2.max() > 1.3) and (cam_poisson2.min() == -1.) + + +def test_clip_gaussian(): + seed = 42 + data = camera() # 512x512 grayscale uint8 + data_signed = img_as_float(data) * 2. - 1. # Same image, on range [-1, 1] + + # Signed and unsigned, clipped + cam_gauss = random_noise(data, mode='gaussian', seed=seed, clip=True) + cam_gauss2 = random_noise(data_signed, mode='gaussian', seed=seed, + clip=True) + assert (cam_gauss.max() == 1.) and (cam_gauss.min() == 0.) + assert (cam_gauss2.max() == 1.) and (cam_gauss2.min() == -1.) + + # Signed and unsigned, unclipped + cam_gauss = random_noise(data, mode='gaussian', seed=seed, clip=False) + cam_gauss2 = random_noise(data_signed, mode='gaussian', seed=seed, + clip=False) + assert (cam_gauss.max() > 1.22) and (cam_gauss.min() < -0.36) + assert (cam_gauss2.max() > 1.219) and (cam_gauss2.min() < -1.337) + + +def test_clip_speckle(): + seed = 42 + data = camera() # 512x512 grayscale uint8 + data_signed = img_as_float(data) * 2. - 1. # Same image, on range [-1, 1] + + # Signed and unsigned, clipped + cam_speckle = random_noise(data, mode='speckle', seed=seed, clip=True) + cam_speckle2 = random_noise(data_signed, mode='speckle', seed=seed, + clip=True) + assert (cam_speckle.max() == 1.) and (cam_speckle.min() == 0.) + assert (cam_speckle2.max() == 1.) and (cam_speckle2.min() == -1.) + + # Signed and unsigned, unclipped + cam_speckle = random_noise(data, mode='speckle', seed=seed, clip=False) + cam_speckle2 = random_noise(data_signed, mode='speckle', seed=seed, + clip=False) + assert (cam_speckle.max() > 1.219) and (cam_speckle.min() == 0.) + assert (cam_speckle2.max() > 1.219) and (cam_speckle2.min() < -1.306) + + def test_bad_mode(): data = np.zeros((64, 64)) assert_raises(KeyError, random_noise, data, 'perlin') diff --git a/skimage/util/unique.py b/skimage/util/unique.py index 16a83f6b..635f6e89 100644 --- a/skimage/util/unique.py +++ b/skimage/util/unique.py @@ -30,9 +30,9 @@ def unique_rows(ar): Examples -------- >>> ar = np.array([[1, 0, 1], - [0, 1, 0], - [1, 0, 1]], np.uint8) - >>> aru = unique_rows(ar) + ... [0, 1, 0], + ... [1, 0, 1]], np.uint8) + >>> unique_rows(ar) array([[0, 1, 0], [1, 0, 1]], dtype=uint8) """ diff --git a/viewer_examples/plugins/watershed_demo.py b/viewer_examples/plugins/watershed_demo.py index 683e8a30..612ec6c9 100644 --- a/viewer_examples/plugins/watershed_demo.py +++ b/viewer_examples/plugins/watershed_demo.py @@ -7,6 +7,7 @@ from skimage.viewer import ImageViewer from skimage.viewer.widgets import history from skimage.viewer.plugins.labelplugin import LabelPainter + class OKCancelButtons(history.OKCancelButtons): def update_original_image(self):