wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-09 04:48:03 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge tag 'v0.9.3' into debian

Browse Source 
* tag 'v0.9.3': (164 commits)
  Set version to 0.9.3
  Merge pull request #796 from ahojnnes/warp-fix
  Set version to 0.9.2 for second try at PyPi upload.
  Set version to 0.9.1.
  Add missing files to MANIFEST for sdist upload.
  Update manifest not to include gh-pages in docs.
  Get rid of that inherited 's' once and for all.
  Update docversions correctly for 0.9.x.
  Mark BRIEF and Censure as experimental in release notes.
  Update gh-pages instructions in RELEASE.txt.
  Fix markup error in marching cubes docs.
  Correctly determine version number from module.
  Update version in docs.
  Update versions for 0.9.0 release.
  Update 0.9 release notes with new features.
  Contrib script now shows PRs and merges.
  Update contributors script to count by date.
  Speed up memory views in watershed function
  Speed up memory views in skeletonize function
  Speed up memory views in line drawing function
  ...
This commit is contained in:
Yaroslav Halchenko
2013-11-22 15:20:16 -05:00
parent 589681e023 af16654230
commit 5ef665a0b4
77 changed files with 2105 additions and 854 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CONTRIBUTORS.txt
+5 -1
View File
@@ -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
MANIFEST.in
+2 -1
View File
@@ -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
RELEASE.txt
+4 -4
View File
@@ -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 <commit>``,
where ``<commit>`` is the first commit since the previous release.
- To show a list of contributors and changes, run
``doc/release/contribs.py <tag of prev release>``.
- 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::
TODO.txt
+8 -8
View File
@@ -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`
bento.info
+1 -1
View File
@@ -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
doc/examples/plot_censure_keypoints.py
-54
View File
@@ -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()
doc/examples/plot_circular_elliptical_hough_transform.py
+21 -20
View File
@@ -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()
doc/examples/plot_join_segmentations.py
+6 -10
View File
@@ -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)
doc/examples/plot_marching_cubes.py
+1 -2
View File
@@ -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
doc/examples/plot_shapes.py
+58 -20
View File
@@ -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()
doc/gh-pages.py
+4
View File
@@ -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:
doc/release/contribs.py
Executable
+48
View File
@@ -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
doc/release/contributors.sh
-2
View File
@@ -1,2 +0,0 @@
git log $1..HEAD --format='- %aN' | sed 's/@/\-at\-/' | sed 's/<>//' | sort -u
doc/release/release_0.9.txt
+130
View File
@@ -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
doc/source/_static/docversions.js
+1 -1
View File
@@ -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++){
doc/source/api_changes.txt
+7
View File
@@ -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
-----------
doc/tools/build_modref_templates.py
+1 -1
View File
@@ -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'
doc/tools/plot_pr.py
+58 -38
View File
@@ -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()
setup.py
+1 -1
View File
@@ -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),
skimage/__init__.py
+1
View File
@@ -91,6 +91,7 @@ test_verbose.__doc__ = test.__doc__
class _Log(Warning):
pass
class _FakeLog(object):
def __init__(self, name):
"""
skimage/color/colorlabel.py
+54 -16
View File
@@ -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
skimage/color/tests/test_colorlabel.py
+22 -1
View File
@@ -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()
skimage/data/__init__.py
-1
View File
@@ -200,4 +200,3 @@ def coffee():
"""
return load("coffee.png")
skimage/draw/__init__.py
+6 -2
View File
@@ -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']
skimage/draw/_draw.pyx
+325 -43
View File
@@ -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 = <int>(sqrt(dx*dx + dy*dy))
x, y = x1, y1
while True:
cc.append(x)
rr.append(y)
val.append(1. * abs(err - dx + dy) / <float>(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) / <float>(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) / <float>(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, <Py_ssize_t>(dy), <Py_ssize_t>(dx),
<Py_ssize_t>(sy), <Py_ssize_t>(sx), cur)
return _bezier_segment(y0, x0, <Py_ssize_t>(dy), <Py_ssize_t>(dx),
<Py_ssize_t>(sy), <Py_ssize_t>(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 = <Py_ssize_t>(xx + x1)
y0 = y2
y2 = <Py_ssize_t>(yy + y1)
if (x0 == x2) or (weight == 1.):
t = <double>(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 = <int>(xx + 0.5)
y = <int>(yy + 0.5)
yy = (xx - x0) * (y1 - y0) / (x1 - x0) + y0
rr, cc = _bezier_segment(y0, x0, <int>(yy + 0.5), x, y, x, ww)
px.extend(rr)
py.extend(cc)
yy = (xx - x2) * (y1 - y2) / (x1 - x2) + y2
y1 = <int>(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 = <int>(xx + 0.5)
y = <int>(yy + 0.5)
xx = (x1 - x0) * (yy - y0) / (y1 - y0) + x0
rr, cc = _bezier_segment(y0, x0, y, <int>(xx + 0.5), y, x, ww)
px.extend(rr)
py.extend(cc)
xx = (x1 - x2) * (yy - y2) / (y1 - y2) + x2
x1 = <int>(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)
skimage/draw/draw.py
+1 -1
View File
@@ -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
skimage/draw/draw3d.py
+12 -14
View File
@@ -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
-------
skimage/draw/tests/test_draw.py
+185 -7
View File
@@ -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()
skimage/draw/tests/test_draw3d.py
+2 -2
View File
@@ -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]
skimage/exposure/exposure.py
+1 -1
View File
@@ -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
skimage/feature/__init__.py
+2 -8
View File
@@ -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']
skimage/feature/_brief.py
+6 -2
View File
@@ -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
skimage/feature/_texture.pyx
+11 -11
View File
@@ -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:
skimage/feature/censure.py
+2 -1
View File
@@ -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.
skimage/feature/tests/test_brief.py → skimage/feature/tests/_test_brief.py
View File
skimage/feature/tests/test_censure.py → skimage/feature/tests/_test_censure.py
View File
skimage/feature/tests/test_peak.py
+131
View File
@@ -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()
skimage/feature/util.py
+3 -1
View File
@@ -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
skimage/filter/__init__.py
+2 -3
View File
@@ -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',
skimage/filter/_denoise.py
+1 -6
View File
@@ -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)
skimage/filter/_rank_order.py
+6 -6
View File
@@ -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)
skimage/filter/ctmf.py
+43 -35
View File
@@ -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
skimage/filter/edges.py
+3 -3
View File
@@ -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
skimage/io/_plugins/gtk_plugin.py
+1 -1
View File
@@ -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:
skimage/io/_plugins/pil_plugin.py
+3 -1
View File
@@ -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())
skimage/io/tests/test_pil.py
+2 -2
View File
@@ -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)
skimage/measure/_marching_cubes.py
+12 -12
View File
@@ -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.
skimage/measure/_marching_cubes_cy.pyx
+36 -36
View File
@@ -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
skimage/measure/_regionprops.py
+28 -11
View File
@@ -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
skimage/measure/tests/test_marching_cubes.py
+5 -5
View File
@@ -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)
skimage/measure/tests/test_regionprops.py
+15 -3
View File
@@ -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()
skimage/morphology/__init__.py
+1 -2
View File
@@ -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',
skimage/morphology/_skeletonize.py
+4 -4
View File
@@ -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)
skimage/morphology/_skeletonize_cy.pyx
+2 -2
View File
@@ -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
skimage/morphology/_watershed.pyx
+3 -3
View File
@@ -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
skimage/morphology/binary.py
+44 -29
View File
@@ -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.
"""
skimage/morphology/convex_hull.py
+1 -1
View File
@@ -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)
skimage/morphology/tests/test_binary.py
+68
View File
@@ -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()
skimage/morphology/tests/test_grey.py
+1 -36
View File
@@ -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()
skimage/morphology/tests/test_watershed.py
+1 -98
View File
@@ -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()
skimage/morphology/watershed.py
+4 -75
View File
@@ -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
skimage/segmentation/__init__.py
+3 -2
View File
@@ -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']
skimage/segmentation/_join.py
+35 -15
View File
@@ -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
skimage/segmentation/random_walker_segmentation.py
+63 -25
View File
@@ -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
skimage/segmentation/tests/test_join.py
+28 -3
View File
@@ -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()
skimage/segmentation/tests/test_random_walker.py
+97 -23
View File
@@ -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()
skimage/transform/__init__.py
+1 -3
View File
@@ -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',
skimage/transform/_geometric.py
+17 -13
View File
@@ -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
skimage/transform/_hough_transform.pyx
+46 -37
View File
@@ -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,
skimage/transform/_warps_cy.pyx
+4 -4
View File
@@ -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)
skimage/transform/hough_transform.py
-24
View File
@@ -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.
skimage/transform/tests/test_hough_transform.py
+208 -33
View File
@@ -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__":
skimage/util/noise.py
+81 -14
View File
@@ -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
skimage/util/shape.py
+1 -3
View File
@@ -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)
skimage/util/tests/test_random_noise.py
+111 -3
View File
@@ -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')
skimage/util/unique.py
+3 -3
View File
@@ -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)
"""
viewer_examples/plugins/watershed_demo.py
+1
View File
@@ -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):