mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-02 15:23:16 +08:00
Merge tag 'v0.6' into debian
* tag 'v0.6': (154 commits) PKG: Update doc versions. PKG: Update version to 0.6. DOC: Add 0.6 release notes. Skip test that fails for PIL < 1.1.7 BUG: Fix structural similarity to use new signature for view_as_windows. Remove bad gradient check. BUG Fix invalid import of structural_similarity. BUG Remove merge artefact. BUG: Remove double import of find contours. PKG: Rename as_windows to view_as_windows. ENH: Promote as_windows to a utility function. STY: Wrap long line. PKG: Rename _ssim to _structural_similarity. ENH: Automatically determine dynamic range in ssim if not specified. DOC: In ssim example, use correct dynamic range in algorithm and when displaying results. DOC: Example of structural similarity. The example is currently broken because structural_similarity does not auto-detect dynamic range. ENH: Rename ssid to structural_similarity to avoid confusion with self-similarity features. ENH: Add SSIM gradient. ENH: Add Structural SIMilarity (SSIM) image comparison. Replace mpltools reference with skimage. ENH: Allow resetting the plugin state. ...
This commit is contained in:
Yaroslav Halchenko
+3
-1
@@ -28,7 +28,7 @@
|
||||
- Tony Yu
|
||||
Reading of paletted images; build, bug and doc fixes.
|
||||
Code to generate skimage logo.
|
||||
Otsu thresholding, histogram equalisation, and more.
|
||||
Otsu thresholding, histogram equalisation, template matching, and more.
|
||||
|
||||
- Zachary Pincus
|
||||
Tracing of low cost paths, FreeImage I/O plugin, iso-contours,
|
||||
@@ -46,6 +46,7 @@
|
||||
|
||||
- Pieter Holtzhausen
|
||||
Incorporating CellProfiler's Sobel edge detector, build and bug fixes.
|
||||
Radon transform, template matching.
|
||||
|
||||
- Emmanuelle Guillart
|
||||
Total variation noise filtering, integration of CellProfiler's
|
||||
@@ -105,3 +106,4 @@
|
||||
- Johannes Schönberger
|
||||
Polygon, circle and ellipse drawing functions
|
||||
Adaptive thresholding
|
||||
Implementation of Matlab's `regionprops`
|
||||
|
||||
+1
-1
@@ -2,7 +2,7 @@ How to make a new release of ``skimage``
|
||||
========================================
|
||||
|
||||
- Update release notes.
|
||||
- Update the version number in setup.py and commit
|
||||
- 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
|
||||
- Build a clean version of the docs. Run "make" in the root dir, then
|
||||
|
||||
+100
@@ -0,0 +1,100 @@
|
||||
Name: scikits-image
|
||||
Version: 0.6
|
||||
Summary: Image processing routines for SciPy
|
||||
Url: http://scikits-image.org
|
||||
DownloadUrl: http://github.com/scikits-image/scikits-image
|
||||
Description: Image Processing SciKit
|
||||
|
||||
Image processing algorithms for SciPy, including IO, morphology, filtering,
|
||||
warping, color manipulation, object detection, etc.
|
||||
|
||||
Please refer to the online documentation at
|
||||
http://scikits-image.org/
|
||||
Maintainer: Stefan van der Walt
|
||||
MaintainerEmail: stefan@sun.ac.za
|
||||
License: Modified BSD
|
||||
Classifiers:
|
||||
Development Status :: 4 - Beta,
|
||||
Environment :: Console,
|
||||
Intended Audience :: Developers,
|
||||
Intended Audience :: Science/Research,
|
||||
License :: OSI Approved :: BSD License,
|
||||
Programming Language :: C,
|
||||
Programming Language :: Python,
|
||||
Programming Language :: Python :: 3,
|
||||
Topic :: Scientific/Engineering,
|
||||
Operating System :: Microsoft :: Windows,
|
||||
Operating System :: POSIX,
|
||||
Operating System :: Unix,
|
||||
Operating System :: MacOS
|
||||
|
||||
HookFile: bscript
|
||||
UseBackends: Waf
|
||||
|
||||
Library:
|
||||
Packages:
|
||||
skimage, skimage.color, skimage.data, skimage.draw, skimage.exposure,
|
||||
skimage.feature, skimage.filter, skimage.graph, skimage.io,
|
||||
skimage.io._plugins, skimage.measure, skimage.morphology,
|
||||
skimage.scripts, skimage.segmentation, skimage.transform, skimage.util
|
||||
Extension: skimage.morphology._pnpoly
|
||||
Sources:
|
||||
skimage/morphology/_pnpoly.pyx
|
||||
Extension: skimage.feature._greycomatrix
|
||||
Sources:
|
||||
skimage/feature/_greycomatrix.pyx
|
||||
Extension: skimage.feature._template
|
||||
Sources:
|
||||
skimage/feature/_template.pyx
|
||||
Extension: skimage.io._plugins._colormixer
|
||||
Sources:
|
||||
skimage/io/_plugins/_colormixer.pyx
|
||||
Extension: skimage.measure._find_contours
|
||||
Sources:
|
||||
skimage/measure/_find_contours.pyx
|
||||
Extension: skimage.measure._moments
|
||||
Sources:
|
||||
skimage/measure/_moments.pyx
|
||||
Extension: skimage.graph._mcp
|
||||
Sources:
|
||||
skimage/graph/_mcp.pyx
|
||||
Extension: skimage.io._plugins._histograms
|
||||
Sources:
|
||||
skimage/io/_plugins/_histograms.pyx
|
||||
Extension: skimage.transform._hough_transform
|
||||
Sources:
|
||||
skimage/transform/_hough_transform.pyx
|
||||
Extension: skimage.filter._ctmf
|
||||
Sources:
|
||||
skimage/filter/_ctmf.pyx
|
||||
Extension: skimage.morphology.ccomp
|
||||
Sources:
|
||||
skimage/morphology/ccomp.pyx
|
||||
Extension: skimage.morphology._watershed
|
||||
Sources:
|
||||
skimage/morphology/_watershed.pyx
|
||||
Extension: skimage.morphology._convex_hull
|
||||
Sources:
|
||||
skimage/morphology/_convex_hull.pyx
|
||||
Extension: skimage.morphology._skeletonize
|
||||
Sources:
|
||||
skimage/morphology/_skeletonize.pyx
|
||||
Extension: skimage.draw._draw
|
||||
Sources:
|
||||
skimage/draw/_draw.pyx
|
||||
Extension: skimage.transform._project
|
||||
Sources:
|
||||
skimage/transform/_project.pyx
|
||||
Extension: skimage.graph._spath
|
||||
Sources:
|
||||
skimage/graph/_spath.pyx
|
||||
Extension: skimage.morphology.cmorph
|
||||
Sources:
|
||||
skimage/morphology/cmorph.pyx
|
||||
Extension: skimage.graph.heap
|
||||
Sources:
|
||||
skimage/graph/heap.pyx
|
||||
|
||||
Executable: skivi
|
||||
Module: skimage.scripts.skivi
|
||||
Function: main
|
||||
@@ -0,0 +1,12 @@
|
||||
import os.path as op
|
||||
|
||||
from numpy.distutils.misc_util \
|
||||
import \
|
||||
get_numpy_include_dirs
|
||||
|
||||
from bento.commands import hooks
|
||||
|
||||
@hooks.post_configure
|
||||
def post_configure(context):
|
||||
conf = context.waf_context
|
||||
conf.env.INCLUDES = get_numpy_include_dirs() + [op.join("skimage", "morphology")]
|
||||
@@ -3,57 +3,40 @@
|
||||
Comparing edge-based segmentation and region-based segmentation
|
||||
===============================================================
|
||||
|
||||
In this example, we will see how to segment objects from a background.
|
||||
We use the ``coins`` image from ``skimage.data``. This image shows
|
||||
several coins outlined against a darker background. The segmentation of
|
||||
the coins cannot be done directly from the histogram of grey values,
|
||||
because the background shares enough grey levels with the coins that a
|
||||
thresholding segmentation is not sufficient. Simply thresholding the image
|
||||
leads either to missing significant parts of the coins, or to merging parts
|
||||
of the background with the coins.
|
||||
|
||||
We first try an edge-based segmentation. We use the Canny detector to
|
||||
delineate the contours of the coins. These contours are filled using
|
||||
mathematical morphology (``scipy.ndimage.binary_fill_holes``). Small spurious
|
||||
objects are easily removed by applying a threshold on the size of
|
||||
unconnected objects. However, this method is not very robust, since contours
|
||||
that are not perfectly closed are not filled correctly. This happens for one
|
||||
of the coins.
|
||||
|
||||
We therefore try a second method, that is region-based. Here we use the
|
||||
watershed transform. An elevation map is provided by the Sobel gradient
|
||||
of the image. Markers of the background and the coins are determined from
|
||||
the extreme parts of the histogram of grey values.
|
||||
|
||||
This second method works even better, and the coins can be segmented and
|
||||
labeled individually.
|
||||
|
||||
In this example, we will see how to segment objects from a background. We use
|
||||
the ``coins`` image from ``skimage.data``, which shows several coins outlined
|
||||
against a darker background.
|
||||
"""
|
||||
|
||||
|
||||
import numpy as np
|
||||
from scipy import ndimage
|
||||
import matplotlib.pyplot as plt
|
||||
import skimage
|
||||
from skimage.filter import canny, sobel
|
||||
from skimage.morphology import watershed
|
||||
|
||||
#------------------ Loading data --------------------------------
|
||||
from skimage import data
|
||||
coins = data.coins()
|
||||
|
||||
#------------ Histogram of grey values ---------------------------
|
||||
histo = np.histogram(coins, bins=np.arange(0, 256))
|
||||
coins = data.coins()
|
||||
hist = np.histogram(coins, bins=np.arange(0, 256))
|
||||
|
||||
plt.figure(figsize=(8, 3))
|
||||
plt.subplot(121)
|
||||
plt.imshow(coins, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
plt.subplot(122)
|
||||
plt.plot(histo[1][:-1], histo[0], lw=2)
|
||||
plt.plot(hist[1][:-1], hist[0], lw=2)
|
||||
plt.title('histogram of grey values')
|
||||
|
||||
#------------------ Tentative thresholding --------------------------------
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Thresholding
|
||||
============
|
||||
|
||||
A simple way to segment the coins is to choose a threshold based on the
|
||||
histogram of grey values. Unfortunately, thresholding this image gives a binary
|
||||
image that either misses significant parts of the coins or merges parts of the
|
||||
background with the coins:
|
||||
|
||||
"""
|
||||
|
||||
plt.figure(figsize=(6, 3))
|
||||
plt.subplot(121)
|
||||
plt.imshow(coins > 100, cmap=plt.cm.gray, interpolation='nearest')
|
||||
@@ -63,73 +46,126 @@ plt.subplot(122)
|
||||
plt.imshow(coins > 150, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.title('coins > 150')
|
||||
plt.axis('off')
|
||||
margins = dict(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
|
||||
plt.subplots_adjust(**margins)
|
||||
|
||||
plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
|
||||
right=1)
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
|
||||
#------------------ Edge-based segmentation --------------------------------
|
||||
Edge-based segmentation
|
||||
=======================
|
||||
|
||||
Next, we try to delineate the contours of the coins using edge-based
|
||||
segmentation. To do this, we first get the edges of features using the Canny
|
||||
edge-detector.
|
||||
"""
|
||||
|
||||
from skimage.filter import canny
|
||||
edges = canny(coins/255.)
|
||||
|
||||
plt.figure(figsize=(4, 3))
|
||||
plt.imshow(edges, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
plt.title('Canny detector')
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
These contours are then filled using mathematical morphology.
|
||||
"""
|
||||
|
||||
from scipy import ndimage
|
||||
|
||||
fill_coins = ndimage.binary_fill_holes(edges)
|
||||
|
||||
plt.figure(figsize=(4, 3))
|
||||
plt.imshow(fill_coins, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
plt.title('Filling the holes')
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Small spurious objects are easily removed by setting a minimum size for valid
|
||||
objects.
|
||||
"""
|
||||
|
||||
label_objects, nb_labels = ndimage.label(fill_coins)
|
||||
sizes = np.bincount(label_objects.ravel())
|
||||
mask_sizes = sizes > 20
|
||||
mask_sizes[0] = 0
|
||||
coins_cleaned = mask_sizes[label_objects]
|
||||
|
||||
plt.figure(figsize=(7, 3.))
|
||||
plt.subplot(131)
|
||||
plt.imshow(edges, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
plt.title('Canny detector')
|
||||
plt.subplot(132)
|
||||
plt.imshow(fill_coins, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
plt.title('Filling the holes')
|
||||
plt.subplot(133)
|
||||
plt.figure(figsize=(4, 3))
|
||||
plt.imshow(coins_cleaned, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
plt.title('Removing small objects')
|
||||
plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
|
||||
right=1)
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
However, this method is not very robust, since contours that are not perfectly
|
||||
closed are not filled correctly, as is the case for one unfilled coin above.
|
||||
|
||||
|
||||
#------------------ Region-based segmentation --------------------------------
|
||||
Region-based segmentation
|
||||
=========================
|
||||
|
||||
We therefore try a region-based method using the
|
||||
watershed transform. First, we find an elevation map using the Sobel gradient of the
|
||||
image.
|
||||
|
||||
"""
|
||||
|
||||
from skimage.filter import sobel
|
||||
|
||||
elevation_map = sobel(coins)
|
||||
|
||||
plt.figure(figsize=(4, 3))
|
||||
plt.imshow(elevation_map, cmap=plt.cm.jet, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
plt.title('elevation_map')
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Next we find markers of the background and the coins based on the extreme parts
|
||||
of the histogram of grey values.
|
||||
"""
|
||||
|
||||
markers = np.zeros_like(coins)
|
||||
markers[coins < 30] = 1
|
||||
markers[coins > 150] = 2
|
||||
|
||||
elevation_map = sobel(coins)
|
||||
|
||||
|
||||
segmentation = watershed(elevation_map, markers)
|
||||
|
||||
|
||||
plt.figure(figsize=(7, 3))
|
||||
plt.subplot(131)
|
||||
plt.figure(figsize=(4, 3))
|
||||
plt.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
plt.title('markers')
|
||||
plt.subplot(132)
|
||||
plt.imshow(elevation_map, cmap=plt.cm.jet, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
plt.title('elevation_map')
|
||||
plt.subplot(133)
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Finally, we use the watershed transform to fill regions of the elevation map starting from the markers determined above:
|
||||
|
||||
"""
|
||||
from skimage.morphology import watershed
|
||||
segmentation = watershed(elevation_map, markers)
|
||||
|
||||
plt.figure(figsize=(4, 3))
|
||||
plt.imshow(segmentation, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
plt.title('segmentation')
|
||||
|
||||
plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
|
||||
right=1)
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
# ------------------- Removing a few small holes ---------------------
|
||||
This last method works even better, and the coins can be segmented and
|
||||
labeled individually.
|
||||
|
||||
"""
|
||||
|
||||
segmentation = ndimage.binary_fill_holes(segmentation - 1)
|
||||
|
||||
#------------------ Labeling the coins --------------------------------
|
||||
labeled_coins, _ = ndimage.label(segmentation)
|
||||
|
||||
plt.figure(figsize=(6, 3))
|
||||
@@ -141,8 +177,11 @@ plt.subplot(122)
|
||||
plt.imshow(labeled_coins, cmap=plt.cm.spectral, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
|
||||
right=1)
|
||||
plt.subplots_adjust(**margins)
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
"""
|
||||
|
||||
plt.show()
|
||||
|
||||
@@ -0,0 +1,64 @@
|
||||
"""
|
||||
=========================
|
||||
Measure region properties
|
||||
=========================
|
||||
|
||||
This example shows how to measure properties of labelled image regions.
|
||||
"""
|
||||
|
||||
import math
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
|
||||
from skimage.draw import ellipse
|
||||
from skimage.morphology import label
|
||||
from skimage.measure import regionprops
|
||||
from scipy.ndimage import geometric_transform
|
||||
|
||||
|
||||
ANGLE = 0.2
|
||||
|
||||
def rotate(xy):
|
||||
x, y = xy
|
||||
out_x = math.cos(ANGLE) * x - math.sin(ANGLE) * y
|
||||
out_y = math.sin(ANGLE) * x + math.cos(ANGLE) * y
|
||||
return (out_x, out_y)
|
||||
|
||||
image = np.zeros((600, 600), 'int')
|
||||
|
||||
rr, cc = ellipse(300, 350, 100, 220)
|
||||
image[rr,cc] = 1
|
||||
|
||||
image = geometric_transform(image, rotate)
|
||||
|
||||
label_img = label(image)
|
||||
props = regionprops(label_img, [
|
||||
'BoundingBox',
|
||||
'Centroid',
|
||||
'Orientation',
|
||||
'MajorAxisLength',
|
||||
'MinorAxisLength'
|
||||
])
|
||||
|
||||
plt.imshow(image)
|
||||
|
||||
for prop in props:
|
||||
x0 = prop['Centroid'][1]
|
||||
y0 = prop['Centroid'][0]
|
||||
x1 = x0 + math.cos(prop['Orientation']) * 0.5 * prop['MajorAxisLength']
|
||||
y1 = y0 - math.sin(prop['Orientation']) * 0.5 * prop['MajorAxisLength']
|
||||
x2 = x0 - math.sin(prop['Orientation']) * 0.5 * prop['MinorAxisLength']
|
||||
y2 = y0 - math.cos(prop['Orientation']) * 0.5 * prop['MinorAxisLength']
|
||||
|
||||
plt.plot((x0, x1), (y0, y1), '-r', linewidth=2.5)
|
||||
plt.plot((x0, x2), (y0, y2), '-r', linewidth=2.5)
|
||||
plt.plot(x0, y0, '.g', markersize=15)
|
||||
|
||||
minr, minc, maxr, maxc = prop['BoundingBox']
|
||||
bx = (minc, maxc, maxc, minc, minc)
|
||||
by = (minr, minr, maxr, maxr, minr)
|
||||
plt.plot(bx, by, '-b', linewidth=2.5)
|
||||
|
||||
plt.gray()
|
||||
plt.axis((0, 600, 600, 0))
|
||||
plt.show()
|
||||
@@ -0,0 +1,68 @@
|
||||
'''
|
||||
===========================
|
||||
Structural similarity index
|
||||
===========================
|
||||
|
||||
When comparing images, the mean squared error (MSE)--while simple to
|
||||
implement--is not highly indicative of perceived similarity. Structural
|
||||
similarity aims to address this shortcoming by taking texture into account
|
||||
[1]_, [2]_.
|
||||
|
||||
The example shows two modifications of the input image, each with the same MSE,
|
||||
but with very different mean structural similarity indices.
|
||||
|
||||
.. [1] Zhou Wang; Bovik, A.C.; ,"Mean squared error: Love it or leave it? A new
|
||||
look at Signal Fidelity Measures," Signal Processing Magazine, IEEE,
|
||||
vol. 26, no. 1, pp. 98-117, Jan. 2009.
|
||||
|
||||
.. [2] Z. Wang, A. C. Bovik, H. R. Sheikh and E. P. Simoncelli, "Image quality
|
||||
assessment: From error visibility to structural similarity," IEEE
|
||||
Transactions on Image Processing, vol. 13, no. 4, pp. 600-612,
|
||||
Apr. 2004.
|
||||
|
||||
'''
|
||||
|
||||
from skimage import data, color, io, exposure, img_as_float
|
||||
from skimage.measure import structural_similarity as ssim
|
||||
|
||||
import numpy as np
|
||||
|
||||
img = img_as_float(data.camera())
|
||||
rows, cols = img.shape
|
||||
|
||||
noise = np.ones_like(img) * 0.2 * (img.max() - img.min())
|
||||
noise[np.random.random(size=noise.shape) > 0.5] *= -1
|
||||
|
||||
def mse(x, y):
|
||||
return np.linalg.norm(x - y)
|
||||
|
||||
img_noise = img + noise
|
||||
img_const = img + abs(noise)
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
f, (ax0, ax1, ax2) = plt.subplots(1, 3)
|
||||
|
||||
mse_none = mse(img, img)
|
||||
ssim_none = ssim(img, img, dynamic_range=img.max() - img.min())
|
||||
|
||||
mse_noise = mse(img, img_noise)
|
||||
ssim_noise = ssim(img, img_noise, dynamic_range=img_const.max() - img_const.min())
|
||||
|
||||
mse_const = mse(img, img_const)
|
||||
ssim_const = ssim(img, img_const, dynamic_range=img_noise.max() - img_noise.min())
|
||||
|
||||
label = 'MSE: %2.f, SSIM: %.2f'
|
||||
|
||||
ax0.imshow(img, cmap=plt.cm.gray, vmin=0, vmax=1)
|
||||
ax0.set_xlabel(label % (mse_none, ssim_none))
|
||||
ax0.set_title('Original image')
|
||||
|
||||
ax1.imshow(img_noise, cmap=plt.cm.gray, vmin=0, vmax=1)
|
||||
ax1.set_xlabel(label % (mse_noise, ssim_noise))
|
||||
ax1.set_title('Image with noise')
|
||||
|
||||
ax2.imshow(img_const, cmap=plt.cm.gray, vmin=0, vmax=1)
|
||||
ax2.set_xlabel(label % (mse_const, ssim_const))
|
||||
ax2.set_title('Image plus constant')
|
||||
|
||||
plt.show()
|
||||
@@ -0,0 +1,83 @@
|
||||
r"""
|
||||
=====
|
||||
Swirl
|
||||
=====
|
||||
|
||||
Image swirling is a non-linear image deformation that creates a whirlpool
|
||||
effect. This example describes the implementation of this transform in
|
||||
``skimage``, as well as the underlying warp mechanism.
|
||||
|
||||
Image warping
|
||||
`````````````
|
||||
When applying a geometric transformation on an image, we typically make use of
|
||||
a reverse mapping, i.e., for each pixel in the output image, we compute its
|
||||
corresponding position in the input. The reason is that, if we were to do it
|
||||
the other way around (map each input pixel to its new output position), some
|
||||
pixels in the output may be left empty. On the other hand, each output
|
||||
coordinate has exactly one corresponding location in (or outside) the input
|
||||
image, and even if that position is non-integer, we may use interpolation to
|
||||
compute the corresponding image value.
|
||||
|
||||
Performing a reverse mapping
|
||||
````````````````````````````
|
||||
To perform a geometric warp in ``skimage``, you simply need to provide the
|
||||
reverse mapping to the ``skimage.transform.warp`` function. E.g., consider the
|
||||
case where we would like to shift an image 50 pixels to the left. The reverse
|
||||
mapping for such a shift would be::
|
||||
|
||||
def shift_left(xy):
|
||||
xy[:, 0] += 50
|
||||
return xy
|
||||
|
||||
The corresponding call to warp is::
|
||||
|
||||
from skimage.transform import warp
|
||||
warp(image, shift_left)
|
||||
|
||||
The swirl transformation
|
||||
````````````````````````
|
||||
Consider the coordinate :math:`(x, y)` in the output image. The reverse
|
||||
mapping for the swirl transformation first computes, relative to a center
|
||||
:math:`(x_0, y_0)`, its polar coordinates,
|
||||
|
||||
.. math::
|
||||
|
||||
\theta = \arctan(y/x)
|
||||
|
||||
\rho = \sqrt{(x - x_0)^2 + (y - y_0)^2},
|
||||
|
||||
and then transforms them according to
|
||||
|
||||
.. math::
|
||||
|
||||
r = \ln(2) \, \mathtt{radius} / 5
|
||||
|
||||
\phi = \mathtt{rotation}
|
||||
|
||||
s = \mathtt{strength}
|
||||
|
||||
\theta' = \phi + s \, e^{-\rho / r + \theta}
|
||||
|
||||
where ``strength`` is a parameter for the amount of swirl, ``radius`` indicates
|
||||
the swirl extent in pixels, and ``rotation`` adds a rotation angle. The
|
||||
transformation of ``radius`` into :math:`r` is to ensure that the
|
||||
transformation decays to :math:`\approx 1/1000^{\mathsf{th}}` within the
|
||||
specified radius.
|
||||
"""
|
||||
|
||||
from skimage import data
|
||||
from skimage.transform import swirl
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
image = data.checkerboard()
|
||||
swirled = swirl(image, rotation=0, strength=10, radius=120, order=2)
|
||||
|
||||
f, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 3))
|
||||
|
||||
ax0.imshow(image, cmap=plt.cm.gray, interpolation='none')
|
||||
ax0.axis('off')
|
||||
ax1.imshow(swirled, cmap=plt.cm.gray, interpolation='none')
|
||||
ax1.axis('off')
|
||||
|
||||
plt.show()
|
||||
@@ -0,0 +1,56 @@
|
||||
"""
|
||||
=================
|
||||
Template Matching
|
||||
=================
|
||||
|
||||
In this example, we use template matching to identify the occurrence of an
|
||||
image patch (in this case, a sub-image centered on a single coin). Here, we
|
||||
return a single match (the exact same coin), so the maximum value in the
|
||||
``match_template`` result corresponds to the coin location. The other coins
|
||||
look similar, and thus have local maxima; if you expect multiple matches, you
|
||||
should use a proper peak-finding function.
|
||||
|
||||
The ``match_template`` function uses fast, normalized cross-correlation [1]_
|
||||
to find instances of the template in the image. Note that the peaks in the
|
||||
output of ``match_template`` correspond to the origin (i.e. top-left corner) of
|
||||
the template.
|
||||
|
||||
.. [1] J. P. Lewis, "Fast Normalized Cross-Correlation", Industrial Light and
|
||||
Magic.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from skimage import data
|
||||
from skimage.feature import match_template
|
||||
|
||||
image = data.coins()
|
||||
coin = image[170:220, 75:130]
|
||||
|
||||
result = match_template(image, coin)
|
||||
ij = np.unravel_index(np.argmax(result), result.shape)
|
||||
x, y = ij[::-1]
|
||||
|
||||
fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 3))
|
||||
|
||||
ax1.imshow(coin)
|
||||
ax1.set_axis_off()
|
||||
ax1.set_title('template')
|
||||
|
||||
ax2.imshow(image)
|
||||
ax2.set_axis_off()
|
||||
ax2.set_title('image')
|
||||
# highlight matched region
|
||||
hcoin, wcoin = coin.shape
|
||||
rect = plt.Rectangle((x, y), wcoin, hcoin, edgecolor='r', facecolor='none')
|
||||
ax2.add_patch(rect)
|
||||
|
||||
ax3.imshow(result)
|
||||
ax3.set_axis_off()
|
||||
ax3.set_title('`match_template`\nresult')
|
||||
# highlight matched region
|
||||
ax3.autoscale(False)
|
||||
ax3.plot(x, y, 'o', markeredgecolor='r', markerfacecolor='none', markersize=10)
|
||||
|
||||
plt.show()
|
||||
|
||||
@@ -21,7 +21,7 @@ line.
|
||||
|
||||
See Wikipedia_ for more details on the algorithm.
|
||||
|
||||
.. _Wikipedia: <http://en.wikipedia.org/wiki/Watershed_(image_processing)>
|
||||
.. _Wikipedia: http://en.wikipedia.org/wiki/Watershed_(image_processing)
|
||||
|
||||
"""
|
||||
|
||||
|
||||
@@ -1,299 +0,0 @@
|
||||
"""
|
||||
Example generation for the scikit image.
|
||||
|
||||
Generate the rst files for the examples by iterating over the python
|
||||
example files.
|
||||
|
||||
Files that generate images should start with 'plot'.
|
||||
|
||||
This code was taken from the scikit-learn.
|
||||
|
||||
"""
|
||||
import os
|
||||
import shutil
|
||||
import traceback
|
||||
import glob
|
||||
|
||||
import matplotlib
|
||||
matplotlib.use('Agg')
|
||||
|
||||
import token, tokenize
|
||||
|
||||
rst_template = """
|
||||
|
||||
.. _example_%(short_fname)s:
|
||||
|
||||
%(docstring)s
|
||||
|
||||
**Python source code:** :download:`%(fname)s <%(fname)s>`
|
||||
(generated using ``skimage`` |version|)
|
||||
|
||||
.. literalinclude:: %(fname)s
|
||||
:lines: %(end_row)s-
|
||||
"""
|
||||
|
||||
plot_rst_template = """
|
||||
|
||||
.. _example_%(short_fname)s:
|
||||
|
||||
%(docstring)s
|
||||
|
||||
%(image_list)s
|
||||
|
||||
**Python source code:** :download:`%(fname)s <%(fname)s>`
|
||||
(generated using ``skimage`` |version|)
|
||||
|
||||
.. literalinclude:: %(fname)s
|
||||
:lines: %(end_row)s-
|
||||
"""
|
||||
|
||||
|
||||
toctree_template = """
|
||||
|
||||
.. toctree::
|
||||
:hidden:
|
||||
|
||||
%s
|
||||
|
||||
"""
|
||||
|
||||
# The following strings are used when we have several pictures: we use
|
||||
# an html div tag that our CSS uses to turn the lists into horizontal
|
||||
# lists.
|
||||
HLIST_HEADER = """
|
||||
.. rst-class:: horizontal
|
||||
|
||||
"""
|
||||
|
||||
HLIST_IMAGE_TEMPLATE = """
|
||||
*
|
||||
|
||||
.. image:: images/%s
|
||||
:scale: 50
|
||||
"""
|
||||
|
||||
SINGLE_IMAGE = """
|
||||
.. image:: images/%s
|
||||
:align: center
|
||||
"""
|
||||
|
||||
def extract_docstring(filename):
|
||||
""" Extract a module-level docstring, if any
|
||||
"""
|
||||
lines = file(filename).readlines()
|
||||
start_row = 0
|
||||
if lines[0].startswith('#!'):
|
||||
lines.pop(0)
|
||||
start_row = 1
|
||||
|
||||
docstring = ''
|
||||
first_par = ''
|
||||
tokens = tokenize.generate_tokens(lines.__iter__().next)
|
||||
for tok_type, tok_content, _, (erow, _), _ in tokens:
|
||||
tok_type = token.tok_name[tok_type]
|
||||
if tok_type in ('NEWLINE', 'COMMENT', 'NL', 'INDENT', 'DEDENT'):
|
||||
continue
|
||||
elif tok_type == 'STRING':
|
||||
docstring = eval(tok_content)
|
||||
# If the docstring is formatted with several paragraphs, extract
|
||||
# the first one:
|
||||
paragraphs = '\n'.join(line.rstrip()
|
||||
for line in docstring.split('\n')).split('\n\n')
|
||||
if len(paragraphs) > 0:
|
||||
first_par = paragraphs[0]
|
||||
break
|
||||
return docstring, first_par, erow+1+start_row
|
||||
|
||||
|
||||
def generate_example_rst(app):
|
||||
""" Generate the list of examples, as well as the contents of
|
||||
examples.
|
||||
"""
|
||||
root_dir = os.path.join(app.builder.srcdir, 'auto_examples')
|
||||
example_dir = os.path.abspath(app.builder.srcdir + '/../' + 'examples')
|
||||
try:
|
||||
plot_gallery = eval(app.builder.config.plot_gallery)
|
||||
except TypeError:
|
||||
plot_gallery = bool(app.builder.config.plot_gallery)
|
||||
if not os.path.exists(example_dir):
|
||||
os.makedirs(example_dir)
|
||||
if not os.path.exists(root_dir):
|
||||
os.makedirs(root_dir)
|
||||
|
||||
# we create an index.rst with all examples
|
||||
fhindex = file(os.path.join(root_dir, 'index.txt'), 'w')
|
||||
fhindex.write("""\
|
||||
|
||||
Examples
|
||||
========
|
||||
|
||||
.. _examples-index:
|
||||
""")
|
||||
# Here we don't use an os.walk, but we recurse only twice: flat is
|
||||
# better than nested.
|
||||
generate_dir_rst('.', fhindex, example_dir, root_dir, plot_gallery)
|
||||
for dir in sorted(os.listdir(example_dir)):
|
||||
if os.path.isdir(os.path.join(example_dir, dir)):
|
||||
generate_dir_rst(dir, fhindex, example_dir, root_dir, plot_gallery)
|
||||
fhindex.flush()
|
||||
|
||||
|
||||
def generate_dir_rst(dir, fhindex, example_dir, root_dir, plot_gallery):
|
||||
""" Generate the rst file for an example directory.
|
||||
"""
|
||||
if not dir == '.':
|
||||
target_dir = os.path.join(root_dir, dir)
|
||||
src_dir = os.path.join(example_dir, dir)
|
||||
else:
|
||||
target_dir = root_dir
|
||||
src_dir = example_dir
|
||||
if not os.path.exists(os.path.join(src_dir, 'README.txt')):
|
||||
print 80*'_'
|
||||
print ('Example directory %s does not have a README.txt file'
|
||||
% src_dir)
|
||||
print 'Skipping this directory'
|
||||
print 80*'_'
|
||||
return
|
||||
|
||||
example_description = file(os.path.join(src_dir, 'README.txt')).read()
|
||||
fhindex.write("""\n\n%s\n\n""" % example_description)
|
||||
|
||||
if not os.path.exists(target_dir):
|
||||
os.makedirs(target_dir)
|
||||
|
||||
def sort_key(a):
|
||||
# put last elements without a plot
|
||||
if not a.startswith('plot') and a.endswith('.py'):
|
||||
return 'zz' + a
|
||||
return a
|
||||
|
||||
examples = [fname for fname in sorted(os.listdir(src_dir), key=sort_key)
|
||||
if fname.endswith('py')]
|
||||
ex_names = [ex[:-3] for ex in examples] # strip '.py' extension
|
||||
fhindex.write(toctree_template % '\n '.join(ex_names))
|
||||
|
||||
for fname in examples:
|
||||
generate_file_rst(fname, target_dir, src_dir, plot_gallery)
|
||||
thumb = os.path.join(dir, 'images', 'thumb', fname[:-3] + '.png')
|
||||
link_name = os.path.join(dir, fname).replace(os.path.sep, '_')
|
||||
fhindex.write('.. figure:: %s\n' % thumb)
|
||||
if link_name.startswith('._'):
|
||||
link_name = link_name[2:]
|
||||
fhindex.write(' :figclass: gallery\n')
|
||||
if dir != '.':
|
||||
fhindex.write(' :target: ./%s/%s.html\n\n' % (dir, fname[:-3]))
|
||||
else:
|
||||
fhindex.write(' :target: ./%s.html\n\n' % link_name[:-3])
|
||||
fhindex.write(' :ref:`example_%s`\n\n' % link_name)
|
||||
fhindex.write("""
|
||||
.. raw:: html
|
||||
|
||||
<div style="clear: both"></div>
|
||||
""") # clear at the end of the section
|
||||
|
||||
|
||||
def generate_file_rst(fname, target_dir, src_dir, plot_gallery):
|
||||
""" Generate the rst file for a given example.
|
||||
"""
|
||||
base_image_name = os.path.splitext(fname)[0]
|
||||
image_fname = '%s_%%s.png' % base_image_name
|
||||
|
||||
this_template = rst_template
|
||||
last_dir = os.path.split(src_dir)[-1]
|
||||
# to avoid leading . in file names, and wrong names in links
|
||||
if last_dir == '.' or last_dir == 'examples':
|
||||
last_dir = ''
|
||||
else:
|
||||
last_dir += '_'
|
||||
short_fname = last_dir + fname
|
||||
src_file = os.path.join(src_dir, fname)
|
||||
example_file = os.path.join(target_dir, fname)
|
||||
shutil.copyfile(src_file, example_file)
|
||||
|
||||
# The following is a list containing all the figure names
|
||||
figure_list = []
|
||||
|
||||
image_dir = os.path.join(target_dir, 'images')
|
||||
thumb_dir = os.path.join(image_dir, 'thumb')
|
||||
if not os.path.exists(image_dir):
|
||||
os.makedirs(image_dir)
|
||||
if not os.path.exists(thumb_dir):
|
||||
os.makedirs(thumb_dir)
|
||||
image_path = os.path.join(image_dir, image_fname)
|
||||
thumb_file = os.path.join(thumb_dir, fname[:-3] + '.png')
|
||||
if plot_gallery and fname.startswith('plot'):
|
||||
# generate the plot as png image if file name
|
||||
# starts with plot and if it is more recent than an
|
||||
# existing image.
|
||||
first_image_file = image_path % 1
|
||||
|
||||
if (not os.path.exists(first_image_file) or
|
||||
os.stat(first_image_file).st_mtime <=
|
||||
os.stat(src_file).st_mtime):
|
||||
# We need to execute the code
|
||||
print 'plotting %s' % fname
|
||||
import matplotlib.pyplot as plt
|
||||
plt.close('all')
|
||||
cwd = os.getcwd()
|
||||
try:
|
||||
# First CD in the original example dir, so that any file created
|
||||
# by the example get created in this directory
|
||||
os.chdir(os.path.dirname(src_file))
|
||||
execfile(os.path.basename(src_file), {'pl' : plt})
|
||||
os.chdir(cwd)
|
||||
|
||||
# In order to save every figure we have two solutions :
|
||||
# * iterate from 1 to infinity and call plt.fignum_exists(n)
|
||||
# (this requires the figures to be numbered
|
||||
# incrementally: 1, 2, 3 and not 1, 2, 5)
|
||||
# * iterate over [fig_mngr.num for fig_mngr in
|
||||
# matplotlib._pylab_helpers.Gcf.get_all_fig_managers()]
|
||||
for fig_num in (fig_mngr.num for fig_mngr in
|
||||
matplotlib._pylab_helpers.Gcf.get_all_fig_managers()):
|
||||
# Set the fig_num figure as the current figure as we can't
|
||||
# save a figure that's not the current figure.
|
||||
plt.figure(fig_num)
|
||||
plt.savefig(image_path % fig_num)
|
||||
figure_list.append(image_fname % fig_num)
|
||||
except:
|
||||
print 80*'_'
|
||||
print '%s is not compiling:' % fname
|
||||
traceback.print_exc()
|
||||
print 80*'_'
|
||||
finally:
|
||||
os.chdir(cwd)
|
||||
else:
|
||||
figure_list = [f[len(image_dir):]
|
||||
for f in glob.glob(image_path % '[1-9]')]
|
||||
#for f in glob.glob(image_path % '*')]
|
||||
|
||||
# generate thumb file
|
||||
this_template = plot_rst_template
|
||||
from matplotlib import image
|
||||
if os.path.exists(first_image_file):
|
||||
image.thumbnail(first_image_file, thumb_file, 0.25)
|
||||
|
||||
if not os.path.exists(thumb_file):
|
||||
# create something not to replace the thumbnail
|
||||
shutil.copy('source/auto_examples/images/blank_image.png', thumb_file)
|
||||
|
||||
docstring, short_desc, end_row = extract_docstring(example_file)
|
||||
|
||||
# Depending on whether we have one or more figures, we're using a
|
||||
# horizontal list or a single rst call to 'image'.
|
||||
if len(figure_list) == 1:
|
||||
figure_name = figure_list[0]
|
||||
image_list = SINGLE_IMAGE % figure_name.lstrip('/')
|
||||
else:
|
||||
image_list = HLIST_HEADER
|
||||
for figure_name in figure_list:
|
||||
image_list += HLIST_IMAGE_TEMPLATE % figure_name.lstrip('/')
|
||||
|
||||
f = open(os.path.join(target_dir, fname[:-2] + 'txt'),'w')
|
||||
f.write(this_template % locals())
|
||||
f.flush()
|
||||
|
||||
|
||||
def setup(app):
|
||||
app.connect('builder-inited', generate_example_rst)
|
||||
app.add_config_value('plot_gallery', True, 'html')
|
||||
@@ -0,0 +1,502 @@
|
||||
"""
|
||||
Example generation from python files.
|
||||
|
||||
Generate the rst files for the examples by iterating over the python
|
||||
example files. Files that generate images should start with 'plot'.
|
||||
|
||||
To generate your own examples, add this extension to the list of
|
||||
``extensions``in your Sphinx configuration file. In addition, make sure the
|
||||
example directory(ies) in `plot2rst_paths` (see below) points to a directory
|
||||
with examples named `plot_*.py` and include an `index.rst` file.
|
||||
|
||||
This code was adapted from scikits-image, which took it from scikits-learn.
|
||||
|
||||
Options
|
||||
-------
|
||||
The ``plot2rst`` extension accepts the following options:
|
||||
|
||||
plot2rst_paths : length-2 tuple, or list of tuples
|
||||
Tuple or list of tuples of paths to (python plot, generated rst) files,
|
||||
i.e. (source, destination). Note that both paths are relative to Sphinx
|
||||
'source' directory. Defaults to ('../examples', 'auto_examples')
|
||||
|
||||
plot2rst_rcparams : dict
|
||||
Matplotlib configuration parameters. See
|
||||
http://matplotlib.sourceforge.net/users/customizing.html for details.
|
||||
|
||||
plot2rst_default_thumb : str
|
||||
Path (relative to doc root) of default thumbnail image.
|
||||
|
||||
plot2rst_thumb_scale : float
|
||||
Scale factor for thumbnail (e.g., 0.2 to scale plot to 1/5th the
|
||||
original size).
|
||||
|
||||
plot2rst_plot_tag : str
|
||||
When this tag is found in the example file, the current plot is saved and
|
||||
tag is replaced with plot path. Defaults to 'PLOT2RST.current_figure'.
|
||||
|
||||
|
||||
Suggested CSS definitions
|
||||
-------------------------
|
||||
|
||||
div.body h2 {
|
||||
border-bottom: 1px solid #BBB;
|
||||
clear: left;
|
||||
}
|
||||
|
||||
/*---- example gallery ----*/
|
||||
|
||||
.gallery.figure {
|
||||
float: left;
|
||||
margin: 1em;
|
||||
}
|
||||
|
||||
.gallery.figure img{
|
||||
display: block;
|
||||
margin-left: auto;
|
||||
margin-right: auto;
|
||||
width: 200px;
|
||||
}
|
||||
|
||||
.gallery.figure .caption {
|
||||
width: 200px;
|
||||
text-align: center !important;
|
||||
}
|
||||
|
||||
"""
|
||||
import os
|
||||
import shutil
|
||||
import token
|
||||
import tokenize
|
||||
|
||||
import numpy as np
|
||||
import matplotlib
|
||||
matplotlib.use('Agg')
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib import image
|
||||
|
||||
|
||||
LITERALINCLUDE = """
|
||||
.. literalinclude:: {src_name}
|
||||
:lines: {code_start}-
|
||||
|
||||
"""
|
||||
|
||||
CODE_LINK = """
|
||||
|
||||
**Python source code:** :download:`download <{0}>`
|
||||
(generated using ``skimage`` |version|)
|
||||
|
||||
"""
|
||||
|
||||
TOCTREE_TEMPLATE = """
|
||||
.. toctree::
|
||||
:hidden:
|
||||
|
||||
%s
|
||||
|
||||
"""
|
||||
|
||||
IMAGE_TEMPLATE = """
|
||||
.. image:: images/%s
|
||||
:align: center
|
||||
|
||||
"""
|
||||
|
||||
GALLERY_IMAGE_TEMPLATE = """
|
||||
.. figure:: %(thumb)s
|
||||
:figclass: gallery
|
||||
:target: ./%(source)s.html
|
||||
|
||||
:ref:`example_%(link_name)s`
|
||||
|
||||
"""
|
||||
|
||||
|
||||
class Path(str):
|
||||
"""Path object for manipulating directory and file paths."""
|
||||
|
||||
def __init__(self, path):
|
||||
super(Path, self).__init__(path)
|
||||
|
||||
@property
|
||||
def isdir(self):
|
||||
return os.path.isdir(self)
|
||||
|
||||
@property
|
||||
def exists(self):
|
||||
"""Return True if path exists"""
|
||||
return os.path.exists(self)
|
||||
|
||||
def pjoin(self, *args):
|
||||
"""Join paths. `p` prefix prevents confusion with string method."""
|
||||
return self.__class__(os.path.join(self, *args))
|
||||
|
||||
def psplit(self):
|
||||
"""Split paths. `p` prefix prevents confusion with string method."""
|
||||
return [self.__class__(p) for p in os.path.split(self)]
|
||||
|
||||
def makedirs(self):
|
||||
if not self.exists:
|
||||
os.makedirs(self)
|
||||
|
||||
def listdir(self):
|
||||
return os.listdir(self)
|
||||
|
||||
def format(self, *args, **kwargs):
|
||||
return self.__class__(super(Path, self).format(*args, **kwargs))
|
||||
|
||||
def __add__(self, other):
|
||||
return self.__class__(super(Path, self).__add__(other))
|
||||
|
||||
def __iadd__(self, other):
|
||||
return self.__add__(other)
|
||||
|
||||
|
||||
def setup(app):
|
||||
app.connect('builder-inited', generate_example_galleries)
|
||||
|
||||
app.add_config_value('plot2rst_paths',
|
||||
('../examples', 'auto_examples'), True)
|
||||
app.add_config_value('plot2rst_rcparams', {}, True)
|
||||
app.add_config_value('plot2rst_default_thumb', None, True)
|
||||
app.add_config_value('plot2rst_thumb_scale', 0.25, True)
|
||||
app.add_config_value('plot2rst_plot_tag', 'PLOT2RST.current_figure', True)
|
||||
app.add_config_value('plot2rst_index_name', 'index', True)
|
||||
|
||||
|
||||
def generate_example_galleries(app):
|
||||
cfg = app.builder.config
|
||||
|
||||
doc_src = Path(os.path.abspath(app.builder.srcdir)) # path/to/doc/source
|
||||
|
||||
if isinstance(cfg.plot2rst_paths, tuple):
|
||||
cfg.plot2rst_paths = [cfg.plot2rst_paths]
|
||||
for src_dest in cfg.plot2rst_paths:
|
||||
plot_path, rst_path = [Path(p) for p in src_dest]
|
||||
example_dir = doc_src.pjoin(plot_path)
|
||||
rst_dir = doc_src.pjoin(rst_path)
|
||||
generate_examples_and_gallery(example_dir, rst_dir, cfg)
|
||||
|
||||
|
||||
def generate_examples_and_gallery(example_dir, rst_dir, cfg):
|
||||
"""Generate rst from examples and create gallery to showcase examples."""
|
||||
if not example_dir.exists:
|
||||
print "No example directory found at", example_dir
|
||||
return
|
||||
rst_dir.makedirs()
|
||||
|
||||
# we create an index.rst with all examples
|
||||
gallery_index = file(rst_dir.pjoin('index'+cfg.source_suffix), 'w')
|
||||
|
||||
# Here we don't use an os.walk, but we recurse only twice: flat is
|
||||
# better than nested.
|
||||
write_gallery(gallery_index, example_dir, rst_dir, cfg)
|
||||
for d in sorted(example_dir.listdir()):
|
||||
example_sub = example_dir.pjoin(d)
|
||||
if example_sub.isdir:
|
||||
rst_sub = rst_dir.pjoin(d)
|
||||
rst_sub.makedirs()
|
||||
write_gallery(gallery_index, example_sub, rst_sub, cfg, depth=1)
|
||||
gallery_index.flush()
|
||||
|
||||
|
||||
def write_gallery(gallery_index, src_dir, rst_dir, cfg, depth=0):
|
||||
"""Generate the rst files for an example directory, i.e. gallery.
|
||||
|
||||
Write rst files from python examples and add example links to gallery.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
gallery_index : file
|
||||
Index file for plot gallery.
|
||||
src_dir : 'str'
|
||||
Source directory for python examples.
|
||||
rst_dir : 'str'
|
||||
Destination directory for rst files generated from python examples.
|
||||
cfg : config object
|
||||
Sphinx config object created by Sphinx.
|
||||
"""
|
||||
index_name = cfg.plot2rst_index_name + cfg.source_suffix
|
||||
gallery_template = src_dir.pjoin(index_name)
|
||||
if not os.path.exists(gallery_template):
|
||||
print src_dir
|
||||
print 80*'_'
|
||||
print ('Example directory %s does not have a %s file'
|
||||
% (src_dir, index_name))
|
||||
print 'Skipping this directory'
|
||||
print 80*'_'
|
||||
return
|
||||
|
||||
gallery_description = file(gallery_template).read()
|
||||
gallery_index.write('\n\n%s\n\n' % gallery_description)
|
||||
|
||||
rst_dir.makedirs()
|
||||
examples = [fname for fname in sorted(src_dir.listdir(), key=_plots_first)
|
||||
if fname.endswith('py')]
|
||||
ex_names = [ex[:-3] for ex in examples] # strip '.py' extension
|
||||
if depth == 0:
|
||||
sub_dir = Path('')
|
||||
else:
|
||||
sub_dir_list = src_dir.psplit()[-depth:]
|
||||
sub_dir = Path('/'.join(sub_dir_list) + '/')
|
||||
gallery_index.write(TOCTREE_TEMPLATE % (sub_dir + '\n '.join(ex_names)))
|
||||
|
||||
for src_name in examples:
|
||||
write_example(src_name, src_dir, rst_dir, cfg)
|
||||
|
||||
link_name = sub_dir.pjoin(src_name)
|
||||
link_name = link_name.replace(os.path.sep, '_')
|
||||
if link_name.startswith('._'):
|
||||
link_name = link_name[2:]
|
||||
|
||||
info = {}
|
||||
info['thumb'] = sub_dir.pjoin('images/thumb', src_name[:-3] + '.png')
|
||||
info['source'] = sub_dir + src_name[:-3]
|
||||
info['link_name'] = link_name
|
||||
gallery_index.write(GALLERY_IMAGE_TEMPLATE % info)
|
||||
|
||||
|
||||
def _plots_first(fname):
|
||||
"""Decorate filename so that examples with plots are displayed first."""
|
||||
if not (fname.startswith('plot') and fname.endswith('.py')):
|
||||
return 'zz' + fname
|
||||
return fname
|
||||
|
||||
|
||||
def write_example(src_name, src_dir, rst_dir, cfg):
|
||||
"""Write rst file from a given python example.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
src_name : str
|
||||
Name of example file.
|
||||
src_dir : 'str'
|
||||
Source directory for python examples.
|
||||
rst_dir : 'str'
|
||||
Destination directory for rst files generated from python examples.
|
||||
cfg : config object
|
||||
Sphinx config object created by Sphinx.
|
||||
"""
|
||||
last_dir = src_dir.psplit()[-1]
|
||||
# to avoid leading . in file names, and wrong names in links
|
||||
if last_dir == '.' or last_dir == 'examples':
|
||||
last_dir = Path('')
|
||||
else:
|
||||
last_dir += '_'
|
||||
|
||||
src_path = src_dir.pjoin(src_name)
|
||||
example_file = rst_dir.pjoin(src_name)
|
||||
shutil.copyfile(src_path, example_file)
|
||||
|
||||
image_dir = rst_dir.pjoin('images')
|
||||
thumb_dir = image_dir.pjoin('thumb')
|
||||
image_dir.makedirs()
|
||||
thumb_dir.makedirs()
|
||||
|
||||
base_image_name = os.path.splitext(src_name)[0]
|
||||
image_path = image_dir.pjoin(base_image_name + '_{0}.png')
|
||||
|
||||
basename, py_ext = os.path.splitext(src_name)
|
||||
rst_path = rst_dir.pjoin(basename + cfg.source_suffix)
|
||||
|
||||
if _plots_are_current(src_path, image_path) and rst_path.exists:
|
||||
return
|
||||
|
||||
blocks = split_code_and_text_blocks(example_file)
|
||||
if blocks[0][2].startswith('#!'):
|
||||
blocks.pop(0) # don't add shebang line to rst file.
|
||||
|
||||
rst_link = '.. _example_%s:\n\n' % (last_dir + src_name)
|
||||
figure_list, rst = process_blocks(blocks, src_path, image_path, cfg)
|
||||
|
||||
has_inline_plots = any(cfg.plot2rst_plot_tag in b[2] for b in blocks)
|
||||
if has_inline_plots:
|
||||
example_rst = ''.join([rst_link, rst])
|
||||
else:
|
||||
# print first block of text, display all plots, then display code.
|
||||
first_text_block = [b for b in blocks if b[0] == 'text'][0]
|
||||
label, (start, end), content = first_text_block
|
||||
figure_list = save_all_figures(image_path)
|
||||
rst_blocks = [IMAGE_TEMPLATE % f.lstrip('/') for f in figure_list]
|
||||
|
||||
example_rst = rst_link
|
||||
example_rst += eval(content)
|
||||
example_rst += ''.join(rst_blocks)
|
||||
code_info = dict(src_name=src_name, code_start=end)
|
||||
example_rst += LITERALINCLUDE.format(**code_info)
|
||||
|
||||
example_rst += CODE_LINK.format(src_name)
|
||||
|
||||
f = open(rst_path,'w')
|
||||
f.write(example_rst)
|
||||
f.flush()
|
||||
|
||||
thumb_path = thumb_dir.pjoin(src_name[:-3] + '.png')
|
||||
first_image_file = image_dir.pjoin(figure_list[0].lstrip('/'))
|
||||
if first_image_file.exists:
|
||||
image.thumbnail(first_image_file, thumb_path, cfg.plot2rst_thumb_scale)
|
||||
|
||||
if not thumb_path.exists:
|
||||
if cfg.plot2rst_default_thumb is None:
|
||||
print "WARNING: No plots found and default thumbnail not defined."
|
||||
print "Specify 'plot2rst_default_thumb' in Sphinx config file."
|
||||
else:
|
||||
shutil.copy(cfg.plot2rst_default_thumb, thumb_path)
|
||||
|
||||
|
||||
def _plots_are_current(src_path, image_path):
|
||||
first_image_file = Path(image_path.format(1))
|
||||
needs_replot = (not first_image_file.exists or
|
||||
_mod_time(first_image_file) <= _mod_time(src_path))
|
||||
return not needs_replot
|
||||
|
||||
|
||||
def _mod_time(file_path):
|
||||
return os.stat(file_path).st_mtime
|
||||
|
||||
|
||||
def split_code_and_text_blocks(source_file):
|
||||
"""Return list with source file separated into code and text blocks.
|
||||
|
||||
Returns
|
||||
-------
|
||||
blocks : list of (label, (start, end+1), content)
|
||||
List where each element is a tuple with the label ('text' or 'code'),
|
||||
the (start, end+1) line numbers, and content string of block.
|
||||
"""
|
||||
block_edges, idx_first_text_block = get_block_edges(source_file)
|
||||
|
||||
with open(source_file) as f:
|
||||
source_lines = f.readlines()
|
||||
|
||||
# Every other block should be a text block
|
||||
idx_text_block = np.arange(idx_first_text_block, len(block_edges), 2)
|
||||
blocks = []
|
||||
slice_ranges = zip(block_edges[:-1], block_edges[1:])
|
||||
for i, (start, end) in enumerate(slice_ranges):
|
||||
block_label = 'text' if i in idx_text_block else 'code'
|
||||
# subtract 1 from indices b/c line numbers start at 1, not 0
|
||||
content = ''.join(source_lines[start-1:end-1])
|
||||
blocks.append((block_label, (start, end), content))
|
||||
return blocks
|
||||
|
||||
|
||||
def get_block_edges(source_file):
|
||||
"""Return starting line numbers of code and text blocks
|
||||
|
||||
Returns
|
||||
-------
|
||||
block_edges : list of int
|
||||
Line number for the start of each block. Note the
|
||||
idx_first_text_block : {0 | 1}
|
||||
0 if first block is text then, else 1 (second block better be text).
|
||||
"""
|
||||
block_edges = []
|
||||
with open(source_file) as f:
|
||||
token_iter = tokenize.generate_tokens(f.readline)
|
||||
for token_tuple in token_iter:
|
||||
t_id, t_str, (srow, scol), (erow, ecol), src_line = token_tuple
|
||||
if (token.tok_name[t_id] == 'STRING' and scol == 0):
|
||||
# Add one point to line after text (for later slicing)
|
||||
block_edges.extend((srow, erow+1))
|
||||
idx_first_text_block = 0
|
||||
# when example doesn't start with text block.
|
||||
if not block_edges[0] == 1:
|
||||
block_edges.insert(0, 1)
|
||||
idx_first_text_block = 1
|
||||
# when example doesn't end with text block.
|
||||
if not block_edges[-1] == erow: # iffy: I'm using end state of loop
|
||||
block_edges.append(erow)
|
||||
return block_edges, idx_first_text_block
|
||||
|
||||
|
||||
def process_blocks(blocks, src_path, image_path, cfg):
|
||||
"""Run source, save plots as images, and convert blocks to rst.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
blocks : list of block tuples
|
||||
Code and text blocks from example. See `split_code_and_text_blocks`.
|
||||
src_path : str
|
||||
Path to example file.
|
||||
image_path : str
|
||||
Path where plots are saved (format string which accepts figure number).
|
||||
cfg : config object
|
||||
Sphinx config object created by Sphinx.
|
||||
|
||||
Returns
|
||||
-------
|
||||
figure_list : list
|
||||
List of figure names saved by the example.
|
||||
rst_text : str
|
||||
Text with code wrapped code-block directives.
|
||||
"""
|
||||
src_dir, src_name = src_path.psplit()
|
||||
if not src_name.startswith('plot'):
|
||||
return [], ''
|
||||
|
||||
# index of blocks which have inline plots
|
||||
inline_tag = cfg.plot2rst_plot_tag
|
||||
idx_inline_plot = [i for i, b in enumerate(blocks)
|
||||
if inline_tag in b[2]]
|
||||
|
||||
image_dir, image_fmt_str = image_path.psplit()
|
||||
|
||||
figure_list = []
|
||||
plt.rcdefaults()
|
||||
plt.rcParams.update(cfg.plot2rst_rcparams)
|
||||
plt.close('all')
|
||||
|
||||
example_globals = {}
|
||||
rst_blocks = []
|
||||
fig_num = 1
|
||||
for i, (blabel, brange, bcontent) in enumerate(blocks):
|
||||
if blabel == 'code':
|
||||
exec(bcontent, example_globals)
|
||||
rst_blocks.append(codestr2rst(bcontent))
|
||||
else:
|
||||
if i in idx_inline_plot:
|
||||
plt.savefig(image_path.format(fig_num))
|
||||
figure_name = image_fmt_str.format(fig_num)
|
||||
fig_num += 1
|
||||
figure_list.append(figure_name)
|
||||
figure_link = os.path.join('images', figure_name)
|
||||
bcontent = bcontent.replace(inline_tag, figure_link)
|
||||
rst_blocks.append(docstr2rst(bcontent))
|
||||
return figure_list, '\n'.join(rst_blocks)
|
||||
|
||||
|
||||
def codestr2rst(codestr):
|
||||
"""Return reStructuredText code block from code string"""
|
||||
code_directive = ".. code-block:: python\n\n"
|
||||
indented_block = '\t' + codestr.replace('\n', '\n\t')
|
||||
return code_directive + indented_block
|
||||
|
||||
|
||||
def docstr2rst(docstr):
|
||||
"""Return reStructuredText from docstring"""
|
||||
idx_whitespace = len(docstr.rstrip()) - len(docstr)
|
||||
whitespace = docstr[idx_whitespace:]
|
||||
return eval(docstr) + whitespace
|
||||
|
||||
|
||||
def save_all_figures(image_path):
|
||||
"""Save all matplotlib figures.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image_path : str
|
||||
Path where plots are saved (format string which accepts figure number).
|
||||
"""
|
||||
figure_list = []
|
||||
image_dir, image_fmt_str = image_path.psplit()
|
||||
fig_mngr = matplotlib._pylab_helpers.Gcf.get_all_fig_managers()
|
||||
for fig_num in (m.num for m in fig_mngr):
|
||||
# Set the fig_num figure as the current figure as we can't
|
||||
# save a figure that's not the current figure.
|
||||
plt.figure(fig_num)
|
||||
plt.savefig(image_path.format(fig_num))
|
||||
figure_list.append(image_fmt_str.format(fig_num))
|
||||
return figure_list
|
||||
|
||||
@@ -0,0 +1,40 @@
|
||||
Announcement: scikits-image 0.6
|
||||
===============================
|
||||
|
||||
We're happy to announce the 6th version of scikits-image!
|
||||
|
||||
Scikits-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://skimage.org
|
||||
|
||||
New Features
|
||||
------------
|
||||
- Packaged in Debian as ``python-skimage``
|
||||
- Template matching
|
||||
- Fast user-defined image warping
|
||||
- Adaptive thresholding
|
||||
- Structural similarity index
|
||||
- Polygon, circle and ellipse drawing
|
||||
- Peak detection
|
||||
- Region properties
|
||||
- TiffFile I/O plugin
|
||||
|
||||
... along with some bug fixes and performance tweaks.
|
||||
|
||||
Contributors to this release
|
||||
----------------------------
|
||||
- Vincent Albufera
|
||||
- David Cournapeau
|
||||
- Christoph Gohlke
|
||||
- Emmanuelle Gouillart
|
||||
- Pieter Holtzhausen
|
||||
- Zachary Pincus
|
||||
- Johannes Schönberger
|
||||
- Tom (tangofoxtrotmike)
|
||||
- James Turner
|
||||
- Stefan van der Walt
|
||||
- Tony S Yu
|
||||
+6
-8
@@ -23,17 +23,10 @@ sys.path.append(os.path.join(curpath, '..', 'ext'))
|
||||
|
||||
# -- General configuration -----------------------------------------------------
|
||||
|
||||
try:
|
||||
import gen_rst
|
||||
except:
|
||||
pass
|
||||
|
||||
|
||||
# Add any Sphinx extension module names here, as strings. They can be extensions
|
||||
# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
|
||||
extensions = ['sphinx.ext.autodoc', 'sphinx.ext.pngmath', 'numpydoc',
|
||||
'sphinx.ext.autosummary', 'sphinx.ext.inheritance_diagram',
|
||||
'plot_directive', 'gen_rst']
|
||||
'sphinx.ext.autosummary', 'plot_directive', 'plot2rst']
|
||||
|
||||
# Add any paths that contain templates here, relative to this directory.
|
||||
templates_path = ['_templates']
|
||||
@@ -258,3 +251,8 @@ matplotlib.rcParams.update({
|
||||
"""
|
||||
plot_include_source = True
|
||||
plot_formats = [('png', 100)]
|
||||
|
||||
plot2rst_index_name = 'README'
|
||||
plot2rst_rcparams = {'image.cmap' : 'gray',
|
||||
'image.interpolation' : 'none'}
|
||||
|
||||
|
||||
@@ -79,6 +79,7 @@ h1, h2, h3, h4 {
|
||||
font-weight: normal;
|
||||
color: {{ theme_headercolor2 }};
|
||||
margin-bottom: .8em;
|
||||
clear: left;
|
||||
}
|
||||
|
||||
h1 {
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
function insert_version_links() {
|
||||
var labels = ['dev', '0.5', '0.4', '0.3'];
|
||||
var labels = ['dev', '0.6', '0.5', '0.4', '0.3'];
|
||||
|
||||
document.write('<ul class="versions">\n');
|
||||
|
||||
|
||||
@@ -17,7 +17,7 @@ MAINTAINER_EMAIL = 'stefan@sun.ac.za'
|
||||
URL = 'http://scikits-image.org'
|
||||
LICENSE = 'Modified BSD'
|
||||
DOWNLOAD_URL = 'http://github.com/scikits-image/scikits-image'
|
||||
VERSION = '0.6dev'
|
||||
VERSION = '0.6'
|
||||
|
||||
import os
|
||||
import setuptools
|
||||
|
||||
+5
-3
@@ -38,15 +38,17 @@ def cython(pyx_files, working_path=''):
|
||||
cmd = 'cython -o %s %s' % (c_file, pyxfile)
|
||||
print(cmd)
|
||||
|
||||
if platform.system() == 'Windows':
|
||||
try:
|
||||
status = subprocess.call(['cython', '-o', c_file, pyxfile])
|
||||
except WindowsError:
|
||||
# On Windows cython.exe may be missing if Cython was installed
|
||||
# via distutils. Run the cython.py script instead.
|
||||
status = subprocess.call(
|
||||
[sys.executable,
|
||||
os.path.join(os.path.dirname(sys.executable),
|
||||
'Scripts', 'cython.py'),
|
||||
'-o', c_file, pyxfile],
|
||||
shell=True)
|
||||
else:
|
||||
status = subprocess.call(['cython', '-o', c_file, pyxfile])
|
||||
|
||||
def _md5sum(f):
|
||||
m = hashlib.new('md5')
|
||||
|
||||
@@ -174,6 +174,8 @@ def rescale_intensity(image, in_range=None, out_range=None):
|
||||
|
||||
if out_range is None:
|
||||
omin, omax = dtype_range[dtype]
|
||||
if imin >= 0:
|
||||
omin = 0
|
||||
else:
|
||||
omin, omax = out_range
|
||||
|
||||
|
||||
@@ -2,3 +2,4 @@ from .hog import hog
|
||||
from .greycomatrix import greycomatrix, greycoprops
|
||||
from .peak import peak_local_max
|
||||
from .harris import harris
|
||||
from .template import match_template
|
||||
|
||||
@@ -0,0 +1,129 @@
|
||||
"""
|
||||
Template matching using normalized cross-correlation.
|
||||
|
||||
We use fast normalized cross-correlation algorithm (see [1]_ and [2]_) to
|
||||
compute match probability. This algorithm calculates the normalized
|
||||
cross-correlation of an image, `I`, with a template `T` according to the
|
||||
following equation::
|
||||
|
||||
sum{ I(x, y) [T(x, y) - <T>] }
|
||||
-------------------------------------------------------
|
||||
sqrt(sum{ [I(x, y) - <I>]^2 } sum{ [T(x, y) - <T>]^2 })
|
||||
|
||||
where `<T>` is the average of the template, and `<I>` is the average of the
|
||||
image *coincident with the template*, and sums are over the template and the
|
||||
image window coincident with the template. Note that the numerator is simply
|
||||
the cross-correlation of the image and the zero-mean template.
|
||||
|
||||
To speed up calculations, we use summed-area tables (a.k.a. integral images) to
|
||||
quickly calculate sums of image windows inside the loop. This step relies on
|
||||
the following relation (see Eq. 10 of [1])::
|
||||
|
||||
sum{ [I(x, y) - <I>]^2 } =
|
||||
sum{ I^2(x, y) } - [sum{ I(x, y) }]^2 / N_x N_y
|
||||
|
||||
(Without this relation, you would need to subtract each image-window mean from
|
||||
the image window *before* squaring.)
|
||||
|
||||
.. [1] Briechle and Hanebeck, "Template Matching using Fast Normalized
|
||||
Cross Correlation", Proceedings of the SPIE (2001).
|
||||
.. [2] J. P. Lewis, "Fast Normalized Cross-Correlation", Industrial Light and
|
||||
Magic.
|
||||
"""
|
||||
import cython
|
||||
cimport numpy as np
|
||||
import numpy as np
|
||||
from scipy.signal import fftconvolve
|
||||
from skimage.transform import integral
|
||||
|
||||
|
||||
cdef extern from "math.h":
|
||||
float sqrt(float x)
|
||||
float fabs(float x)
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
cdef float integrate(np.ndarray[float, ndim=2, mode="c"] sat,
|
||||
int r0, int c0, int r1, int c1):
|
||||
"""
|
||||
Using a summed area table / integral image, calculate the sum
|
||||
over a given window.
|
||||
|
||||
This function is the same as the `integrate` function in
|
||||
`skimage.transform.integrate`, but this Cython version significantly
|
||||
speeds up the code.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
sat : ndarray of float
|
||||
Summed area table / integral image.
|
||||
r0, c0 : int
|
||||
Top-left corner of block to be summed.
|
||||
r1, c1 : int
|
||||
Bottom-right corner of block to be summed.
|
||||
|
||||
Returns
|
||||
-------
|
||||
S : int
|
||||
Sum over the given window.
|
||||
"""
|
||||
cdef float S = 0
|
||||
|
||||
S += sat[r1, c1]
|
||||
|
||||
if (r0 - 1 >= 0) and (c0 - 1 >= 0):
|
||||
S += sat[r0 - 1, c0 - 1]
|
||||
|
||||
if (r0 - 1 >= 0):
|
||||
S -= sat[r0 - 1, c1]
|
||||
|
||||
if (c0 - 1 >= 0):
|
||||
S -= sat[r1, c0 - 1]
|
||||
return S
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
def match_template(np.ndarray[float, ndim=2, mode="c"] image,
|
||||
np.ndarray[float, ndim=2, mode="c"] template):
|
||||
cdef np.ndarray[float, ndim=2, mode="c"] corr
|
||||
cdef np.ndarray[float, ndim=2, mode="c"] image_sat
|
||||
cdef np.ndarray[float, ndim=2, mode="c"] image_sqr_sat
|
||||
cdef float template_mean = np.mean(template)
|
||||
cdef float template_ssd
|
||||
cdef float inv_area
|
||||
|
||||
image_sat = integral.integral_image(image)
|
||||
image_sqr_sat = integral.integral_image(image**2)
|
||||
|
||||
template -= template_mean
|
||||
template_ssd = np.sum(template**2)
|
||||
# use inversed area for accuracy
|
||||
inv_area = 1.0 / (template.shape[0] * template.shape[1])
|
||||
|
||||
# when `dtype=float` is used, ascontiguousarray returns ``double``.
|
||||
corr = np.ascontiguousarray(fftconvolve(image,
|
||||
template[::-1, ::-1],
|
||||
mode="valid"),
|
||||
dtype=np.float32)
|
||||
|
||||
cdef int i, j
|
||||
cdef float den, window_sqr_sum, window_mean_sqr, window_sum,
|
||||
# move window through convolution results, normalizing in the process
|
||||
for i in range(corr.shape[0]):
|
||||
for j in range(corr.shape[1]):
|
||||
# subtract 1 because `i_end` and `j_end` are used for indexing into
|
||||
# summed-area table, instead of slicing windows of the image.
|
||||
i_end = i + template.shape[0] - 1
|
||||
j_end = j + template.shape[1] - 1
|
||||
|
||||
window_sum = integrate(image_sat, i, j, i_end, j_end)
|
||||
window_mean_sqr = window_sum * window_sum * inv_area
|
||||
window_sqr_sum = integrate(image_sqr_sat, i, j, i_end, j_end)
|
||||
if window_sqr_sum <= window_mean_sqr:
|
||||
corr[i, j] = 0
|
||||
continue
|
||||
|
||||
den = sqrt((window_sqr_sum - window_mean_sqr) * template_ssd)
|
||||
corr[i, j] /= den
|
||||
return corr
|
||||
|
||||
@@ -97,7 +97,7 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8),
|
||||
magnitude = sqrt(gx ** 2 + gy ** 2)
|
||||
orientation = arctan2(gy, (gx + 1e-15)) * (180 / pi) + 90
|
||||
|
||||
sx, sy = image.shape
|
||||
sy, sx = image.shape
|
||||
cx, cy = pixels_per_cell
|
||||
bx, by = cells_per_block
|
||||
|
||||
@@ -105,7 +105,7 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8),
|
||||
n_cellsy = int(np.floor(sy // cy)) # number of cells in y
|
||||
|
||||
# compute orientations integral images
|
||||
orientation_histogram = np.zeros((n_cellsx, n_cellsy, orientations))
|
||||
orientation_histogram = np.zeros((n_cellsy, n_cellsx, orientations))
|
||||
for i in range(orientations):
|
||||
#create new integral image for this orientation
|
||||
# isolate orientations in this range
|
||||
@@ -118,7 +118,7 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8),
|
||||
cond2 = temp_ori > 0
|
||||
temp_mag = np.where(cond2, magnitude, 0)
|
||||
|
||||
orientation_histogram[:,:,i] = uniform_filter(temp_mag, size=(cx, cy))[cx/2::cx, cy/2::cy].T
|
||||
orientation_histogram[:,:,i] = uniform_filter(temp_mag, size=(cy, cx))[cy/2::cy, cx/2::cx]
|
||||
|
||||
|
||||
# now for each cell, compute the histogram
|
||||
@@ -140,7 +140,7 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8),
|
||||
dy = radius * sin(float(o) / orientations * np.pi)
|
||||
rr, cc = draw.bresenham(centre[0] - dx, centre[1] - dy,
|
||||
centre[0] + dx, centre[1] + dy)
|
||||
hog_image[rr, cc] += orientation_histogram[x, y, o]
|
||||
hog_image[rr, cc] += orientation_histogram[y, x, o]
|
||||
|
||||
"""
|
||||
The fourth stage computes normalisation, which takes local groups of
|
||||
@@ -159,14 +159,14 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8),
|
||||
|
||||
n_blocksx = (n_cellsx - bx) + 1
|
||||
n_blocksy = (n_cellsy - by) + 1
|
||||
normalised_blocks = np.zeros((n_blocksx, n_blocksy,
|
||||
bx, by, orientations))
|
||||
normalised_blocks = np.zeros((n_blocksy, n_blocksx,
|
||||
by, bx, orientations))
|
||||
|
||||
for x in range(n_blocksx):
|
||||
for y in range(n_blocksy):
|
||||
block = orientation_histogram[x:x + bx, y:y + by, :]
|
||||
block = orientation_histogram[y:y + by, x:x + bx, :]
|
||||
eps = 1e-5
|
||||
normalised_blocks[x, y, :] = block / sqrt(block.sum() ** 2 + eps)
|
||||
normalised_blocks[y, x, :] = block / sqrt(block.sum() ** 2 + eps)
|
||||
|
||||
"""
|
||||
The final step collects the HOG descriptors from all blocks of a dense
|
||||
|
||||
@@ -12,9 +12,12 @@ def configuration(parent_package='', top_path=None):
|
||||
config.add_data_dir('tests')
|
||||
|
||||
cython(['_greycomatrix.pyx'], working_path=base_path)
|
||||
cython(['_template.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('_greycomatrix', sources=['_greycomatrix.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_template', sources=['_template.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
|
||||
@@ -0,0 +1,84 @@
|
||||
"""template.py - Template matching
|
||||
"""
|
||||
import numpy as np
|
||||
from . import _template
|
||||
|
||||
from skimage.util.dtype import convert
|
||||
|
||||
|
||||
def match_template(image, template, pad_input=False):
|
||||
"""Match a template to an image using normalized correlation.
|
||||
|
||||
The output is an array with values between -1.0 and 1.0, which correspond
|
||||
to the probability that the template is found at that position.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : array_like
|
||||
Image to process.
|
||||
template : array_like
|
||||
Template to locate.
|
||||
pad_input : bool
|
||||
If True, pad `image` with image mean so that output is the same size as
|
||||
the image, and output values correspond to the template center.
|
||||
Otherwise, the output is an array with shape `(M - m + 1, N - n + 1)`
|
||||
for an `(M, N)` image and an `(m, n)` template, and matches correspond
|
||||
to origin (top-left corner) of the template.
|
||||
|
||||
Returns
|
||||
-------
|
||||
output : ndarray
|
||||
Correlation results between -1.0 and 1.0. For an `(M, N)` image and an
|
||||
`(m, n)` template, the `output` is `(M - m + 1, N - n + 1)` when
|
||||
`pad_input = False` and `(M, N)` when `pad_input = True`.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> template = np.zeros((3, 3))
|
||||
>>> template[1, 1] = 1
|
||||
>>> print template
|
||||
[[ 0. 0. 0.]
|
||||
[ 0. 1. 0.]
|
||||
[ 0. 0. 0.]]
|
||||
>>> image = np.zeros((6, 6))
|
||||
>>> image[1, 1] = 1
|
||||
>>> image[4, 4] = -1
|
||||
>>> print image
|
||||
[[ 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. 0. 0. -1. 0.]
|
||||
[ 0. 0. 0. 0. 0. 0.]]
|
||||
>>> result = match_template(image, template)
|
||||
>>> print np.round(result, 3)
|
||||
[[ 1. -0.125 0. 0. ]
|
||||
[-0.125 -0.125 0. 0. ]
|
||||
[ 0. 0. 0.125 0.125]
|
||||
[ 0. 0. 0.125 -1. ]]
|
||||
>>> result = match_template(image, template, pad_input=True)
|
||||
>>> print np.round(result, 3)
|
||||
[[-0.125 -0.125 -0.125 0. 0. 0. ]
|
||||
[-0.125 1. -0.125 0. 0. 0. ]
|
||||
[-0.125 -0.125 -0.125 0. 0. 0. ]
|
||||
[ 0. 0. 0. 0.125 0.125 0.125]
|
||||
[ 0. 0. 0. 0.125 -1. 0.125]
|
||||
[ 0. 0. 0. 0.125 0.125 0.125]]
|
||||
"""
|
||||
if np.any(np.less(image.shape, template.shape)):
|
||||
raise ValueError("Image must be larger than template.")
|
||||
image = convert(image, np.float32)
|
||||
template = convert(template, np.float32)
|
||||
|
||||
if pad_input:
|
||||
pad_size = tuple(np.array(image.shape) + np.array(template.shape) - 1)
|
||||
pad_image = np.mean(image) * np.ones(pad_size, dtype=np.float32)
|
||||
h, w = image.shape
|
||||
i0, j0 = template.shape
|
||||
i0 /= 2
|
||||
j0 /= 2
|
||||
pad_image[i0:i0+h, j0:j0+w] = image
|
||||
image = pad_image
|
||||
result = _template.match_template(image, template)
|
||||
return result
|
||||
|
||||
@@ -5,12 +5,12 @@ from skimage.feature import hog
|
||||
|
||||
def test_histogram_of_oriented_gradients():
|
||||
# Replace with skimage.data.lena() after merge
|
||||
img = scipy.misc.lena().astype(np.int8)
|
||||
img = scipy.misc.lena()[:256,:].astype(np.int8)
|
||||
|
||||
fd = hog(img, orientations=9, pixels_per_cell=(8, 8),
|
||||
cells_per_block=(1, 1))
|
||||
|
||||
assert len(fd) == 9 * (512//8) ** 2
|
||||
assert len(fd) == 9 * (256//8) * (512//8)
|
||||
|
||||
if __name__ == '__main__':
|
||||
from numpy.testing import run_module_suite
|
||||
|
||||
@@ -0,0 +1,124 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_array_almost_equal as assert_close
|
||||
|
||||
from skimage.morphology import diamond
|
||||
from skimage.feature import match_template, peak_local_max
|
||||
|
||||
|
||||
def test_template():
|
||||
size = 100
|
||||
# Float prefactors ensure that image range is between 0 and 1
|
||||
image = 0.5 * np.ones((400, 400))
|
||||
target = 0.1 * (np.tri(size) + np.tri(size)[::-1])
|
||||
target_positions = [(50, 50), (200, 200)]
|
||||
for x, y in target_positions:
|
||||
image[x:x + size, y:y + size] = target
|
||||
np.random.seed(1)
|
||||
image += 0.1 * np.random.uniform(size=(400, 400))
|
||||
|
||||
result = match_template(image, target)
|
||||
delta = 5
|
||||
|
||||
positions = peak_local_max(result, min_distance=delta)
|
||||
|
||||
if len(positions) > 2:
|
||||
# Keep the two maximum peaks.
|
||||
intensities = result[tuple(positions.T)]
|
||||
i_maxsort = np.argsort(intensities)[::-1]
|
||||
positions = positions[i_maxsort][:2]
|
||||
|
||||
# Sort so that order matches `target_positions`.
|
||||
positions = positions[np.argsort(positions[:, 0])]
|
||||
|
||||
for xy_target, xy in zip(target_positions, positions):
|
||||
yield assert_close, xy, xy_target
|
||||
|
||||
|
||||
def test_normalization():
|
||||
"""Test that `match_template` gives the correct normalization.
|
||||
|
||||
Normalization gives 1 for a perfect match and -1 for an inverted-match.
|
||||
This test adds positive and negative squares to a zero-array and matches
|
||||
the array with a positive template.
|
||||
"""
|
||||
n = 5
|
||||
N = 20
|
||||
ipos, jpos = (2, 3)
|
||||
ineg, jneg = (12, 11)
|
||||
image = 0.5 * np.ones((N, N))
|
||||
image[ipos:ipos + n, jpos:jpos + n] = 1
|
||||
image[ineg:ineg + n, jneg:jneg + n] = 0
|
||||
|
||||
# white square with a black border
|
||||
template = np.zeros((n+2, n+2))
|
||||
template[1:1+n, 1:1+n] = 1
|
||||
|
||||
result = match_template(image, template)
|
||||
|
||||
# get the max and min results.
|
||||
sorted_result = np.argsort(result.flat)
|
||||
iflat_min = sorted_result[0]
|
||||
iflat_max = sorted_result[-1]
|
||||
min_result = np.unravel_index(iflat_min, result.shape)
|
||||
max_result = np.unravel_index(iflat_max, result.shape)
|
||||
|
||||
# shift result by 1 because of template border
|
||||
assert np.all((np.array(min_result) + 1) == (ineg, jneg))
|
||||
assert np.all((np.array(max_result) + 1) == (ipos, jpos))
|
||||
|
||||
assert np.allclose(result.flat[iflat_min], -1)
|
||||
assert np.allclose(result.flat[iflat_max], 1)
|
||||
|
||||
|
||||
def test_no_nans():
|
||||
"""Test that `match_template` doesn't return NaN values.
|
||||
|
||||
When image values are only slightly different, floating-point errors can
|
||||
cause a subtraction inside of a square root to go negative (without an
|
||||
explicit check that was added to `match_template`).
|
||||
"""
|
||||
np.random.seed(1)
|
||||
image = 0.5 + 1e-9 * np.random.normal(size=(20, 20))
|
||||
template = np.ones((6, 6))
|
||||
template[:3, :] = 0
|
||||
result = match_template(image, template)
|
||||
assert not np.any(np.isnan(result))
|
||||
|
||||
|
||||
def test_switched_arguments():
|
||||
image = np.ones((5, 5))
|
||||
template = np.ones((3, 3))
|
||||
np.testing.assert_raises(ValueError, match_template, template, image)
|
||||
|
||||
|
||||
def test_pad_input():
|
||||
"""Test `match_template` when `pad_input=True`.
|
||||
|
||||
This test places two full templates (one with values lower than the image
|
||||
mean, the other higher) and two half templates, which are on the edges of
|
||||
the image. The two full templates should score the top (positive and
|
||||
negative) matches and the centers of the half templates should score 2nd.
|
||||
"""
|
||||
# Float prefactors ensure that image range is between 0 and 1
|
||||
template = 0.5 * diamond(2)
|
||||
image = 0.5 * np.ones((9, 19))
|
||||
mid = slice(2, 7)
|
||||
image[mid, :3] -= template[:, -3:] # half min template centered at 0
|
||||
image[mid, 4:9] += template # full max template centered at 6
|
||||
image[mid, -9:-4] -= template # full min template centered at 12
|
||||
image[mid, -3:] += template[:, :3] # half max template centered at 18
|
||||
|
||||
result = match_template(image, template, pad_input=True)
|
||||
|
||||
# get the max and min results.
|
||||
sorted_result = np.argsort(result.flat)
|
||||
i, j = np.unravel_index(sorted_result[:2], result.shape)
|
||||
assert_close(j, (12, 0))
|
||||
i, j = np.unravel_index(sorted_result[-2:], result.shape)
|
||||
assert_close(j, (18, 6))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
from numpy import testing
|
||||
testing.run_module_suite()
|
||||
|
||||
@@ -17,7 +17,8 @@ cdef class MCP:
|
||||
cdef object flat_costs
|
||||
cdef object flat_cumulative_costs
|
||||
cdef object traceback_offsets
|
||||
cdef object flat_edge_map
|
||||
cdef object flat_pos_edge_map
|
||||
cdef object flat_neg_edge_map
|
||||
cdef readonly object offsets
|
||||
cdef object flat_offsets
|
||||
cdef object offset_lengths
|
||||
|
||||
+40
-42
@@ -77,47 +77,43 @@ def _offset_edge_map(shape, offsets):
|
||||
"""Return an array with positions marked where offsets will step
|
||||
out of bounds.
|
||||
|
||||
Given a shape (of length n) and a list of n-d offsets, return a shape + (n,)
|
||||
sized edge_map, where, for each dimension edge_map[...,dim] has zeros at
|
||||
indices at which none of the given offsets (in that dimension) will step
|
||||
out of bounds. If the value is nonzero, it gives the largest offset (in
|
||||
terms of absolute value) that will step out of bounds in that direction.
|
||||
Given a shape (of length n) and a list of n-d offsets, return a two arrays
|
||||
of (n,) + shape: pos_edge_map and neg_edge_map.
|
||||
For each dimension xxx_edge_map[dim, ...] has zeros at indices at which
|
||||
none of the given offsets (in that dimension) of the given sign (positive
|
||||
or negative, respectively) will step out of bounds. If the value is
|
||||
nonzero, it gives the largest offset (in terms of absolute value) that
|
||||
will step out of bounds in that direction.
|
||||
|
||||
An example will be explanatory:
|
||||
>>> offsets = [[-2,0], [1,1], [0,2]]
|
||||
>>> edge_map = _offset_edge_map((4,4), offsets)
|
||||
>>> edge_map[...,0]
|
||||
array([[-2, -2, -2, -2],
|
||||
>>> pos_edge_map, neg_edge_map = _offset_edge_map((4,4), offsets)
|
||||
>>> neg_edge_map[0]
|
||||
array([[-1, -1, -1, -1],
|
||||
[-2, -2, -2, -2],
|
||||
[ 0, 0, 0, 0],
|
||||
[ 1, 1, 1, 1]], dtype=int8)
|
||||
[ 0, 0, 0, 0]], dtype=int8)
|
||||
|
||||
>>> edge_map[...,1]
|
||||
>>> pos_edge_map[1]
|
||||
array([[0, 0, 2, 1],
|
||||
[0, 0, 2, 1],
|
||||
[0, 0, 2, 1],
|
||||
[0, 0, 2, 1]], dtype=int8)
|
||||
|
||||
"""
|
||||
d = len(shape)
|
||||
edges = np.zeros(shape+(d,), order='F', dtype=EDGE_D)
|
||||
indices = np.indices(shape) # indices.shape = (n,)+shape
|
||||
|
||||
#get the distance from each index to the upper or lower edge in each dim
|
||||
pos_edges = (shape - indices.T).T
|
||||
neg_edges = -1 - indices
|
||||
# now set the distances to zero if none of the given offsets could reach
|
||||
offsets = np.asarray(offsets)
|
||||
for i in range(d):
|
||||
slices = [slice(None)] * (d+1)
|
||||
slices[d] = i
|
||||
distinct_offsets = set(offsets[:,i])
|
||||
if 0 in distinct_offsets:
|
||||
distinct_offsets.remove(0)
|
||||
for offset in sorted(distinct_offsets, key=np.absolute, reverse=True):
|
||||
# process offsets with larger absolute values first, so that smaller
|
||||
# offsets will overwrite the correct region of the offsets array.
|
||||
slice_stop = -offset
|
||||
if offset > 0:
|
||||
slice_stop -= 1
|
||||
slice_step = -np.sign(offset)
|
||||
slices[i] = slice(None, slice_stop, slice_step)
|
||||
edges[tuple(slices)] = offset
|
||||
return edges
|
||||
maxes = offsets.max(axis=0)
|
||||
mins = offsets.min(axis=0)
|
||||
for pos, neg, mx, mn in zip(pos_edges, neg_edges, maxes, mins):
|
||||
pos[pos > mx] = 0
|
||||
neg[neg < mn] = 0
|
||||
return pos_edges.astype(EDGE_D), neg_edges.astype(EDGE_D)
|
||||
|
||||
def make_offsets(d, fully_connected):
|
||||
"""Make a list of offsets from a center point defining a n-dim
|
||||
@@ -296,9 +292,10 @@ cdef class MCP:
|
||||
# The edge map stores more than a boolean "on some edge" flag so as to
|
||||
# allow us to examine the non-out-of-bounds neighbors for a given edge
|
||||
# point while excluding the neighbors which are outside the array.
|
||||
self.flat_edge_map = \
|
||||
_offset_edge_map(costs.shape, self.offsets).reshape(
|
||||
(size, self.dim), order='F')
|
||||
pos, neg = _offset_edge_map(costs.shape, self.offsets)
|
||||
self.flat_pos_edge_map = pos.reshape((self.dim, size), order='F')
|
||||
self.flat_neg_edge_map = neg.reshape((self.dim, size), order='F')
|
||||
|
||||
|
||||
# The offset lengths are the distances traveled along each offset
|
||||
self.offset_lengths = np.sqrt(
|
||||
@@ -393,7 +390,10 @@ cdef class MCP:
|
||||
self.flat_cumulative_costs
|
||||
cdef np.ndarray[OFFSETS_INDEX_T, ndim=1] traceback_offsets = \
|
||||
self.traceback_offsets
|
||||
cdef np.ndarray[EDGE_T, ndim=2] flat_edge_map = self.flat_edge_map
|
||||
cdef np.ndarray[EDGE_T, ndim=2] flat_pos_edge_map = \
|
||||
self.flat_pos_edge_map
|
||||
cdef np.ndarray[EDGE_T, ndim=2] flat_neg_edge_map = \
|
||||
self.flat_neg_edge_map
|
||||
cdef np.ndarray[OFFSET_T, ndim=2] offsets = self.offsets
|
||||
cdef np.ndarray[INDEX_T, ndim=1] flat_offsets = self.flat_offsets
|
||||
cdef np.ndarray[FLOAT_T, ndim=1] offset_lengths = self.offset_lengths
|
||||
@@ -413,7 +413,7 @@ cdef class MCP:
|
||||
cdef BOOL_T is_at_edge, use_offset
|
||||
cdef INDEX_T d, i
|
||||
cdef OFFSET_T offset
|
||||
cdef EDGE_T edge_val
|
||||
cdef EDGE_T pos_edge_val, neg_edge_val
|
||||
cdef int num_ends_found = 0
|
||||
cdef FLOAT_T inf = np.inf
|
||||
cdef FLOAT_T travel_cost
|
||||
@@ -449,7 +449,8 @@ cdef class MCP:
|
||||
# edge along any axis
|
||||
is_at_edge = 0
|
||||
for d in range(dim):
|
||||
if flat_edge_map[index, d] != 0:
|
||||
if (flat_pos_edge_map[d, index] != 0 or
|
||||
flat_neg_edge_map[d, index] != 0):
|
||||
is_at_edge = 1
|
||||
break
|
||||
|
||||
@@ -466,14 +467,11 @@ cdef class MCP:
|
||||
if is_at_edge:
|
||||
for d in range(dim):
|
||||
offset = offsets[i, d]
|
||||
edge_val = flat_edge_map[index, d]
|
||||
if (offset < 0 and
|
||||
edge_val < 0 and
|
||||
offset <= edge_val) or \
|
||||
(offset > 0 and
|
||||
edge_val > 0 and
|
||||
offset >= edge_val):
|
||||
|
||||
pos_edge_val = flat_pos_edge_map[d, index]
|
||||
neg_edge_val = flat_neg_edge_map[d, index]
|
||||
if (pos_edge_val > 0 and offset >= pos_edge_val) or \
|
||||
(neg_edge_val < 0 and offset <= neg_edge_val):
|
||||
# the offset puts us out of bounds...
|
||||
use_offset = 0
|
||||
break
|
||||
# If not at an edge, or the specific offset doesn't
|
||||
|
||||
@@ -114,8 +114,23 @@ def test_no_diagonal():
|
||||
(7, 2)])
|
||||
|
||||
|
||||
def test_offsets():
|
||||
offsets = [(1,i) for i in range(10)] + [(1, -i) for i in range(1,10)]
|
||||
m = mcp.MCP(a, offsets=offsets)
|
||||
costs, traceback = m.find_costs([(1,6)])
|
||||
assert_array_equal(traceback,
|
||||
[[-1, -1, -1, -1, -1, -1, -1, -1],
|
||||
[-1, -1, -1, -1, -1, -1, -1, -1],
|
||||
[15, 14, 13, 12, 11, 10, 0, 1],
|
||||
[10, 0, 1, 2, 3, 4, 5, 6],
|
||||
[10, 0, 1, 2, 3, 4, 5, 6],
|
||||
[10, 0, 1, 2, 3, 4, 5, 6],
|
||||
[10, 0, 1, 2, 3, 4, 5, 6],
|
||||
[10, 0, 1, 2, 3, 4, 5, 6]])
|
||||
|
||||
|
||||
def test_crashing():
|
||||
for shape in [(100, 100), (5, 8, 13, 17)]:
|
||||
for shape in [(100, 100), (5, 8, 13, 17)]*5:
|
||||
yield _test_random, shape
|
||||
|
||||
def _test_random(shape):
|
||||
|
||||
+31
-16
@@ -8,22 +8,7 @@ from ._plugins import use as use_plugin
|
||||
from ._plugins import available as plugins
|
||||
from ._plugins import info as plugin_info
|
||||
from ._plugins import configuration as plugin_order
|
||||
available_plugins = plugins()
|
||||
|
||||
for preferred_plugin in ['matplotlib', 'pil', 'qt', 'freeimage', 'null']:
|
||||
if preferred_plugin in available_plugins:
|
||||
try:
|
||||
use_plugin(preferred_plugin)
|
||||
break
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
# Use PIL as the default imread plugin, since matplotlib (1.2.x)
|
||||
# is buggy (flips PNGs around, returns bytes as floats, etc.)
|
||||
try:
|
||||
use_plugin('pil', 'imread')
|
||||
except ImportError:
|
||||
pass
|
||||
from ._plugins import reset_plugins as _reset_plugins
|
||||
|
||||
from .sift import *
|
||||
from .collection import *
|
||||
@@ -32,6 +17,34 @@ from ._io import *
|
||||
from .video import *
|
||||
|
||||
|
||||
available_plugins = plugins()
|
||||
|
||||
|
||||
def _load_preferred_plugins():
|
||||
# Load preferred plugin for each io function.
|
||||
io_funcs = ['imsave', 'imshow', 'imread_collection', 'imread']
|
||||
preferred_plugins = ['matplotlib', 'pil', 'qt', 'freeimage', 'null']
|
||||
for func in io_funcs:
|
||||
for plugin in preferred_plugins:
|
||||
if plugin not in available_plugins:
|
||||
continue
|
||||
try:
|
||||
use_plugin(plugin, kind=func)
|
||||
break
|
||||
except (ImportError, RuntimeError):
|
||||
pass
|
||||
|
||||
# Use PIL as the default imread plugin, since matplotlib (1.2.x)
|
||||
# is buggy (flips PNGs around, returns bytes as floats, etc.)
|
||||
try:
|
||||
use_plugin('pil', 'imread')
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
def reset_plugins():
|
||||
_reset_plugins()
|
||||
_load_preferred_plugins()
|
||||
|
||||
def _update_doc(doc):
|
||||
"""Add a list of plugins to the module docstring, formatted as
|
||||
a ReStructuredText table.
|
||||
@@ -61,3 +74,5 @@ def _update_doc(doc):
|
||||
return doc
|
||||
|
||||
__doc__ = _update_doc(__doc__)
|
||||
|
||||
reset_plugins()
|
||||
|
||||
@@ -95,7 +95,11 @@ def imread_collection(load_pattern, conserve_memory=True):
|
||||
isinstance(hdu, pyfits.PrimaryHDU):
|
||||
# Ignore (primary) header units with no data (use '.size'
|
||||
# rather than '.data' to avoid actually loading the image):
|
||||
if hdu.size() > 0:
|
||||
try:
|
||||
data_size = hdu.size()
|
||||
except TypeError: # (size changed to int in PyFITS 3.1)
|
||||
data_size = hdu.size
|
||||
if data_size > 0:
|
||||
ext_list.append((filename, n))
|
||||
hdulist.close()
|
||||
|
||||
|
||||
@@ -55,19 +55,22 @@ def load_freeimage():
|
||||
try:
|
||||
freeimage = loader.LoadLibrary(lib)
|
||||
break
|
||||
except Exception, e:
|
||||
except Exception:
|
||||
if lib not in bare_libs:
|
||||
# Don't record errors when it couldn't load the library from
|
||||
# a bare name -- this fails often, and doesn't provide any
|
||||
# useful debugging information anyway, beyond "couldn't find
|
||||
# library..."
|
||||
errors.append((lib, e))
|
||||
# Get exception instance in Python 2.x/3.x compatible manner
|
||||
e_type, e_value, e_tb = sys.exc_info()
|
||||
del e_tb
|
||||
errors.append((lib, e_value))
|
||||
|
||||
if freeimage is None:
|
||||
if errors:
|
||||
# No freeimage library loaded, and load-errors reported for some
|
||||
# candidate libs
|
||||
err_txt = ['%s:\n%s'%(l, e.message) for l, e in errors]
|
||||
err_txt = ['%s:\n%s'%(l, str(e.message)) for l, e in errors]
|
||||
raise OSError('One or more FreeImage libraries were found, but '
|
||||
'could not be loaded due to the following errors:\n'+
|
||||
'\n\n'.join(err_txt))
|
||||
@@ -340,6 +343,9 @@ class METADATA_DATATYPE(object):
|
||||
FIDT_DOUBLE = 12 # 64-bit IEEE floating point
|
||||
FIDT_IFD = 13 # 32-bit unsigned integer (offset)
|
||||
FIDT_PALETTE = 14 # 32-bit RGBQUAD
|
||||
FIDT_LONG8 = 16 # 64-bit unsigned integer
|
||||
FIDT_SLONG8 = 17 # 64-bit signed integer
|
||||
FIDT_IFD8 = 18 # 64-bit unsigned integer (offset)
|
||||
|
||||
dtypes = {
|
||||
FIDT_BYTE: numpy.uint8,
|
||||
@@ -357,7 +363,10 @@ class METADATA_DATATYPE(object):
|
||||
FIDT_DOUBLE: numpy.float64,
|
||||
FIDT_IFD: numpy.uint32,
|
||||
FIDT_PALETTE: [('R', numpy.uint8), ('G', numpy.uint8),
|
||||
('B', numpy.uint8), ('A', numpy.uint8)]
|
||||
('B', numpy.uint8), ('A', numpy.uint8)],
|
||||
FIDT_LONG8: numpy.uint64,
|
||||
FIDT_SLONG8: numpy.int64,
|
||||
FIDT_IFD8: numpy.uint64
|
||||
}
|
||||
|
||||
|
||||
|
||||
@@ -2,23 +2,32 @@
|
||||
|
||||
"""
|
||||
|
||||
__all__ = ['use', 'available', 'call', 'info', 'configuration']
|
||||
__all__ = ['use', 'available', 'call', 'info', 'configuration', 'reset_plugins']
|
||||
|
||||
import warnings
|
||||
from ConfigParser import ConfigParser
|
||||
import os.path
|
||||
from glob import glob
|
||||
|
||||
plugin_store = {'imread': [],
|
||||
'imsave': [],
|
||||
'imshow': [],
|
||||
'imread_collection': [],
|
||||
'_app_show': []}
|
||||
plugin_store = None
|
||||
|
||||
plugin_provides = {}
|
||||
plugin_module_name = {}
|
||||
plugin_meta_data = {}
|
||||
|
||||
def reset_plugins():
|
||||
"""Clear the plugin state to the default, i.e., where no plugins are loaded.
|
||||
|
||||
"""
|
||||
global plugin_store
|
||||
plugin_store = {'imread': [],
|
||||
'imsave': [],
|
||||
'imshow': [],
|
||||
'imread_collection': [],
|
||||
'_app_show': []}
|
||||
|
||||
reset_plugins()
|
||||
|
||||
def _scan_plugins():
|
||||
"""Scan the plugins directory for .ini files and parse them
|
||||
to gather plugin meta-data.
|
||||
|
||||
@@ -3,4 +3,4 @@ try:
|
||||
except ImportError:
|
||||
raise ImportError("The tifffile module could not be found.\n"
|
||||
"It can be obtained at "
|
||||
"<http://www.lfd.uci.edu/~gohlke/code/tifffile.py.html>\n")
|
||||
"<http://www.lfd.uci.edu/~gohlke/code/tifffile.py>\n")
|
||||
|
||||
@@ -9,7 +9,7 @@ from skimage.util import img_as_ubyte
|
||||
try:
|
||||
import multiprocessing
|
||||
CPU_COUNT = multiprocessing.cpu_count()
|
||||
except ImportError:
|
||||
except:
|
||||
CPU_COUNT = 2
|
||||
|
||||
class GuiLockError(Exception):
|
||||
|
||||
@@ -26,6 +26,11 @@ def test_fits_plugin_import():
|
||||
else:
|
||||
assert pyfits_available == True
|
||||
|
||||
|
||||
def teardown():
|
||||
io.reset_plugins()
|
||||
|
||||
|
||||
@skipif(not pyfits_available)
|
||||
def test_imread_MEF():
|
||||
io.use_plugin('fits')
|
||||
|
||||
@@ -27,6 +27,10 @@ def setup_module(self):
|
||||
pass
|
||||
|
||||
|
||||
def teardown():
|
||||
sio.reset_plugins()
|
||||
|
||||
|
||||
@skipif(not FI_available)
|
||||
def test_imread():
|
||||
img = sio.imread(os.path.join(si.data_dir, 'color.png'))
|
||||
|
||||
@@ -6,7 +6,7 @@ from numpy.testing.decorators import skipif
|
||||
from tempfile import NamedTemporaryFile
|
||||
|
||||
from skimage import data_dir
|
||||
from skimage.io import imread, imsave, use_plugin
|
||||
from skimage.io import imread, imsave, use_plugin, reset_plugins
|
||||
|
||||
try:
|
||||
from PIL import Image
|
||||
@@ -18,6 +18,10 @@ else:
|
||||
PIL_available = True
|
||||
|
||||
|
||||
def teardown():
|
||||
reset_plugins()
|
||||
|
||||
|
||||
def setup_module(self):
|
||||
"""The effect of the `plugin.use` call may be overridden by later imports.
|
||||
Call `use_plugin` directly before the tests to ensure that PIL is used.
|
||||
@@ -65,10 +69,12 @@ def test_bilevel():
|
||||
def test_imread_uint16():
|
||||
expected = np.load(os.path.join(data_dir, 'chessboard_GRAY_U8.npy'))
|
||||
img = imread(os.path.join(data_dir, 'chessboard_GRAY_U16.tif'))
|
||||
assert img.dtype == np.uint16
|
||||
assert np.issubdtype(img.dtype, np.uint16)
|
||||
assert_array_almost_equal(img, expected)
|
||||
|
||||
@skipif(not PIL_available)
|
||||
# Big endian images not correctly loaded for PIL < 1.1.7
|
||||
# Renable test when PIL 1.1.7 is more common.
|
||||
@skipif(True)
|
||||
def test_imread_uint16_big_endian():
|
||||
expected = np.load(os.path.join(data_dir, 'chessboard_GRAY_U8.npy'))
|
||||
img = imread(os.path.join(data_dir, 'chessboard_GRAY_U16B.tif'))
|
||||
@@ -97,3 +103,6 @@ class TestSave:
|
||||
else:
|
||||
x = (x * 255).astype(dtype)
|
||||
yield self.roundtrip, dtype, x
|
||||
|
||||
if __name__ == "__main__":
|
||||
run_module_suite()
|
||||
|
||||
@@ -4,8 +4,6 @@ from skimage import io
|
||||
from skimage.io._plugins import plugin
|
||||
from numpy.testing.decorators import skipif
|
||||
|
||||
from copy import deepcopy
|
||||
|
||||
try:
|
||||
io.use_plugin('pil')
|
||||
PIL_available = True
|
||||
@@ -22,11 +20,10 @@ except OSError:
|
||||
|
||||
|
||||
def setup_module(self):
|
||||
self.backup_plugin_store = deepcopy(plugin.plugin_store)
|
||||
plugin.use('test') # see ../_plugins/test_plugin.py
|
||||
|
||||
def teardown_module(self):
|
||||
plugin.plugin_store = self.backup_plugin_store
|
||||
io.reset_plugins()
|
||||
|
||||
class TestPlugin:
|
||||
def test_read(self):
|
||||
|
||||
@@ -17,9 +17,7 @@ except ImportError:
|
||||
|
||||
|
||||
def teardown():
|
||||
if TF_available:
|
||||
for k, v in _plugins.items():
|
||||
sio.use_plugin(v[0], k)
|
||||
sio.reset_plugins()
|
||||
|
||||
|
||||
@skipif(not TF_available)
|
||||
|
||||
@@ -1 +1,3 @@
|
||||
from .find_contours import find_contours
|
||||
from ._regionprops import regionprops
|
||||
from ._structural_similarity import structural_similarity
|
||||
|
||||
@@ -0,0 +1,52 @@
|
||||
#cython: boundscheck=False
|
||||
#cython: wraparound=False
|
||||
#cython: cdivision=True
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
|
||||
|
||||
def central_moments(np.ndarray[np.double_t, ndim=2] array, double cr, double cc,
|
||||
int order):
|
||||
cdef int p, q, r, c
|
||||
cdef np.ndarray[np.double_t, ndim=2] mu
|
||||
mu = np.zeros((order + 1, order + 1), 'double')
|
||||
for p in range(order + 1):
|
||||
for q in range(order + 1):
|
||||
for r in range(array.shape[0]):
|
||||
for c in range(array.shape[1]):
|
||||
mu[p,q] += array[r,c] * (r - cr) ** q * (c - cc) ** p
|
||||
return mu
|
||||
|
||||
def normalized_moments(np.ndarray[np.double_t, ndim=2] mu, int order):
|
||||
cdef int p, q
|
||||
cdef np.ndarray[np.double_t, ndim=2] nu
|
||||
nu = np.zeros((order + 1, order + 1), 'double')
|
||||
for p in range(order + 1):
|
||||
for q in range(order + 1):
|
||||
if p + q >= 2:
|
||||
nu[p,q] = mu[p,q] / mu[0,0]**(<double>(p + q) / 2 + 1)
|
||||
else:
|
||||
nu[p,q] = np.nan
|
||||
return nu
|
||||
|
||||
def hu_moments(np.ndarray[np.double_t, ndim=2] nu):
|
||||
cdef np.ndarray[np.double_t, ndim=1] hu = np.zeros((7,), 'double')
|
||||
cdef double t0 = nu[3,0] + nu[1,2]
|
||||
cdef double t1 = nu[2,1] + nu[0,3]
|
||||
cdef double q0 = t0 * t0
|
||||
cdef double q1 = t1 * t1
|
||||
cdef double n4 = 4 * nu[1,1]
|
||||
cdef double s = nu[2,0] + nu[0,2]
|
||||
cdef double d = nu[2,0] - nu[0,2]
|
||||
hu[0] = s
|
||||
hu[1] = d * d + n4 * nu[1,1]
|
||||
hu[3] = q0 + q1
|
||||
hu[5] = d * (q0 - q1) + n4 * t0 * t1
|
||||
t0 *= q0 - 3 * q1
|
||||
t1 *= 3 * q0 - q1
|
||||
q0 = nu[3,0]- 3 * nu[1,2]
|
||||
q1 = 3 * nu[2,1] - nu[0,3]
|
||||
hu[2] = q0 * q0 + q1 * q1
|
||||
hu[4] = q0 * t0 + q1 * t1
|
||||
hu[6] = q1 * t0 - q0 * t1
|
||||
return hu
|
||||
@@ -0,0 +1,353 @@
|
||||
# coding: utf-8
|
||||
from math import sqrt, atan, pi as PI
|
||||
import numpy as np
|
||||
from scipy import ndimage
|
||||
|
||||
from skimage.morphology import convex_hull_image
|
||||
from . import _moments
|
||||
|
||||
|
||||
__all__ = ['regionprops']
|
||||
|
||||
|
||||
STREL_8 = np.ones((3, 3), 'int8')
|
||||
PROPS = (
|
||||
'Area',
|
||||
'BoundingBox',
|
||||
'CentralMoments',
|
||||
'Centroid',
|
||||
'ConvexArea',
|
||||
# 'ConvexHull',
|
||||
'ConvexImage',
|
||||
'Eccentricity',
|
||||
'EquivDiameter',
|
||||
'EulerNumber',
|
||||
'Extent',
|
||||
# 'Extrema',
|
||||
'FilledArea',
|
||||
'FilledImage',
|
||||
'HuMoments',
|
||||
'Image',
|
||||
'MajorAxisLength',
|
||||
'MaxIntensity',
|
||||
'MeanIntensity',
|
||||
'MinIntensity',
|
||||
'MinorAxisLength',
|
||||
'Moments',
|
||||
'NormalizedMoments',
|
||||
'Orientation',
|
||||
# 'Perimeter',
|
||||
# 'PixelIdxList',
|
||||
# 'PixelList',
|
||||
'Solidity',
|
||||
# 'SubarrayIdx'
|
||||
'WeightedCentralMoments',
|
||||
'WeightedCentroid',
|
||||
'WeightedHuMoments',
|
||||
'WeightedMoments',
|
||||
'WeightedNormalizedMoments'
|
||||
)
|
||||
|
||||
|
||||
def regionprops(label_image, properties=['Area', 'Centroid'],
|
||||
intensity_image=None):
|
||||
"""Measure properties of labelled image regions.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
label_image : N x M ndarray
|
||||
Labelled input image.
|
||||
properties : {'all', list}
|
||||
Shape measurements to be determined for each labelled image region.
|
||||
Default is `['Area', 'Centroid']`. The following properties can be
|
||||
determined:
|
||||
* Area : int
|
||||
Number of pixels of region.
|
||||
* BoundingBox : tuple
|
||||
Bounding box `(min_row, min_col, max_row, max_col)`
|
||||
* CentralMoments : 3 x 3 ndarray
|
||||
Central moments (translation invariant) up to 3rd order.
|
||||
mu_ji = sum{ array(x, y) * (x - x_c)^j * (y - y_c)^i }
|
||||
where the sum is over the `x`, `y` coordinates of the region,
|
||||
and `x_c` and `y_c` are the coordinates of the region's centroid.
|
||||
* Centroid : array
|
||||
Centroid coordinate tuple `(row, col)`.
|
||||
* ConvexArea : int
|
||||
Number of pixels of convex hull image.
|
||||
* ConvexImage : H x J ndarray
|
||||
Binary convex hull image which has the same size as bounding box.
|
||||
* Eccentricity : float
|
||||
Eccentricity of the ellipse that has the same second-moments as the
|
||||
region. The eccentricity is the ratio of the distance between its
|
||||
minor and major axis length. The value is between 0 and 1.
|
||||
* EquivDiameter : float
|
||||
The diameter of a circle with the same area as the region.
|
||||
* EulerNumber : int
|
||||
Euler number of region. Computed as number of objects (= 1)
|
||||
subtracted by number of holes (8-connectivity).
|
||||
* Extent : float
|
||||
Ratio of pixels in the region to pixels in the total bounding box.
|
||||
Computed as `Area / (rows*cols)`
|
||||
* FilledArea : int
|
||||
Number of pixels of filled region.
|
||||
* FilledImage : H x J ndarray
|
||||
Binary region image with filled holes which has the same size as
|
||||
bounding box.
|
||||
* HuMoments : tuple
|
||||
Hu moments (translation, scale and rotation invariant).
|
||||
* Image : H x J ndarray
|
||||
Sliced binary region image which has the same size as bounding box.
|
||||
* MajorAxisLength : float
|
||||
The length of the major axis of the ellipse that has the same
|
||||
normalized second central moments as the region.
|
||||
* MaxIntensity: float
|
||||
Value with the greatest intensity in the region.
|
||||
* MeanIntensity: float
|
||||
Value with the mean intensity in the region.
|
||||
* MinIntensity: float
|
||||
Value with the least intensity in the region.
|
||||
* MinorAxisLength : float
|
||||
The length of the minor axis of the ellipse that has the same
|
||||
normalized second central moments as the region.
|
||||
* Moments : 3 x 3 ndarray
|
||||
Spatial moments up to 3rd order.
|
||||
m_ji = sum{ array(x, y) * x^j * y^i }
|
||||
where the sum is over the `x`, `y` coordinates of the region.
|
||||
* NormalizedMoments : 3 x 3 ndarray
|
||||
Normalized moments (translation and scale invariant) up to 3rd
|
||||
order.
|
||||
nu_ji = mu_ji / m_00^[(i+j)/2 + 1]
|
||||
where `m_00` is the zeroth spatial moment.
|
||||
* Orientation : float
|
||||
Angle between the X-axis and the major axis of the ellipse that has
|
||||
the same second-moments as the region. Ranging from `-pi/2` to
|
||||
`-pi/2` in counter-clockwise direction.
|
||||
* Solidity : float
|
||||
Ratio of pixels in the region to pixels of the convex hull image.
|
||||
* WeightedCentralMoments : 3 x 3 ndarray
|
||||
Central moments (translation invariant) of intensity image up to 3rd
|
||||
order.
|
||||
wmu_ji = sum{ array(x, y) * (x - x_c)^j * (y - y_c)^i }
|
||||
where the sum is over the `x`, `y` coordinates of the region,
|
||||
and `x_c` and `y_c` are the coordinates of the region's centroid.
|
||||
* WeightedCentroid : array
|
||||
Centroid coordinate tuple `(row, col)` weighted with intensity
|
||||
image.
|
||||
* WeightedHuMoments : tuple
|
||||
Hu moments (translation, scale and rotation invariant) of intensity
|
||||
image.
|
||||
* WeightedMoments : 3 x 3 ndarray
|
||||
Spatial moments of intensity image up to 3rd order.
|
||||
wm_ji = sum{ array(x, y) * x^j * y^i }
|
||||
where the sum is over the `x`, `y` coordinates of the region.
|
||||
* WeightedNormalizedMoments : 3 x 3 ndarray
|
||||
Normalized moments (translation and scale invariant) of intensity
|
||||
image up to 3rd order.
|
||||
wnu_ji = wmu_ji / wm_00^[(i+j)/2 + 1]
|
||||
where `wm_00` is the zeroth spatial moment (intensity-weighted
|
||||
area).
|
||||
intensity_image : N x M ndarray, optional
|
||||
Intensity image with same size as labelled image. Default is None.
|
||||
|
||||
Returns
|
||||
-------
|
||||
properties : list of dicts
|
||||
List containing a property dict for each region. The property dicts
|
||||
contain all the specified properties plus a 'Label' field.
|
||||
|
||||
References
|
||||
----------
|
||||
Wilhelm Burger, Mark Burge. Principles of Digital Image Processing: Core
|
||||
Algorithms. Springer-Verlag, London, 2009.
|
||||
B. Jähne. Digital Image Processing. Springer-Verlag,
|
||||
Berlin-Heidelberg, 6. edition, 2005.
|
||||
T. H. Reiss. Recognizing Planar Objects Using Invariant Image Features,
|
||||
from Lecture notes in computer science, p. 676. Springer, Berlin, 1993.
|
||||
http://en.wikipedia.org/wiki/Image_moment
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from skimage.data import coins
|
||||
>>> from skimage.morphology import label
|
||||
>>> img = coins() > 110
|
||||
>>> label_img = label(img)
|
||||
>>> props = regionprops(label_img)
|
||||
>>> props[0]['Centroid'] # centroid of first labelled object
|
||||
"""
|
||||
if not np.issubdtype(label_image.dtype, 'int'):
|
||||
raise TypeError('labelled image must be of integer dtype')
|
||||
|
||||
# determine all properties if nothing specified
|
||||
if properties == 'all':
|
||||
properties = PROPS
|
||||
|
||||
props = []
|
||||
|
||||
objects = ndimage.find_objects(label_image)
|
||||
for i, sl in enumerate(objects):
|
||||
label = i + 1
|
||||
|
||||
# create property dict for current label
|
||||
obj_props = {}
|
||||
props.append(obj_props)
|
||||
|
||||
obj_props['Label'] = label
|
||||
|
||||
array = (label_image[sl] == label).astype('double')
|
||||
|
||||
# upper left corner of object bbox
|
||||
r0 = sl[0].start
|
||||
c0 = sl[1].start
|
||||
|
||||
m = _moments.central_moments(array, 0, 0, 3)
|
||||
# centroid
|
||||
cr = m[0,1] / m[0,0]
|
||||
cc = m[1,0] / m[0,0]
|
||||
mu = _moments.central_moments(array, cr, cc, 3)
|
||||
|
||||
#: elements of the inertia tensor [a b; b c]
|
||||
a = mu[2,0] / mu[0,0]
|
||||
b = mu[1,1] / mu[0,0]
|
||||
c = mu[0,2] / mu[0,0]
|
||||
#: eigen values of inertia tensor
|
||||
l1 = (a + c) / 2 + sqrt(4 * b ** 2 + (a - c) ** 2) / 2
|
||||
l2 = (a + c) / 2 - sqrt(4 * b ** 2 + (a - c) ** 2) / 2
|
||||
|
||||
# cached results which are used by several properties
|
||||
_filled_image = None
|
||||
_convex_image = None
|
||||
_nu = None
|
||||
|
||||
if 'Area' in properties:
|
||||
obj_props['Area'] = m[0,0]
|
||||
|
||||
if 'BoundingBox' in properties:
|
||||
obj_props['BoundingBox'] = (r0, c0, sl[0].stop, sl[1].stop)
|
||||
|
||||
if 'Centroid' in properties:
|
||||
obj_props['Centroid'] = cr + r0, cc + c0
|
||||
|
||||
if 'CentralMoments' in properties:
|
||||
obj_props['CentralMoments'] = mu
|
||||
|
||||
if 'ConvexArea' in properties:
|
||||
if _convex_image is None:
|
||||
_convex_image = convex_hull_image(array)
|
||||
obj_props['ConvexArea'] = np.sum(_convex_image)
|
||||
|
||||
if 'ConvexImage' in properties:
|
||||
if _convex_image is None:
|
||||
_convex_image = convex_hull_image(array)
|
||||
obj_props['ConvexImage'] = _convex_image
|
||||
|
||||
if 'Eccentricity' in properties:
|
||||
if l1 == 0:
|
||||
obj_props['Eccentricity'] = 0
|
||||
else:
|
||||
obj_props['Eccentricity'] = sqrt(1 - l2 / l1)
|
||||
|
||||
if 'EquivDiameter' in properties:
|
||||
obj_props['EquivDiameter'] = sqrt(4 * m[0,0] / PI)
|
||||
|
||||
if 'EulerNumber' in properties:
|
||||
if _filled_image is None:
|
||||
_filled_image = ndimage.binary_fill_holes(array, STREL_8)
|
||||
euler_array = _filled_image != array
|
||||
_, num = ndimage.label(euler_array, STREL_8)
|
||||
obj_props['EulerNumber'] = - num
|
||||
|
||||
if 'Extent' in properties:
|
||||
obj_props['Extent'] = m[0,0] / (array.shape[0] * array.shape[1])
|
||||
|
||||
if 'HuMoments' in properties:
|
||||
if _nu is None:
|
||||
_nu = _moments.normalized_moments(mu, 3)
|
||||
obj_props['HuMoments'] = _moments.hu_moments(_nu)
|
||||
|
||||
if 'Image' in properties:
|
||||
obj_props['Image'] = array
|
||||
|
||||
if 'FilledArea' in properties:
|
||||
if _filled_image is None:
|
||||
_filled_image = ndimage.binary_fill_holes(array, STREL_8)
|
||||
obj_props['FilledArea'] = np.sum(_filled_image)
|
||||
|
||||
if 'FilledImage' in properties:
|
||||
if _filled_image is None:
|
||||
_filled_image = ndimage.binary_fill_holes(array, STREL_8)
|
||||
obj_props['FilledImage'] = _filled_image
|
||||
|
||||
if 'MajorAxisLength' in properties:
|
||||
obj_props['MajorAxisLength'] = 4 * sqrt(l1)
|
||||
|
||||
if 'MinorAxisLength' in properties:
|
||||
obj_props['MinorAxisLength'] = 4 * sqrt(l2)
|
||||
|
||||
if 'Moments' in properties:
|
||||
obj_props['Moments'] = m
|
||||
|
||||
if 'NormalizedMoments' in properties:
|
||||
if _nu is None:
|
||||
_nu = _moments.normalized_moments(mu, 3)
|
||||
obj_props['NormalizedMoments'] = _nu
|
||||
|
||||
if 'Orientation' in properties:
|
||||
if a - c == 0:
|
||||
obj_props['Orientation'] = PI / 2
|
||||
else:
|
||||
obj_props['Orientation'] = - 0.5 * atan(2 * b / (a - c))
|
||||
|
||||
if 'Solidity' in properties:
|
||||
if _convex_image is None:
|
||||
_convex_image = convex_hull_image(array)
|
||||
obj_props['Solidity'] = m[0,0] / np.sum(_convex_image)
|
||||
|
||||
|
||||
if intensity_image is not None:
|
||||
weighted_array = array * intensity_image[sl]
|
||||
|
||||
wm = _moments.central_moments(weighted_array, 0, 0, 3)
|
||||
# weighted centroid
|
||||
wcr = wm[0,1] / wm[0,0]
|
||||
wcc = wm[1,0] / wm[0,0]
|
||||
wmu = _moments.central_moments(weighted_array, wcr, wcc, 3)
|
||||
|
||||
# cached results which are used by several properties
|
||||
_wnu = None
|
||||
_vals = None
|
||||
|
||||
if 'MaxIntensity' in properties:
|
||||
if _vals is None:
|
||||
_vals = weighted_array[array.astype('bool')]
|
||||
obj_props['MaxIntensity'] = np.max(_vals)
|
||||
|
||||
if 'MeanIntensity' in properties:
|
||||
if _vals is None:
|
||||
_vals = weighted_array[array.astype('bool')]
|
||||
obj_props['MeanIntensity'] = np.mean(_vals)
|
||||
|
||||
if 'MinIntensity' in properties:
|
||||
if _vals is None:
|
||||
_vals = weighted_array[array.astype('bool')]
|
||||
obj_props['MinIntensity'] = np.min(_vals)
|
||||
|
||||
if 'WeightedCentralMoments' in properties:
|
||||
obj_props['WeightedCentralMoments'] = wmu
|
||||
|
||||
if 'WeightedCentroid' in properties:
|
||||
obj_props['WeightedCentroid'] = wcr + r0, wcc + c0
|
||||
|
||||
if 'WeightedHuMoments' in properties:
|
||||
if _wnu is None:
|
||||
_wnu = _moments.normalized_moments(wmu, 3)
|
||||
obj_props['WeightedHuMoments'] = _moments.hu_moments(_wnu)
|
||||
|
||||
if 'WeightedMoments' in properties:
|
||||
obj_props['WeightedMoments'] = wm
|
||||
|
||||
if 'WeightedNormalizedMoments' in properties:
|
||||
if _wnu is None:
|
||||
_wnu = _moments.normalized_moments(wmu, 3)
|
||||
obj_props['WeightedNormalizedMoments'] = _wnu
|
||||
|
||||
return props
|
||||
@@ -0,0 +1,105 @@
|
||||
from __future__ import division
|
||||
|
||||
__all__ = ['structural_similarity']
|
||||
|
||||
import numpy as np
|
||||
|
||||
from ..util.dtype import dtype_range
|
||||
from ..util.shape import view_as_windows
|
||||
|
||||
def structural_similarity(X, Y, win_size=7,
|
||||
gradient=False, dynamic_range=None):
|
||||
"""Compute the mean structural similarity index between two images.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
X, Y : (N,N) ndarray
|
||||
Images.
|
||||
win_size : int
|
||||
The side-length of the sliding window used in comparison. Must
|
||||
be an odd value.
|
||||
gradient : bool
|
||||
If True, also return the gradient.
|
||||
dynamic_range : int
|
||||
Dynamic range of the input image (distance between minimum and
|
||||
maximum possible values). By default, this is estimated from
|
||||
the image data-type.
|
||||
|
||||
Returns
|
||||
-------
|
||||
s : float
|
||||
Strucutural similarity.
|
||||
grad : (N * N,) ndarray
|
||||
Gradient of the structural similarity index between X and Y.
|
||||
This is only returned if `gradient` is set to True.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P.
|
||||
(2004). Image quality assessment: From error visibility to
|
||||
structural similarity. IEEE Transactions on Image Processing,
|
||||
13, 600-612.
|
||||
|
||||
"""
|
||||
if not X.dtype == Y.dtype:
|
||||
raise ValueError('Input images must have the same dtype.')
|
||||
|
||||
if not X.shape == Y.shape:
|
||||
raise ValueError('Input images must have the same dimensions.')
|
||||
|
||||
if not (win_size % 2 == 1):
|
||||
raise ValueError('Window size must be odd.')
|
||||
|
||||
if dynamic_range is None:
|
||||
dmin, dmax = dtype_range[X.dtype.type]
|
||||
dynamic_range = dmax - dmin
|
||||
|
||||
XW = view_as_windows(X, (win_size, win_size))
|
||||
YW = view_as_windows(Y, (win_size, win_size))
|
||||
|
||||
NS = len(XW)
|
||||
NP = win_size * win_size
|
||||
|
||||
ux = np.mean(np.mean(XW, axis=2), axis=2)
|
||||
uy = np.mean(np.mean(YW, axis=2), axis=2)
|
||||
|
||||
# Compute variances var(X), var(Y) and var(X, Y)
|
||||
cov_norm = 1 / (win_size**2 - 1)
|
||||
XWM = XW - ux[..., None, None]
|
||||
YWM = YW - uy[..., None, None]
|
||||
vx = cov_norm * np.sum(np.sum(XWM**2, axis=2), axis=2)
|
||||
vy = cov_norm * np.sum(np.sum(YWM**2, axis=2), axis=2)
|
||||
vxy = cov_norm * np.sum(np.sum(XWM * YWM, axis=2), axis=2)
|
||||
|
||||
R = dynamic_range
|
||||
K1 = 0.01
|
||||
K2 = 0.03
|
||||
C1 = (K1 * R)**2
|
||||
C2 = (K2 * R)**2
|
||||
|
||||
A1, A2, B1, B2 = (v[..., None, None] for v in
|
||||
(2 * ux * uy + C1,
|
||||
2 * vxy + C2,
|
||||
ux**2 + uy**2 + C1,
|
||||
vx + vy + C2))
|
||||
|
||||
S = np.mean((A1 * A2) / (B1 * B2))
|
||||
|
||||
if gradient:
|
||||
local_grad = 2 / (NP * B1**2 * B2**2) * \
|
||||
(
|
||||
A1 * B1 * (B2 * XW - A2 * YW) - \
|
||||
B1 * B2 * (A2 - A1) * ux[..., None, None] + \
|
||||
A1 * A2 * (B1 - B2) * uy[..., None, None]
|
||||
)
|
||||
|
||||
grad = np.zeros_like(X, dtype=float)
|
||||
OW = view_as_windows(grad, (win_size, win_size))
|
||||
|
||||
OW += local_grad
|
||||
grad /= NS
|
||||
|
||||
return S, grad
|
||||
|
||||
else:
|
||||
return S
|
||||
@@ -12,9 +12,12 @@ def configuration(parent_package='', top_path=None):
|
||||
config.add_data_dir('tests')
|
||||
|
||||
cython(['_find_contours.pyx'], working_path=base_path)
|
||||
cython(['_moments.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('_find_contours', sources=['_find_contours.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_moments', sources=['_moments.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
|
||||
@@ -0,0 +1,249 @@
|
||||
from numpy.testing import assert_array_equal, assert_almost_equal, \
|
||||
assert_array_almost_equal
|
||||
import numpy as np
|
||||
|
||||
from skimage.measure import regionprops
|
||||
|
||||
|
||||
SAMPLE = np.array(
|
||||
[[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
|
||||
[0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
|
||||
[1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0],
|
||||
[0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1],
|
||||
[0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1]]
|
||||
)
|
||||
INTENSITY_SAMPLE = SAMPLE.copy()
|
||||
INTENSITY_SAMPLE[1,9:11] = 2
|
||||
|
||||
|
||||
def test_area():
|
||||
area = regionprops(SAMPLE, ['Area'])[0]['Area']
|
||||
assert area == np.sum(SAMPLE)
|
||||
|
||||
def test_bbox():
|
||||
bbox = regionprops(SAMPLE, ['BoundingBox'])[0]['BoundingBox']
|
||||
assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1]))
|
||||
|
||||
SAMPLE_mod = SAMPLE.copy()
|
||||
SAMPLE_mod[:,-1] = 0
|
||||
bbox = regionprops(SAMPLE_mod, ['BoundingBox'])[0]['BoundingBox']
|
||||
assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1]-1))
|
||||
|
||||
def test_central_moments():
|
||||
mu = regionprops(SAMPLE, ['CentralMoments'])[0]['CentralMoments']
|
||||
#: determined with OpenCV
|
||||
assert_almost_equal(mu[0,2], 436.00000000000045)
|
||||
# different from OpenCV results, bug in OpenCV
|
||||
assert_almost_equal(mu[0,3], -737.333333333333)
|
||||
assert_almost_equal(mu[1,1], -87.33333333333303)
|
||||
assert_almost_equal(mu[1,2], -127.5555555555593)
|
||||
assert_almost_equal(mu[2,0], 1259.7777777777774)
|
||||
assert_almost_equal(mu[2,1], 2000.296296296291)
|
||||
assert_almost_equal(mu[3,0], -760.0246913580195)
|
||||
|
||||
def test_centroid():
|
||||
centroid = regionprops(SAMPLE, ['Centroid'])[0]['Centroid']
|
||||
# determined with MATLAB
|
||||
assert_array_almost_equal(centroid, (5.66666666666666, 9.444444444444444))
|
||||
|
||||
def test_convex_area():
|
||||
area = regionprops(SAMPLE, ['ConvexArea'])[0]['ConvexArea']
|
||||
# determined with MATLAB
|
||||
assert area == 124
|
||||
|
||||
def test_convex_image():
|
||||
img = regionprops(SAMPLE, ['ConvexImage'])[0]['ConvexImage']
|
||||
# determined with MATLAB
|
||||
ref = np.array(
|
||||
[[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
|
||||
[0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
|
||||
[0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
|
||||
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
|
||||
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
|
||||
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
|
||||
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]]
|
||||
)
|
||||
assert_array_equal(img, ref)
|
||||
|
||||
def test_eccentricity():
|
||||
eps = regionprops(SAMPLE, ['Eccentricity'])[0]['Eccentricity']
|
||||
assert_almost_equal(eps, 0.814629313427)
|
||||
|
||||
def test_equiv_diameter():
|
||||
diameter = regionprops(SAMPLE, ['EquivDiameter'])[0]['EquivDiameter']
|
||||
# determined with MATLAB
|
||||
assert_almost_equal(diameter, 9.57461472963)
|
||||
|
||||
def test_euler_number():
|
||||
en = regionprops(SAMPLE, ['EulerNumber'])[0]['EulerNumber']
|
||||
assert en == 0
|
||||
|
||||
SAMPLE_mod = SAMPLE.copy()
|
||||
SAMPLE_mod[7,-3] = 0
|
||||
en = regionprops(SAMPLE_mod, ['EulerNumber'])[0]['EulerNumber']
|
||||
assert en == -1
|
||||
|
||||
def test_extent():
|
||||
extent = regionprops(SAMPLE, ['Extent'])[0]['Extent']
|
||||
assert_almost_equal(extent, 0.4)
|
||||
|
||||
def test_hu_moments():
|
||||
hu = regionprops(SAMPLE, ['HuMoments'])[0]['HuMoments']
|
||||
ref = np.array([
|
||||
3.27117627e-01,
|
||||
2.63869194e-02,
|
||||
2.35390060e-02,
|
||||
1.23151193e-03,
|
||||
1.38882330e-06,
|
||||
-2.72586158e-05,
|
||||
6.48350653e-06
|
||||
])
|
||||
# bug in OpenCV caused in Central Moments calculation?
|
||||
assert_array_almost_equal(hu, ref)
|
||||
|
||||
def test_image():
|
||||
img = regionprops(SAMPLE, ['Image'])[0]['Image']
|
||||
assert_array_equal(img, SAMPLE)
|
||||
|
||||
def test_filled_area():
|
||||
area = regionprops(SAMPLE, ['FilledArea'])[0]['FilledArea']
|
||||
assert area == np.sum(SAMPLE)
|
||||
|
||||
SAMPLE_mod = SAMPLE.copy()
|
||||
SAMPLE_mod[7,-3] = 0
|
||||
area = regionprops(SAMPLE_mod, ['FilledArea'])[0]['FilledArea']
|
||||
assert area == np.sum(SAMPLE)
|
||||
|
||||
def test_major_axis_length():
|
||||
length = regionprops(SAMPLE, ['MajorAxisLength'])[0]['MajorAxisLength']
|
||||
# MATLAB has different interpretation of ellipse than found in literature,
|
||||
# here implemented as found in literature
|
||||
assert_almost_equal(length, 16.7924234999)
|
||||
|
||||
def test_max_intensity():
|
||||
intensity = regionprops(SAMPLE, ['MaxIntensity'], INTENSITY_SAMPLE
|
||||
)[0]['MaxIntensity']
|
||||
assert_almost_equal(intensity, 2)
|
||||
|
||||
def test_mean_intensity():
|
||||
intensity = regionprops(SAMPLE, ['MeanIntensity'], INTENSITY_SAMPLE
|
||||
)[0]['MeanIntensity']
|
||||
assert_almost_equal(intensity, 1.02777777777777)
|
||||
|
||||
def test_min_intensity():
|
||||
intensity = regionprops(SAMPLE, ['MinIntensity'], INTENSITY_SAMPLE
|
||||
)[0]['MinIntensity']
|
||||
assert_almost_equal(intensity, 1)
|
||||
|
||||
def test_minor_axis_length():
|
||||
length = regionprops(SAMPLE, ['MinorAxisLength'])[0]['MinorAxisLength']
|
||||
# MATLAB has different interpretation of ellipse than found in literature,
|
||||
# here implemented as found in literature
|
||||
assert_almost_equal(length, 9.739302807263)
|
||||
|
||||
def test_moments():
|
||||
m = regionprops(SAMPLE, ['Moments'])[0]['Moments']
|
||||
#: determined with OpenCV
|
||||
assert_almost_equal(m[0,0], 72.0)
|
||||
assert_almost_equal(m[0,1], 408.0)
|
||||
assert_almost_equal(m[0,2], 2748.0)
|
||||
assert_almost_equal(m[0,3], 19776.0)
|
||||
assert_almost_equal(m[1,0], 680.0)
|
||||
assert_almost_equal(m[1,1], 3766.0)
|
||||
assert_almost_equal(m[1,2], 24836.0)
|
||||
assert_almost_equal(m[2,0], 7682.0)
|
||||
assert_almost_equal(m[2,1], 43882.0)
|
||||
assert_almost_equal(m[3,0], 95588.0)
|
||||
|
||||
def test_normalized_moments():
|
||||
nu = regionprops(SAMPLE, ['NormalizedMoments'])[0]['NormalizedMoments']
|
||||
#: determined with OpenCV
|
||||
assert_almost_equal(nu[0,2], 0.08410493827160502)
|
||||
assert_almost_equal(nu[1,1], -0.016846707818929982)
|
||||
assert_almost_equal(nu[1,2], -0.002899800614433943)
|
||||
assert_almost_equal(nu[2,0], 0.24301268861454037)
|
||||
assert_almost_equal(nu[2,1], 0.045473992910668816)
|
||||
assert_almost_equal(nu[3,0], -0.017278118992041805)
|
||||
|
||||
def test_orientation():
|
||||
orientation = regionprops(SAMPLE, ['Orientation'])[0]['Orientation']
|
||||
# determined with MATLAB
|
||||
assert_almost_equal(orientation, 0.10446844651921)
|
||||
|
||||
def test_solidity():
|
||||
solidity = regionprops(SAMPLE, ['Solidity'])[0]['Solidity']
|
||||
# determined with MATLAB
|
||||
assert_almost_equal(solidity, 0.580645161290323)
|
||||
|
||||
def test_weighted_central_moments():
|
||||
wmu = regionprops(SAMPLE, ['WeightedCentralMoments'], INTENSITY_SAMPLE
|
||||
)[0]['WeightedCentralMoments']
|
||||
ref = np.array(
|
||||
[[ 7.4000000000e+01, -2.1316282073e-13, 4.7837837838e+02,
|
||||
-7.5943608473e+02],
|
||||
[ 3.7303493627e-14, -8.7837837838e+01, -1.4801314828e+02,
|
||||
-1.2714707125e+03],
|
||||
[ 1.2602837838e+03, 2.1571526662e+03, 6.6989799420e+03,
|
||||
1.5304076361e+04],
|
||||
[ -7.6561796932e+02, -4.2385971907e+03, -9.9501164076e+03,
|
||||
-3.3156729271e+04]]
|
||||
)
|
||||
np.set_printoptions(precision=10)
|
||||
assert_array_almost_equal(wmu, ref)
|
||||
|
||||
def test_weighted_centroid():
|
||||
centroid = regionprops(SAMPLE, ['WeightedCentroid'], INTENSITY_SAMPLE
|
||||
)[0]['WeightedCentroid']
|
||||
assert_array_almost_equal(centroid, (5.540540540540, 9.445945945945))
|
||||
|
||||
def test_weighted_hu_moments():
|
||||
whu = regionprops(SAMPLE, ['WeightedHuMoments'], INTENSITY_SAMPLE
|
||||
)[0]['WeightedHuMoments']
|
||||
ref = np.array([
|
||||
3.1750587329e-01,
|
||||
2.1417517159e-02,
|
||||
2.3609322038e-02,
|
||||
1.2565683360e-03,
|
||||
8.3014209421e-07,
|
||||
-3.5073773473e-05,
|
||||
6.7936409056e-06
|
||||
])
|
||||
assert_array_almost_equal(whu, ref)
|
||||
|
||||
def test_weighted_moments():
|
||||
wm = regionprops(SAMPLE, ['WeightedMoments'], INTENSITY_SAMPLE
|
||||
)[0]['WeightedMoments']
|
||||
ref = np.array(
|
||||
[[ 7.4000000000e+01, 4.1000000000e+02, 2.7500000000e+03,
|
||||
1.9778000000e+04],
|
||||
[ 6.9900000000e+02, 3.7850000000e+03, 2.4855000000e+04,
|
||||
1.7500100000e+05],
|
||||
[ 7.8630000000e+03, 4.4063000000e+04, 2.9347700000e+05,
|
||||
2.0810510000e+06],
|
||||
[ 9.7317000000e+04, 5.7256700000e+05, 3.9007170000e+06,
|
||||
2.8078871000e+07]]
|
||||
)
|
||||
assert_array_almost_equal(wm, ref)
|
||||
|
||||
def test_weighted_normalized_moments():
|
||||
wnu = regionprops(SAMPLE, ['WeightedNormalizedMoments'], INTENSITY_SAMPLE
|
||||
)[0]['WeightedNormalizedMoments']
|
||||
ref = np.array(
|
||||
[[ np.nan, np.nan, 0.0873590903, -0.0161217406],
|
||||
[ np.nan, -0.0160405109, -0.0031421072, -0.0031376984],
|
||||
[ 0.230146783, 0.0457932622, 0.0165315478, 0.0043903193],
|
||||
[-0.0162529732, -0.0104598869, -0.0028544152, -0.0011057191]]
|
||||
)
|
||||
assert_array_almost_equal(wnu, ref)
|
||||
|
||||
if __name__ == "__main__":
|
||||
from numpy.testing import run_module_suite
|
||||
run_module_suite()
|
||||
@@ -0,0 +1,58 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_equal
|
||||
|
||||
from skimage.measure import structural_similarity as ssim
|
||||
import scipy.optimize as opt
|
||||
|
||||
def test_ssim_patch_range():
|
||||
N = 51
|
||||
X = (np.random.random((N, N)) * 255).astype(np.uint8)
|
||||
Y = (np.random.random((N, N)) * 255).astype(np.uint8)
|
||||
|
||||
assert(ssim(X, Y, win_size=N) < 0.1)
|
||||
assert_equal(ssim(X, X, win_size=N), 1)
|
||||
|
||||
def test_ssim_image():
|
||||
N = 100
|
||||
X = (np.random.random((N, N)) * 255).astype(np.uint8)
|
||||
Y = (np.random.random((N, N)) * 255).astype(np.uint8)
|
||||
|
||||
S0 = ssim(X, X, win_size=3)
|
||||
assert_equal(S0, 1)
|
||||
|
||||
S1 = ssim(X, Y, win_size=3)
|
||||
assert(S1 < 0.3)
|
||||
|
||||
## Come up with a better way of testing the gradient
|
||||
##
|
||||
## def test_ssim_grad():
|
||||
## N = 30
|
||||
## X = np.random.random((N, N)) * 255
|
||||
## Y = np.random.random((N, N)) * 255
|
||||
|
||||
## def func(Y):
|
||||
## return ssim(X, Y, dynamic_range=255)
|
||||
|
||||
## def grad(Y):
|
||||
## return ssim(X, Y, dynamic_range=255, gradient=True)[1]
|
||||
|
||||
## assert(np.all(opt.check_grad(func, grad, Y) < 0.05))
|
||||
|
||||
def test_ssim_dtype():
|
||||
N = 30
|
||||
X = np.random.random((N, N))
|
||||
Y = np.random.random((N, N))
|
||||
|
||||
S1 = ssim(X, Y)
|
||||
|
||||
X = (X * 255).astype(np.uint8)
|
||||
Y = (X * 255).astype(np.uint8)
|
||||
|
||||
S2 = ssim(X, Y)
|
||||
|
||||
assert S1 < 0.1
|
||||
assert S2 < 0.1
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
np.testing.run_module_suite()
|
||||
@@ -155,6 +155,9 @@ def label(np.ndarray[DTYPE_t, ndim=2] input,
|
||||
raise ValueError('Neighbors must be either 4 or 8.')
|
||||
|
||||
# Initialize the first row
|
||||
if data[0, 0] == background:
|
||||
link_bg(forest_p, 0, &background_node)
|
||||
|
||||
for j in range(1, cols):
|
||||
if data[0, j] == background:
|
||||
link_bg(forest_p, j, &background_node)
|
||||
|
||||
@@ -61,5 +61,26 @@ class TestConnectedComponents:
|
||||
[0, 0, 1],
|
||||
[-1, -1, -1]])
|
||||
|
||||
def test_background_two_regions(self):
|
||||
x = np.array([[0, 0, 6],
|
||||
[0, 0, 6],
|
||||
[5, 5, 5]])
|
||||
|
||||
assert_array_equal(label(x, background=0),
|
||||
[[-1, -1, 0],
|
||||
[-1, -1, 0],
|
||||
[ 1, 1, 1]])
|
||||
|
||||
def test_background_one_region_center(self):
|
||||
x = np.array([[0, 0, 0],
|
||||
[0, 1, 0],
|
||||
[0, 0, 0]])
|
||||
|
||||
assert_array_equal(label(x, neighbors=4, background=0),
|
||||
[[-1, -1, -1],
|
||||
[-1, 0, -1],
|
||||
[-1, -1, -1]])
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
run_module_suite()
|
||||
|
||||
@@ -4,3 +4,5 @@ from .finite_radon_transform import *
|
||||
from .project import *
|
||||
from ._project import homography as fast_homography
|
||||
from .integral import *
|
||||
from ._warp import warp
|
||||
from ._warp_zoo import swirl
|
||||
|
||||
@@ -0,0 +1,113 @@
|
||||
__all__ = ['warp']
|
||||
|
||||
import numpy as np
|
||||
from scipy import ndimage
|
||||
from skimage.util import img_as_float
|
||||
|
||||
|
||||
def _stackcopy(a, b):
|
||||
"""Copy b into each color layer of a, such that::
|
||||
|
||||
a[:,:,0] = a[:,:,1] = ... = b
|
||||
|
||||
Parameters
|
||||
----------
|
||||
a : (M, N) or (M, N, P) ndarray
|
||||
Target array.
|
||||
b : (M, N)
|
||||
Source array.
|
||||
|
||||
Notes
|
||||
-----
|
||||
Color images are stored as an ``MxNx3`` or ``MxNx4`` arrays.
|
||||
|
||||
"""
|
||||
if a.ndim == 3:
|
||||
a[:] = b[:, :, np.newaxis]
|
||||
else:
|
||||
a[:] = b
|
||||
|
||||
|
||||
def warp(image, reverse_map, map_args={},
|
||||
output_shape=None, order=1, mode='constant', cval=0.):
|
||||
"""Warp an image according to a given coordinate transformation.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : 2-D array
|
||||
Input image.
|
||||
reverse_map : callable xy = f(xy, **kwargs)
|
||||
Reverse coordinate map. A function that transforms a Px2 array of
|
||||
``(x, y)`` coordinates in the *output image* into their corresponding
|
||||
coordinates in the *source image*. Also see examples below.
|
||||
map_args : dict, optional
|
||||
Keyword arguments passed to `reverse_map`.
|
||||
output_shape : tuple (rows, cols)
|
||||
Shape of the output image generated.
|
||||
order : int
|
||||
Order of splines used in interpolation. See
|
||||
`scipy.ndimage.map_coordinates` for detail.
|
||||
mode : string
|
||||
How to handle values outside the image borders. See
|
||||
`scipy.ndimage.map_coordinates` for detail.
|
||||
cval : string
|
||||
Used in conjunction with mode 'constant', the value outside
|
||||
the image boundaries.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Shift an image to the right:
|
||||
|
||||
>>> from skimage import data
|
||||
>>> image = data.camera()
|
||||
>>>
|
||||
>>> def shift_right(xy):
|
||||
... xy[:, 0] -= 10
|
||||
... return xy
|
||||
>>>
|
||||
>>> warp(image, shift_right)
|
||||
|
||||
"""
|
||||
if image.ndim < 2:
|
||||
raise ValueError("Input must have more than 1 dimension.")
|
||||
|
||||
image = np.atleast_3d(img_as_float(image))
|
||||
ishape = np.array(image.shape)
|
||||
bands = ishape[2]
|
||||
|
||||
if output_shape is None:
|
||||
output_shape = ishape
|
||||
|
||||
coords = np.empty(np.r_[3, output_shape], dtype=float)
|
||||
|
||||
## Construct transformed coordinates
|
||||
|
||||
rows, cols = output_shape[:2]
|
||||
|
||||
# Reshape grid coordinates into a (P, 2) array of (x, y) pairs
|
||||
tf_coords = np.indices((cols, rows), dtype=float).reshape(2, -1).T
|
||||
|
||||
# Map each (x, y) pair to the source image according to
|
||||
# the user-provided mapping
|
||||
tf_coords = reverse_map(tf_coords, **map_args)
|
||||
|
||||
# Reshape back to a (2, M, N) coordinate grid
|
||||
tf_coords = tf_coords.T.reshape((-1, cols, rows)).swapaxes(1, 2)
|
||||
|
||||
# Place the y-coordinate mapping
|
||||
_stackcopy(coords[1, ...], tf_coords[0, ...])
|
||||
|
||||
# Place the x-coordinate mapping
|
||||
_stackcopy(coords[0, ...], tf_coords[1, ...])
|
||||
|
||||
# colour-coordinate mapping
|
||||
coords[2, ...] = range(bands)
|
||||
|
||||
# Prefilter not necessary for order 1 interpolation
|
||||
prefilter = order > 1
|
||||
mapped = 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
|
||||
return np.clip(mapped.squeeze(), 0, 1)
|
||||
@@ -0,0 +1,74 @@
|
||||
from __future__ import division
|
||||
import numpy as np
|
||||
|
||||
from ._warp import warp
|
||||
|
||||
|
||||
def _swirl_mapping(xy, center, rotation, strength, radius):
|
||||
x, y = xy.T
|
||||
x0, y0 = center
|
||||
rho = np.sqrt((x - x0)**2 + (y - y0)**2)
|
||||
|
||||
# Ensure that the transformation decays to approximately 1/1000-th
|
||||
# within the specified radius.
|
||||
radius = radius / 5 * np.log(2)
|
||||
|
||||
theta = rotation + strength * \
|
||||
np.exp(-rho / radius) + \
|
||||
np.arctan2(y - y0, x - x0)
|
||||
|
||||
xy[..., 0] = x0 + rho * np.cos(theta)
|
||||
xy[..., 1] = y0 + rho * np.sin(theta)
|
||||
|
||||
return xy
|
||||
|
||||
def swirl(image, center=None, strength=1, radius=100, rotation=0,
|
||||
output_shape=None, order=1, mode='constant', cval=0):
|
||||
"""Perform a swirl transformation.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input image.
|
||||
center : (x,y) tuple or (2,) ndarray
|
||||
Center coordinate of transformation.
|
||||
strength : float
|
||||
The amount of swirling applied.
|
||||
radius : float
|
||||
The extent of the swirl in pixels. The effect dies out
|
||||
rapidly beyond `radius`.
|
||||
rotation : float
|
||||
Additional rotation applied to the image.
|
||||
|
||||
Returns
|
||||
-------
|
||||
swirled : ndarray
|
||||
Swirled version of the input.
|
||||
|
||||
Other parameters
|
||||
----------------
|
||||
output_shape : tuple or ndarray
|
||||
Size of the generated output image.
|
||||
order : int
|
||||
Order of splines used in interpolation. See
|
||||
`scipy.ndimage.map_coordinates` for detail.
|
||||
mode : string
|
||||
How to handle values outside the image borders. See
|
||||
`scipy.ndimage.map_coordinates` for detail.
|
||||
cval : string
|
||||
Used in conjunction with mode 'constant', the value outside
|
||||
the image boundaries.
|
||||
|
||||
"""
|
||||
|
||||
if center is None:
|
||||
center = np.array(image.shape)[:2] / 2
|
||||
|
||||
warp_args = {'center': center,
|
||||
'rotation': rotation,
|
||||
'strength': strength,
|
||||
'radius': radius}
|
||||
|
||||
return warp(image, _swirl_mapping, map_args=warp_args,
|
||||
output_shape=output_shape,
|
||||
order=order, mode=mode, cval=cval)
|
||||
@@ -4,18 +4,12 @@
|
||||
|
||||
import numpy as np
|
||||
from scipy.ndimage import interpolation as ndii
|
||||
from ._warp import _stackcopy
|
||||
|
||||
__all__ = ['homography']
|
||||
|
||||
eps = np.finfo(float).eps
|
||||
|
||||
def _stackcopy(a, b):
|
||||
"""a[:,:,0] = a[:,:,1] = ... = b"""
|
||||
if a.ndim == 3:
|
||||
a.transpose().swapaxes(1, 2)[:] = b
|
||||
else:
|
||||
a[:] = b
|
||||
|
||||
def homography(image, H, output_shape=None, order=1,
|
||||
mode='constant', cval=0.):
|
||||
"""Perform a projective transformation (homography) on an image.
|
||||
@@ -106,6 +100,8 @@ def homography(image, H, output_shape=None, order=1,
|
||||
|
||||
coords = np.empty(np.r_[3, output_shape], dtype=float)
|
||||
|
||||
# TODO: Refactor this method to use transform.warp instead.
|
||||
|
||||
# Construct transformed coordinates
|
||||
rows, cols = output_shape[:2]
|
||||
rows, cols = np.mgrid[:rows, :cols]
|
||||
|
||||
@@ -1,23 +1,23 @@
|
||||
"""
|
||||
radon.py - Radon and inverse radon transforms
|
||||
|
||||
Based on code of Justin K. Romberg
|
||||
Based on code of Justin K. Romberg
|
||||
(http://www.clear.rice.edu/elec431/projects96/DSP/bpanalysis.html)
|
||||
J. Gillam and Chris Griffin.
|
||||
|
||||
References:
|
||||
References:
|
||||
-B.R. Ramesh, N. Srinivasa, K. Rajgopal, "An Algorithm for Computing
|
||||
the Discrete Radon Transform With Some Applications", Proceedings of
|
||||
the Fourth IEEE Region 10 International Conference, TENCON '89, 1989.
|
||||
-A. C. Kak, Malcolm Slaney, "Principles of Computerized Tomographic
|
||||
Imaging", IEEE Press 1988.
|
||||
"""
|
||||
|
||||
from __future__ import division
|
||||
import numpy as np
|
||||
from scipy.fftpack import fftshift, fft, ifft
|
||||
from ._project import homography
|
||||
|
||||
__all__ = ["radon", "iradon"]
|
||||
__all__ = ["radon", "iradon"]
|
||||
|
||||
|
||||
def radon(image, theta=None):
|
||||
@@ -31,7 +31,7 @@ def radon(image, theta=None):
|
||||
Input image.
|
||||
theta : array_like, dtype=float, optional (default np.arange(180))
|
||||
Projection angles (in degrees).
|
||||
|
||||
|
||||
Returns
|
||||
-------
|
||||
output : ndarray
|
||||
@@ -41,7 +41,7 @@ def radon(image, theta=None):
|
||||
if image.ndim != 2:
|
||||
raise ValueError('The input image must be 2-D')
|
||||
if theta == None:
|
||||
theta = np.arange(180)
|
||||
theta = np.arange(180)
|
||||
height, width = image.shape
|
||||
diagonal = np.sqrt(height ** 2 + width ** 2)
|
||||
heightpad = np.ceil(diagonal - height)
|
||||
@@ -57,7 +57,7 @@ def radon(image, theta=None):
|
||||
out = np.zeros((max(padded_image.shape), len(theta)))
|
||||
|
||||
h, w = padded_image.shape
|
||||
dh, dw = h / 2, w / 2
|
||||
dh, dw = h // 2, w // 2
|
||||
shift0 = np.array([[1, 0, -dw],
|
||||
[0, 1, -dh],
|
||||
[0, 0, 1]])
|
||||
@@ -92,7 +92,7 @@ def iradon(radon_image, theta=None, output_size=None,
|
||||
|
||||
Reconstruct an image from the radon transform, using the filtered
|
||||
back projection algorithm.
|
||||
|
||||
|
||||
Parameters
|
||||
----------
|
||||
radon_image : array_like, dtype=float
|
||||
@@ -110,7 +110,7 @@ def iradon(radon_image, theta=None, output_size=None,
|
||||
interpolation : str, optional (default linear)
|
||||
Interpolation method used in reconstruction.
|
||||
Methods available: nearest, linear.
|
||||
|
||||
|
||||
Returns
|
||||
-------
|
||||
output : ndarray
|
||||
@@ -121,19 +121,19 @@ def iradon(radon_image, theta=None, output_size=None,
|
||||
It applies the fourier slice theorem to reconstruct an image by
|
||||
multiplying the frequency domain of the filter with the FFT of the
|
||||
projection data. This algorithm is called filtered back projection.
|
||||
|
||||
|
||||
"""
|
||||
if radon_image.ndim != 2:
|
||||
raise ValueError('The input image must be 2-D')
|
||||
if theta == None:
|
||||
m, n = radon_image.shape
|
||||
theta = np.linspace(0, 180, n, endpoint=False)
|
||||
th = (np.pi / 180.0) * theta
|
||||
th = (np.pi / 180.0) * theta
|
||||
# if output size not specified, estimate from input radon image
|
||||
if not output_size:
|
||||
output_size = int(np.floor(np.sqrt((radon_image.shape[0]) ** 2 / 2.0)))
|
||||
n = radon_image.shape[0]
|
||||
|
||||
|
||||
img = radon_image.copy()
|
||||
# resize image to next power of two for fourier analysis
|
||||
# speeds up fourier and lessens artifacts
|
||||
@@ -142,7 +142,7 @@ def iradon(radon_image, theta=None, output_size=None,
|
||||
img.resize((order, img.shape[1]))
|
||||
# construct the fourier filter
|
||||
freqs = np.zeros((order, 1))
|
||||
|
||||
|
||||
f = fftshift(abs(np.mgrid[-1:1:2 / order])).reshape(-1, 1)
|
||||
w = 2 * np.pi * f
|
||||
# start from first element to avoid divide by zero
|
||||
@@ -160,8 +160,8 @@ def iradon(radon_image, theta=None, output_size=None,
|
||||
f[1:] = 1
|
||||
else:
|
||||
raise ValueError("Unknown filter: %s" % filter)
|
||||
|
||||
filter_ft = np.tile(f, (1, len(theta)))
|
||||
|
||||
filter_ft = np.tile(f, (1, len(theta)))
|
||||
# apply filter in fourier domain
|
||||
projection = fft(img, axis=0) * filter_ft
|
||||
radon_filtered = np.real(ifft(projection, axis=0))
|
||||
@@ -169,15 +169,15 @@ def iradon(radon_image, theta=None, output_size=None,
|
||||
radon_filtered = radon_filtered[:radon_image.shape[0], :]
|
||||
reconstructed = np.zeros((output_size, output_size))
|
||||
mid_index = np.ceil(n / 2.0)
|
||||
|
||||
|
||||
x = output_size
|
||||
y = output_size
|
||||
[X, Y] = np.mgrid[0.0:x, 0.0:y]
|
||||
xpr = X - int(output_size) / 2
|
||||
ypr = Y - int(output_size) / 2
|
||||
xpr = X - int(output_size) // 2
|
||||
ypr = Y - int(output_size) // 2
|
||||
|
||||
# reconstruct image by interpolation
|
||||
if interpolation == "nearest":
|
||||
if interpolation == "nearest":
|
||||
for i in range(len(theta)):
|
||||
k = np.round(mid_index + xpr * np.sin(th[i]) - ypr * np.cos(th[i]))
|
||||
reconstructed += radon_filtered[
|
||||
@@ -186,12 +186,12 @@ def iradon(radon_image, theta=None, output_size=None,
|
||||
for i in range(len(theta)):
|
||||
t = xpr*np.sin(th[i]) - ypr*np.cos(th[i])
|
||||
a = np.floor(t)
|
||||
b = mid_index + a
|
||||
b = mid_index + a
|
||||
b0 = ((((b + 1 > 0) & (b + 1 < n)) * (b + 1)) - 1).astype(np.int)
|
||||
b1 = ((((b > 0) & (b < n)) * b) - 1).astype(np.int)
|
||||
reconstructed += (t - a) * radon_filtered[b0, i] + \
|
||||
(a - t + 1) * radon_filtered[b1, i]
|
||||
else:
|
||||
raise ValueError("Unknown interpolation: %s" % interpolation)
|
||||
|
||||
|
||||
return reconstructed * np.pi / (2 * len(th))
|
||||
|
||||
@@ -1,7 +1,7 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_array_almost_equal
|
||||
|
||||
from skimage.transform.project import _stackcopy
|
||||
from skimage.transform._warp import _stackcopy
|
||||
from skimage.transform import homography, fast_homography
|
||||
from skimage import data
|
||||
from skimage.color import rgb2gray
|
||||
|
||||
@@ -0,0 +1,17 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_array_almost_equal
|
||||
|
||||
from skimage import transform as tf, data, img_as_float
|
||||
|
||||
|
||||
def test_roundtrip():
|
||||
image = img_as_float(data.checkerboard())
|
||||
|
||||
swirl_params = {'radius': 80, 'rotation': 0, 'order': 2, 'mode': 'reflect'}
|
||||
swirled = tf.swirl(image, strength=10, **swirl_params)
|
||||
unswirled = tf.swirl(swirled, strength=-10, **swirl_params)
|
||||
|
||||
assert np.mean(np.abs(image - unswirled)) < 0.01
|
||||
|
||||
if __name__ == "__main__":
|
||||
np.testing.run_module_suite()
|
||||
@@ -1 +1,2 @@
|
||||
from .dtype import *
|
||||
from .shape import *
|
||||
|
||||
+145
-79
@@ -10,8 +10,8 @@ dtype_range = {np.uint8: (0, 255),
|
||||
np.uint16: (0, 65535),
|
||||
np.int8: (-128, 127),
|
||||
np.int16: (-32768, 32767),
|
||||
np.float32: (0, 1),
|
||||
np.float64: (0, 1)}
|
||||
np.float32: (-1, 1),
|
||||
np.float64: (-1, 1)}
|
||||
|
||||
integer_types = (np.uint8, np.uint16, np.int8, np.int16)
|
||||
|
||||
@@ -20,19 +20,24 @@ _supported_types = (np.uint8, np.uint16, np.uint32,
|
||||
np.float32, np.float64)
|
||||
|
||||
if np.__version__ >= "1.6.0":
|
||||
dtype_range[np.float16] = (0, 1)
|
||||
dtype_range[np.float16] = (-1, 1)
|
||||
_supported_types += (np.float16, )
|
||||
|
||||
|
||||
def convert(image, dtype, force_copy=False):
|
||||
def convert(image, dtype, force_copy=False, uniform=False):
|
||||
"""
|
||||
Convert an image to the requested data-type.
|
||||
|
||||
Warnings are issued in case of precision loss, or when
|
||||
negative values have to be scaled into the positive domain.
|
||||
Floating point values must be in the range [0.0, 1.0].
|
||||
Warnings are issued in case of precision loss, or when negative values
|
||||
are clipped during conversion to unsigned integer types (sign loss).
|
||||
|
||||
Floating point values are expected to be normalized and will be clipped
|
||||
to the range [0.0, 1.0] or [-1.0, 1.0] when converting to unsigned or
|
||||
signed integers respectively.
|
||||
|
||||
Numbers are not shifted to the negative side when converting from
|
||||
floating point or unsigned integer types to signed integer types.
|
||||
unsigned to signed integer types. Negative values will be clipped when
|
||||
converting to unsigned integers.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
@@ -42,17 +47,29 @@ def convert(image, dtype, force_copy=False):
|
||||
Target data-type.
|
||||
force_copy : bool
|
||||
Force a copy of the data, irrespective of its current dtype.
|
||||
uniform : bool
|
||||
Uniformly quantize the floating point range to the integer range.
|
||||
By default (uniform=False) floating point values are scaled and
|
||||
rounded to the nearest integers, which minimizes back and forth
|
||||
conversion errors.
|
||||
|
||||
References
|
||||
----------
|
||||
(1) DirectX data conversion rules.
|
||||
http://msdn.microsoft.com/en-us/library/windows/desktop/dd607323%28v=vs.85%29.aspx
|
||||
(2) Data Conversions.
|
||||
In "OpenGL ES 2.0 Specification v2.0.25", pp 7-8. Khronos Group, 2010.
|
||||
(3) Proper treatment of pixels as integers. A.W. Paeth.
|
||||
In "Graphics Gems I", pp 249-256. Morgan Kaufmann, 1990.
|
||||
(4) Dirty Pixels. J. Blinn.
|
||||
In "Jim Blinn's corner: Dirty Pixels", pp 47-57. Morgan Kaufmann, 1998.
|
||||
|
||||
"""
|
||||
image = np.asarray(image)
|
||||
dtype = np.dtype(dtype).type
|
||||
dtype_in = image.dtype.type
|
||||
dtypeobj = np.dtype(dtype)
|
||||
dtypeobj_in = np.dtype(dtype_in)
|
||||
kind = dtypeobj.kind
|
||||
kind_in = dtypeobj_in.kind
|
||||
itemsize = dtypeobj.itemsize
|
||||
itemsize_in = dtypeobj_in.itemsize
|
||||
dtypeobj_in = image.dtype
|
||||
dtype = dtypeobj.type
|
||||
dtype_in = dtypeobj_in.type
|
||||
|
||||
if dtype_in == dtype:
|
||||
if force_copy:
|
||||
@@ -74,93 +91,142 @@ def convert(image, dtype, force_copy=False):
|
||||
# Return first of `dtypes` with itemsize greater than `itemsize`
|
||||
return next(dt for dt in dtypes if itemsize < np.dtype(dt).itemsize)
|
||||
|
||||
def _dtype2(kind, bits, itemsize=1):
|
||||
# Return dtype of `kind` that can store a `bits` wide unsigned int
|
||||
c = lambda x, y: x <= y if kind == 'u' else x < y
|
||||
s = next(i for i in (itemsize, ) + (2, 4, 8) if c(bits, i*8))
|
||||
return np.dtype(kind + str(s))
|
||||
|
||||
def _scale(a, n, m, copy=True):
|
||||
# Scale unsigned/positive integers from n to m bits
|
||||
# Numbers can be represented exactly only if m is a multiple of n
|
||||
# Output array is of same kind as input.
|
||||
kind = a.dtype.kind
|
||||
if n == m:
|
||||
return a.copy() if copy else a
|
||||
elif n > m:
|
||||
# downscale with precision loss
|
||||
prec_loss()
|
||||
if copy:
|
||||
b = np.empty(a.shape, _dtype2(kind, m))
|
||||
np.floor_divide(a, 2**(n - m), out=b, dtype=a.dtype,
|
||||
casting='unsafe')
|
||||
return b
|
||||
else:
|
||||
a //= 2**(n - m)
|
||||
return a
|
||||
elif m % n == 0:
|
||||
# exact upscale to a multiple of n bits
|
||||
if copy:
|
||||
b = np.empty(a.shape, _dtype2(kind, m))
|
||||
np.multiply(a, (2**m - 1) // (2**n - 1), out=b, dtype=b.dtype)
|
||||
return b
|
||||
else:
|
||||
a = np.array(a, _dtype2(kind, m, a.dtype.itemsize), copy=False)
|
||||
a *= (2**m - 1) // (2**n - 1)
|
||||
return a
|
||||
else:
|
||||
# upscale to a multiple of n bits,
|
||||
# then downscale with precision loss
|
||||
prec_loss()
|
||||
o = (m // n + 1) * n
|
||||
if copy:
|
||||
b = np.empty(a.shape, _dtype2(kind, o))
|
||||
np.multiply(a, (2**o - 1) // (2**n - 1), out=b, dtype=b.dtype)
|
||||
b //= 2**(o - m)
|
||||
return b
|
||||
else:
|
||||
a = np.array(a, _dtype2(kind, o, a.dtype.itemsize), copy=False)
|
||||
a *= (2**o - 1) // (2**n - 1)
|
||||
a //= 2**(o - m)
|
||||
return a
|
||||
|
||||
kind = dtypeobj.kind
|
||||
kind_in = dtypeobj_in.kind
|
||||
itemsize = dtypeobj.itemsize
|
||||
itemsize_in = dtypeobj_in.itemsize
|
||||
if kind in 'ui':
|
||||
imin = np.iinfo(dtype).min
|
||||
imax = np.iinfo(dtype).max
|
||||
if kind_in in 'ui':
|
||||
imin_in = np.iinfo(dtype_in).min
|
||||
imax_in = np.iinfo(dtype_in).max
|
||||
|
||||
if kind_in == 'f':
|
||||
if np.min(image) < 0 or np.max(image) > 1:
|
||||
raise ValueError("Images of type float must be between 0 and 1")
|
||||
if np.min(image) < -1.0 or np.max(image) > 1.0:
|
||||
raise ValueError("Images of type float must be between -1 and 1.")
|
||||
if kind == 'f':
|
||||
# floating point -> floating point
|
||||
if itemsize_in > itemsize:
|
||||
prec_loss()
|
||||
return dtype(image)
|
||||
|
||||
# floating point -> integer
|
||||
prec_loss()
|
||||
# use float type that can represent output integer type
|
||||
image = np.array(image, _dtype(itemsize, dtype_in,
|
||||
np.float32, np.float64))
|
||||
image *= np.iinfo(dtype).max + 1
|
||||
np.clip(image, 0, np.iinfo(dtype).max, out=image)
|
||||
if not uniform:
|
||||
if kind == 'u':
|
||||
image *= imax
|
||||
else:
|
||||
image *= imax - imin
|
||||
image -= 1.0
|
||||
image /= 2.0
|
||||
np.rint(image, out=image)
|
||||
np.clip(image, imin, imax, out=image)
|
||||
elif kind == 'u':
|
||||
image *= imax + 1
|
||||
np.clip(image, 0, imax, out=image)
|
||||
else:
|
||||
image *= (imax - imin + 1.0) / 2.0
|
||||
np.floor(image, out=image)
|
||||
np.clip(image, imin, imax, out=image)
|
||||
return dtype(image)
|
||||
|
||||
if kind == 'f':
|
||||
# integer -> floating point
|
||||
if itemsize_in >= itemsize:
|
||||
prec_loss()
|
||||
# use float type that can represent input integers
|
||||
# use float type that can exactly represent input integers
|
||||
image = np.array(image, _dtype(itemsize_in, dtype,
|
||||
np.float32, np.float64))
|
||||
if np.iinfo(dtype_in).min:
|
||||
sign_loss()
|
||||
image -= np.iinfo(dtype_in).min
|
||||
image /= np.iinfo(dtype_in).max - np.iinfo(dtype_in).min
|
||||
return dtype(image)
|
||||
if kind_in == 'u':
|
||||
# unsigned integer -> integer
|
||||
shift = 1 if kind == 'i' else 0
|
||||
if itemsize_in > itemsize:
|
||||
prec_loss()
|
||||
image = image >> 8 * (itemsize_in - itemsize) + shift
|
||||
return dtype(image)
|
||||
result = dtype(image)
|
||||
result <<= 8 * (itemsize - itemsize_in) - shift
|
||||
if itemsize - itemsize_in == 3:
|
||||
# uint8 -> (u)int32
|
||||
# hint: 4294967295 == (255 << 24) + (255 << 16) + (255 << 8) + 255
|
||||
image = dtype(image)
|
||||
image *= 2**16 + 2**8 + 1
|
||||
if shift:
|
||||
result += image >> shift
|
||||
if kind_in == 'u':
|
||||
image /= imax_in
|
||||
# DirectX uses this conversion also for signed ints
|
||||
#if imin_in:
|
||||
# np.maximum(image, -1.0, out=image)
|
||||
else:
|
||||
result += image
|
||||
return dtype(result)
|
||||
image *= 2.0
|
||||
image += 1.0
|
||||
image /= imax_in - imin_in
|
||||
return dtype(image)
|
||||
|
||||
if kind_in == 'u':
|
||||
if kind == 'i':
|
||||
# unsigned integer -> signed integer
|
||||
image = _scale(image, 8*itemsize_in, 8*itemsize-1)
|
||||
return image.view(dtype)
|
||||
else:
|
||||
# unsigned integer -> unsigned integer
|
||||
return _scale(image, 8*itemsize_in, 8*itemsize)
|
||||
|
||||
if kind == 'u':
|
||||
# signed integer -> unsigned integer
|
||||
sign_loss()
|
||||
# use next higher precision signed integer type
|
||||
image = np.array(image, _dtype(itemsize_in,
|
||||
np.int16, np.int32, np.int64))
|
||||
image -= np.iinfo(dtype_in).min
|
||||
if itemsize_in == itemsize:
|
||||
return dtype(image)
|
||||
if itemsize_in > itemsize:
|
||||
prec_loss()
|
||||
image >>= 8 * (itemsize_in - itemsize)
|
||||
return dtype(image)
|
||||
result = dtype(image)
|
||||
result <<= 8 * (itemsize - itemsize_in)
|
||||
if itemsize - itemsize_in == 3:
|
||||
# int8 -> uint32
|
||||
image = dtype(image)
|
||||
image *= 2**16 + 2**8 + 1
|
||||
result += dtype(image)
|
||||
image = _scale(image, 8*itemsize_in-1, 8*itemsize)
|
||||
result = np.empty(image.shape, dtype)
|
||||
np.maximum(image, 0, out=result, dtype=image.dtype, casting='unsafe')
|
||||
return result
|
||||
if kind == 'i':
|
||||
# signed integer -> signed integer
|
||||
if itemsize_in > itemsize:
|
||||
prec_loss()
|
||||
return dtype(image // 2**(8 * (itemsize_in - itemsize)))
|
||||
# use next higher precision signed integer type
|
||||
image = np.array(image, _dtype(itemsize_in,
|
||||
np.int16, np.int32, np.int64))
|
||||
image -= np.iinfo(dtype_in).min
|
||||
# use next higher precision signed integer type
|
||||
result = np.array(image, _dtype(image.itemsize, np.int32, np.int64))
|
||||
result *= 2**(8 * (itemsize - itemsize_in))
|
||||
if itemsize - itemsize_in == 3:
|
||||
# int8 -> int32
|
||||
image = dtype(image)
|
||||
image *= 2**16 + 2**8 + 1
|
||||
result += image
|
||||
result += np.iinfo(dtype).min
|
||||
return dtype(result)
|
||||
|
||||
# signed integer -> signed integer
|
||||
if itemsize_in > itemsize:
|
||||
return _scale(image, 8*itemsize_in-1, 8*itemsize-1)
|
||||
image = image.astype(_dtype2('i', itemsize*8))
|
||||
image -= imin_in
|
||||
image = _scale(image, 8*itemsize_in, 8*itemsize, copy=False)
|
||||
image += imin
|
||||
return dtype(image)
|
||||
|
||||
|
||||
def img_as_float(image, force_copy=False):
|
||||
|
||||
@@ -9,12 +9,13 @@ dtype_range = {np.uint8: (0, 255),
|
||||
np.uint16: (0, 65535),
|
||||
np.int8: (-128, 127),
|
||||
np.int16: (-32768, 32767),
|
||||
np.float32: (0, 1),
|
||||
np.float64: (0, 1)}
|
||||
np.float32: (-1.0, 1.0),
|
||||
np.float64: (-1.0, 1.0)}
|
||||
|
||||
def _verify_range(msg, x, vmin, vmax):
|
||||
def _verify_range(msg, x, vmin, vmax, dtype):
|
||||
assert_equal(x[0], vmin)
|
||||
assert_equal(x[-1], vmax)
|
||||
assert x.dtype == dtype
|
||||
|
||||
def test_range():
|
||||
for dtype in dtype_range:
|
||||
@@ -29,12 +30,13 @@ def test_range():
|
||||
|
||||
omin, omax = dtype_range[dt]
|
||||
|
||||
if imin == 0:
|
||||
if imin == 0 or omin == 0:
|
||||
omin = 0
|
||||
imin = 0
|
||||
|
||||
yield _verify_range, \
|
||||
"From %s to %s" % (np.dtype(dtype), np.dtype(dt)), \
|
||||
y, omin, omax
|
||||
yield (_verify_range,
|
||||
"From %s to %s" % (np.dtype(dtype), np.dtype(dt)),
|
||||
y, omin, omax, np.dtype(dt))
|
||||
|
||||
|
||||
def test_range_extra_dtypes():
|
||||
@@ -57,9 +59,9 @@ def test_range_extra_dtypes():
|
||||
x = np.linspace(imin, imax, 10).astype(dtype_in)
|
||||
y = convert(x, dt)
|
||||
omin, omax = dtype_range_extra[dt]
|
||||
yield _verify_range, \
|
||||
"From %s to %s" % (np.dtype(dtype_in), np.dtype(dt)), \
|
||||
y, omin, omax
|
||||
yield (_verify_range,
|
||||
"From %s to %s" % (np.dtype(dtype_in), np.dtype(dt)),
|
||||
y, omin, omax, np.dtype(dt))
|
||||
|
||||
|
||||
def test_unsupported_dtype():
|
||||
@@ -70,7 +72,7 @@ def test_unsupported_dtype():
|
||||
def test_float_out_of_range():
|
||||
too_high = np.array([2], dtype=np.float32)
|
||||
assert_raises(ValueError, img_as_int, too_high)
|
||||
too_low = np.array([-1], dtype=np.float32)
|
||||
too_low = np.array([-2], dtype=np.float32)
|
||||
assert_raises(ValueError, img_as_int, too_low)
|
||||
|
||||
|
||||
@@ -84,4 +86,3 @@ def test_copy():
|
||||
|
||||
if __name__ == '__main__':
|
||||
np.testing.run_module_suite()
|
||||
|
||||
|
||||
Reference in New Issue
Block a user