wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-13 17:45:20 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' into debian

Browse Source
This commit is contained in:
Stefan van der Walt
2013-10-05 13:18:18 +02:00
parent b6c15cb5b0 547c0bb7d8
commit 512a3900fe
314 changed files with 20136 additions and 6668 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.travis.yml
+15 -2
View File
@@ -1,7 +1,7 @@
# vim ft=yaml
# travis-ci.org definition for skimage build
#
# We pretend to be erlang because we need can't use the python support in
# We pretend to be erlang because we can't use the python support in
# travis-ci; it uses virtualenvs, they do not have numpy, scipy, matplotlib,
# and it is impractical to build them
@@ -19,11 +19,24 @@ install:
- sudo easy_install$PYSUF pip
- sudo pip-$PYVER install cython
- sudo apt-get install libfreeimage3
- if [[ $PYVER == '2.7' ]]; then sudo apt-get install $PYTHON-matplotlib; fi
- if [[ $PYVER == '3.2' ]]; then sudo pip-$PYVER install git+git://github.com/matplotlib/matplotlib.git@v1.2.x; fi
- sudo pip-$PYVER install flake8
- $PYTHON setup.py build
- sudo $PYTHON setup.py install
script:
# Check if setup.py's match bento.info
- $PYTHON check_bento_build.py
# Change into an innocuous directory and find tests from installation
- mkdir $HOME/.matplotlib
- "echo 'backend : Agg' > $HOME/.matplotlib/matplotlibrc"
- "echo 'backend.qt4 : PyQt4' >> $HOME/.matplotlib/matplotlibrc"
- mkdir for_test
- cd for_test
- nosetests-$PYVER --exe -v --cover-package=skimage skimage
# Change back to repository root directory and run all doc examples
- cd ..
- for f in doc/examples/*.py; do $PYTHON "$f"; if [ $? -ne 0 ]; then exit 1; fi done
- for f in doc/examples/applications/*.py; do $PYTHON "$f"; if [ $? -ne 0 ]; then exit 1; fi done
# Run pep8 and flake tests
- flake8 --exit-zero --exclude=test_*,six.py skimage doc/examples viewer_examples
CONTRIBUTING.txt
+180
View File
@@ -0,0 +1,180 @@
Development process
-------------------
Here's the long and short of it:
1. If you are a first-time contributor:
* Go to `https://github.com/scikit-image/scikit-image
<http://github.com/scikit-image/scikit-image>`_ and click the
"fork" button to create your own copy of the project.
* Clone the project to your local computer::
git clone git@github.com:your-username/scikit-image.git
* Add upstream repository::
git remote add upstream git@github.com:scikit-image/scikit-image.git
* Now, you have remote repositories named:
- ``upstream``, which refers to the ``scikit-image`` repository
- ``origin``, which refers to your personal fork
2. Develop your contribution:
* Pull the latest changes from upstream::
git checkout master
git pull upstream master
* Create a branch for the feature you want to work on. Since the
branch name will appear in the merge message, use a sensible name
such as 'transform-speedups'::
git checkout -b transform-speedups
* Commit locally as you progress (``git add`` and ``git commit``)
3. To submit your contribution:
* Push your changes back to your fork on GitHub::
git push origin transform-speedups
* Go to GitHub. The new branch will show up with a Pull Request button -
click it.
* If you want, post on the `mailing list
<http://groups.google.com/group/scikit-image>`_ to explain your changes or
to ask for review.
For a more detailed discussion, read these :doc:`detailed documents
<gitwash/index>` on how to use Git with ``scikit-image``
(`<http://scikit-image.org/docs/dev/gitwash/index.html>`_).
.. note::
To reviewers: add a short explanation of what a branch did to the merge
message and, if closing a bug, also add "Closes gh-123" where 123 is the
bug number.
Divergence between ``upstream master`` and your feature branch
..............................................................
Do *not* ever merge the main branch into yours. If GitHub indicates that the
branch of your Pull Request can no longer be merged automatically, rebase
onto master::
git checkout master
git pull upstream master
git checkout transform-speedups
git rebase master
If any conflicts occur, fix the according files and continue::
git add conflict-file1 conflict-file2
git rebase --continue
However, you should only rebase your own branches and must generally not
rebase any branch which you collaborate on with someone else.
Finally, you must push your rebased branch::
git push --force origin transform-speedups
(If you are curious, here's a further discussion on the
`dangers of rebasing <http://tinyurl.com/lll385>`__.
Also see this `LWN article <http://tinyurl.com/nqcbkj>`__.)
Guidelines
----------
* All code should have tests (see `test coverage`_ below for more details).
* All code should be documented, to the same
`standard <http://projects.scipy.org/numpy/wiki/CodingStyleGuidelines>`_
as NumPy and SciPy.
* For new functionality, always add an example to the
gallery.
* No changes should be committed without review. Ask on the
`mailing list <http://groups.google.com/group/scikit-image>`_ if
you get no response to your pull request.
**Never merge your own pull request.**
* Examples in the gallery should have a maximum figure width of 8 inches.
Stylistic Guidelines
--------------------
* Set up your editor to remove trailing whitespace. Follow `PEP08
<www.python.org/dev/peps/pep-0008/>`__. Check code with pyflakes / flake8.
* Use numpy data types instead of strings (``np.uint8`` instead of
``"uint8"``).
* Use the following import conventions::
import numpy as np
import matplotlib.pyplot as plt
cimport numpy as cnp # in Cython code
* When documenting array parameters, use ``image : (M, N) ndarray``
and then refer to ``M`` and ``N`` in the docstring, if necessary.
* Functions should support all input image dtypes. Use utility functions such
as ``img_as_float`` to help convert to an appropriate type. The output
format can be whatever is most efficient. This allows us to string together
several functions into a pipeline, e.g.::
hough(canny(my_image))
* Use ``Py_ssize_t`` as data type for all indexing, shape and size variables
in C/C++ and Cython code.
Test coverage
-------------
Tests for a module should ideally cover all code in that module,
i.e., statement coverage should be at 100%.
To measure the test coverage, install
`coverage.py <http://nedbatchelder.com/code/coverage/>`__
(using ``easy_install coverage``) and then run::
$ make coverage
This will print a report with one line for each file in `skimage`,
detailing the test coverage::
Name Stmts Exec Cover Missing
------------------------------------------------------------------------------
skimage/color/colorconv 77 77 100%
skimage/filter/__init__ 1 1 100%
...
Activate Travis-CI for your fork (optional)
-------------------------------------------
Travis-CI checks all unittests in the project to prevent breakage.
Before sending a pull request, you may want to check that Travis-CI
successfully passes all tests. To do so,
* Go to `Travis-CI <http://travis-ci.org/>`__ and follow the Sign In link at the top
* Go to your `profile page <https://travis-ci.org/profile>`__ and switch on your
scikit-image fork
It corresponds to steps one and two in
`Travis-CI documentation <http://about.travis-ci.org/docs/user/getting-started/>`__
(Step three is already done in scikit-image).
Thus, as soon as you push your code to your fork, it will trigger Travis-CI,
and you will receive an email notification when the process is done.
Bugs
----
Please `report bugs on GitHub <https://github.com/scikit-image/scikit-image/issues>`_.
CONTRIBUTORS.txt
+22 -4
View File
@@ -48,7 +48,7 @@
Incorporating CellProfiler's Sobel edge detector, build and bug fixes.
Radon transform, template matching.
- Emmanuelle Guillart
- Emmanuelle Gouillart
Total variation noise filtering, integration of CellProfiler's
mathematical morphology tools, random walker segmentation,
tutorials, and more.
@@ -113,7 +113,8 @@
Fixes and tests for Histograms of Oriented Gradients.
- Joshua Warner
Multichannel random walker segmentation.
Multichannel random walker segmentation, unified peak finder backend,
n-dimensional array padding, marching cubes, bug and doc fixes.
- Petter Strandmark
Perimeter calculation in regionprops.
@@ -131,8 +132,25 @@
Dense DAISY feature description, circle perimeter drawing.
- François Boulogne
Andres Method for circle perimeter, ellipse perimeter drawing.
Circular Hough Transform
Drawing: Andres Method for circle perimeter, ellipse perimeter drawing, Bezier curve.
Circular and elliptical Hough Transforms
Various fixes
- Thouis Jones
Vectorized operators for arrays of 16-bit ints.
- Xavier Moles Lopez
Color separation (color deconvolution) for several stainings.
- Jostein Bø Fløystad
Reconstruction circle mode for Radon transform
Simultaneous Algebraic Reconstruction Technique for inverse Radon transform
- Matt Terry
Color difference functions
- Eugene Dvoretsky
Yen threshold implementation.
- Riaan van den Dool
skimage.io plugin: GDAL
DEPENDS.txt
+17 -3
View File
@@ -2,11 +2,15 @@ Build Requirements
------------------
* `Python >= 2.5 <http://python.org>`__
* `Numpy >= 1.6 <http://numpy.scipy.org/>`__
* `Cython >= 0.15 <http://www.cython.org/>`__
* `Cython >= 0.17 <http://www.cython.org/>`__
`Matplotlib >= 1.0 <http://matplotlib.sf.net>`__ is needed to generate the
examples in the documentation.
You can use pip to automatically install the base dependencies as follows::
$ pip install -r requirements.txt
Runtime requirements
--------------------
* `SciPy >= 0.10 <http://scipy.org>`__
@@ -31,10 +35,20 @@ Optional Requirements
You can use this scikit with the basic requirements listed above, but some
functionality is only available with the following installed:
`PyQt4 <http://wiki.python.org/moin/PyQt>`__
* `PyQt4 <http://wiki.python.org/moin/PyQt>`__
The ``qt`` plugin that provides ``imshow(x, fancy=True)`` and `skivi`.
`FreeImage <http://freeimage.sf.net>`__
* `FreeImage <http://freeimage.sf.net>`__
The ``freeimage`` plugin provides support for reading various types of
image file formats, including multi-page TIFFs.
* `PyAMG <http://pyamg.org/>`__
The ``pyamg`` module is used for the fast `cg_mg` mode of random
walker segmentation.
Testing requirements
--------------------
* `Nose <https://nose.readthedocs.org/en/latest/>`__
A Python Unit Testing Framework
* `Coverage.py <http://nedbatchelder.com/code/coverage/>`__
A tool that generates a unit test code coverage report
DEVELOPMENT.txt
-111
View File
@@ -1,111 +0,0 @@
Development process
-------------------
Here's the long and short of it:
* Go to `https://github.com/scikit-image/scikit-image
<http://github.com/scikit-image/scikit-image>`_ and follow the
instructions on making your own fork.
* Create a new branch for the feature you want to work on. Since the
branch name will appear in the merge message, use a sensible name
such as 'transform-speedups'.
* Commit locally as you progress.
* Push your changes back to GitHub and create a Pull Request by
clicking 'Pull Request' in GitHub.
* Optionally, post on the `mailing list <http://groups.google.com/group/scikit-image>`_ to explain your changes.
Read these :doc:`detailed documents <gitwash/index>` on how to use Git with
``scikit-image`` (`<http://scikit-image.org/docs/dev/gitwash/index.html>`_).
.. note::
Do *not* merge the main branch into yours. If GitHub indicates that the
Pull Request can no longer be merged automatically, rebase onto master.
(If you are curious, here's a further discussion on
the `dangers of rebasing <http://tinyurl.com/lll385>`__. Also
see this `LWN article <http://tinyurl.com/nqcbkj>`__.)
* To reviewers: add a short explanation of what a branch did to the merge
message or, if closing a bug, add "Closes gh-XXXX".
You may also read this summary by Fernando Perez of the IPython
project on how they manage to keep review overhead to a minimum:
http://mail.scipy.org/pipermail/ipython-dev/2010-October/006746.html
Guidelines
----------
* All code should have tests (see `test coverage`_ below for more details).
* All code should be documented, to the same
`standard <http://projects.scipy.org/numpy/wiki/CodingStyleGuidelines>`_
as NumPy and SciPy. For new functionality, always add an example to the
gallery.
* Follow the `Python PEPs <http://www.python.org/dev/peps/pep-0008/>`_
where possible.
* No major changes should be committed without review. Ask on the
`mailing list <http://groups.google.com/group/scikit-image>`_ if
you get no response to your pull request.
* Examples in the gallery should have a maximum figure width of 8 inches.
Stylistic Guidelines
--------------------
* Use numpy data types instead of strings (``np.uint8`` instead of
``"uint8"``).
* Use the following import conventions::
import numpy as np
import matplotlib.pyplot as plt
cimport numpy as cnp # in Cython code
* When documenting array parameters, use ``image : (M, N) ndarray``,
``image : (M, N, 3) ndarray`` and then refer to ``M`` and ``N`` in the
docstring.
* Set up your editor to remove trailing whitespace. Follow `PEP08
<www.python.org/dev/peps/pep-0008/>`__. Check code with pyflakes / flake8.
* If a function name, say ``segment(...)``, has the same name as the file in
which it is implemented, name that file ``_segment.py`` so that it can still
be imported. All Cython files start with an underscore, e.g.
``_some_module.pyx``.
* Functions should support all input image dtypes. Use utility functions such
as ``img_as_float`` to help convert to an appropriate type. The output
format can be whatever is most efficient. This allows us to string together
several functions into a pipeline, e.g.::
hough(canny(my_image))
* Use `Py_ssize_t` as data type for all indexing, shape and size variables in
C/C++ and Cython code.
Test coverage
-------------
Tests for a module should ideally cover all code in that module,
i.e., statement coverage should be at 100%.
To measure the test coverage, install
`coverage.py <http://nedbatchelder.com/code/coverage/>`__
(using ``easy_install coverage``) and then run::
$ make coverage
This will print a report with one line for each file in `skimage`,
detailing the test coverage::
Name Stmts Exec Cover Missing
------------------------------------------------------------------------------
skimage/color/colorconv 77 77 100%
skimage/filter/__init__ 1 1 100%
...
Bugs
----
Please `report bugs on GitHub <https://github.com/scikit-image/scikit-image/issues>`_.
RELEASE.txt
+4
View File
@@ -1,6 +1,8 @@
How to make a new release of ``skimage``
========================================
- Check ``TODO.txt`` for any outstanding tasks.
- Update release notes.
- To show a list contributors, run ``doc/release/contributors.sh <commit>``,
@@ -46,6 +48,8 @@ How to make a new release of ``skimage``
- Update stable and development version numbers in
``_templates/sidebar_versions.html``.
- Add release date to ``index.rst`` under "Announcements".
- Add previous stable version documentation path to disallowed paths
in `robots.txt`
- Build using ``make gh-pages``.
- Push upstream: ``git push`` in ``gh-pages``.
TODO.txt
+15
View File
@@ -0,0 +1,15 @@
Version 0.10
------------
* Remove deprecated functions in `skimage.filter.rank.*`
* Remove deprecated parameter `epsilon` of `skimage.viewer.LineProfile`
* Remove backwards-compatability of `skimage.measure.regionprops`
* Remove {`ratio`, `sigma`} deprecation warnings of `skimage.segmentation.slic`
Version 0.9
-----------
* Remove deprecated functions
- `skimage.filter.denoise_tv_chambolle`
- `skimage.morphology.is_local_maximum`
- `skimage.transform.hough`
- `skimage.transform.probabilistic_hough`
- `skimage.transform.hough_peaks`
bento.info
+20 -27
View File
@@ -28,7 +28,6 @@ Classifiers:
Operating System :: Unix,
Operating System :: MacOS
HookFile: bscript
UseBackends: Waf
Library:
@@ -52,6 +51,9 @@ Library:
Extension: skimage.measure._moments
Sources:
skimage/measure/_moments.pyx
Extension: skimage.measure._marching_cubes_cy
Sources:
skimage/measure/_marching_cubes_cy.pyx
Extension: skimage.graph._mcp
Sources:
skimage/graph/_mcp.pyx
@@ -91,6 +93,12 @@ Library:
Extension: skimage.morphology._greyreconstruct
Sources:
skimage/morphology/_greyreconstruct.pyx
Extension: skimage.feature.censure_cy
Sources:
skimage/feature/censure_cy.pyx
Extension: skimage.feature._brief_cy
Sources:
skimage/feature/_brief_cy.pyx
Extension: skimage.feature.corner_cy
Sources:
skimage/feature/corner_cy.pyx
@@ -109,6 +117,9 @@ Library:
Extension: skimage.morphology._skeletonize_cy
Sources:
skimage/morphology/_skeletonize_cy.pyx
Extension: skimage.transform._radon_transform
Sources:
skimage/transform/_radon_transform.pyx
Extension: skimage.transform._warps_cy
Sources:
skimage/transform/_warps_cy.pyx
@@ -121,36 +132,18 @@ Library:
Extension: skimage._shared.geometry
Sources:
skimage/_shared/geometry.pyx
Extension: skimage.filter.rank._core16
Extension: skimage.filter.rank.generic_cy
Sources:
skimage/filter/rank/_core16.pyx
Extension: skimage.filter.rank._crank8
skimage/filter/rank/generic_cy.pyx
Extension: skimage.filter.rank.percentile_cy
Sources:
skimage/filter/rank/_crank8.pyx
Extension: skimage.filter.rank._crank16
skimage/filter/rank/percentile_cy.pyx
Extension: skimage.filter.rank.core_cy
Sources:
skimage/filter/rank/_crank16.pyx
Extension: skimage.filter.rank._core8
skimage/filter/rank/core_cy.pyx
Extension: skimage.filter.rank.bilateral_cy
Sources:
skimage/filter/rank/_core8.pyx
Extension: skimage.filter.rank.rank
Sources:
skimage/filter/rank/rank.pyx
Extension: skimage.filter.rank.bilateral_rank
Sources:
skimage/filter/rank/bilateral_rank.pyx
Extension: skimage.filter.rank._crank16_percentiles
Sources:
skimage/filter/rank/_crank16_percentiles.pyx
Extension: skimage.filter.rank.percentile_rank
Sources:
skimage/filter/rank/percentile_rank.pyx
Extension: skimage.filter.rank._crank8_percentiles
Sources:
skimage/filter/rank/_crank8_percentiles.pyx
Extension: skimage.filter.rank._crank16_bilateral
Sources:
skimage/filter/rank/_crank16_bilateral.pyx
skimage/filter/rank/bilateral_cy.pyx
Executable: skivi
Module: skimage.scripts.skivi
bscript
-12
View File
@@ -1,12 +0,0 @@
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")]
check_bento_build.py
+8 -4
View File
@@ -3,6 +3,7 @@ Check that Cython extensions in setup.py files match those in bento.info.
"""
import os
import re
import sys
RE_CYTHON = re.compile("config.add_extension\(\s*['\"]([\S]+)['\"]")
@@ -62,12 +63,12 @@ def remove_common_extensions(cy_bento, cy_setup):
def print_results(cy_bento, cy_setup):
def info(text):
print
print('')
print(text)
print('-' * len(text))
if not (cy_bento or cy_setup):
print "bento.info and setup.py files match."
print("bento.info and setup.py files match.")
if cy_bento:
info("Extensions found in 'bento.info' but not in any 'setup.py:")
@@ -80,8 +81,8 @@ def print_results(cy_bento, cy_setup):
info("Consider adding the following to the 'bento.info' Library:")
for dir_path in cy_setup:
module_path = dir_path.replace('/', '.')
print BENTO_TEMPLATE.format(module_path=module_path,
dir_path=dir_path)
print(BENTO_TEMPLATE.format(module_path=module_path,
dir_path=dir_path))
if __name__ == '__main__':
@@ -93,3 +94,6 @@ if __name__ == '__main__':
cy_bento, cy_setup = remove_common_extensions(cy_bento, cy_setup)
print_results(cy_bento, cy_setup)
if cy_setup or cy_bento:
sys.exit(1)
doc/Makefile
+10 -9
View File
@@ -2,9 +2,10 @@
#
# You can set these variables from the command line.
SPHINXOPTS =
SPHINXBUILD = sphinx-build
PAPER =
PYTHON ?= python
SPHINXOPTS ?=
SPHINXBUILD ?= sphinx-build
PAPER ?=
# Internal variables.
PAPEROPT_a4 = -D latex_paper_size=a4
@@ -36,14 +37,14 @@ clean:
-find ./source/auto_examples/* -type f | grep -v blank | xargs rm -f
api:
@mkdir -p source/api
python tools/build_modref_templates.py
$(PYTHON) tools/build_modref_templates.py
@echo "Build API docs...done."
random_gallery:
@cd source && python random_gallery.py
@cd source && $(PYTHON) random_gallery.py
coveragetable:
@cd source && python coverage_generator.py
@cd source && $(PYTHON) coverage_generator.py
html: api coveragetable random_gallery
$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(DEST)/html
@@ -90,7 +91,7 @@ devhelp:
@echo "# ln -s build/devhelp $$HOME/.local/share/devhelp/scikitimage"
@echo "# devhelp"
latex:
latex: api
$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(DEST)/latex
@echo
@echo "Build finished; the LaTeX files are in $(DEST)/latex."
@@ -120,10 +121,10 @@ doctest:
"results in build/doctest/output.txt."
gh-pages:
python gh-pages.py
$(PYTHON) gh-pages.py
gitwash:
python tools/gitwash/gitwash_dumper.py source scikit-image \
$(PYTHON) tools/gitwash/gitwash_dumper.py source scikit-image \
--project-url=http://scikit-image.org \
--project-ml-url=http://groups.google.com/group/scikit-image \
--repo-name=scikit-image \
doc/examples/applications/plot_coins_segmentation.py
+11 -8
View File
@@ -3,7 +3,7 @@
Comparing edge-based segmentation and region-based segmentation
===============================================================
In this example, we will see how to segment objects from a background. We use
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.
"""
@@ -108,9 +108,8 @@ closed are not filled correctly, as is the case for one unfilled coin above.
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.
We therefore try a region-based method using the watershed transform. First, we
find an elevation map using the Sobel gradient of the image.
"""
@@ -142,7 +141,8 @@ plt.title('markers')
"""
.. image:: PLOT2RST.current_figure
Finally, we use the watershed transform to fill regions of the elevation map starting from the markers determined above:
Finally, we use the watershed transform to fill regions of the elevation map
starting from the markers determined above:
"""
segmentation = morphology.watershed(elevation_map, markers)
@@ -155,13 +155,16 @@ plt.title('segmentation')
"""
.. image:: PLOT2RST.current_figure
This last method works even better, and the coins can be segmented and
labeled individually.
This last method works even better, and the coins can be segmented and labeled
individually.
"""
from skimage.color import label2rgb
segmentation = ndimage.binary_fill_holes(segmentation - 1)
labeled_coins, _ = ndimage.label(segmentation)
image_label_overlay = label2rgb(labeled_coins, image=coins)
plt.figure(figsize=(6, 3))
plt.subplot(121)
@@ -169,7 +172,7 @@ plt.imshow(coins, cmap=plt.cm.gray, interpolation='nearest')
plt.contour(segmentation, [0.5], linewidths=1.2, colors='y')
plt.axis('off')
plt.subplot(122)
plt.imshow(labeled_coins, cmap=plt.cm.spectral, interpolation='nearest')
plt.imshow(image_label_overlay, interpolation='nearest')
plt.axis('off')
plt.subplots_adjust(**margins)
doc/examples/applications/plot_geometric.py
+5 -3
View File
@@ -7,6 +7,8 @@ In this example, we will see how to use geometric transformations in the context
of image processing.
"""
from __future__ import print_function
import math
import numpy as np
import matplotlib.pyplot as plt
@@ -31,7 +33,7 @@ First we create a transformation using explicit parameters:
tform = tf.SimilarityTransform(scale=1, rotation=math.pi / 2,
translation=(0, 1))
print tform._matrix
print(tform._matrix)
"""
Alternatively you can define a transformation by the transformation matrix
@@ -49,8 +51,8 @@ systems:
"""
coord = [1, 0]
print tform2(coord)
print tform2.inverse(tform(coord))
print(tform2(coord))
print(tform2.inverse(tform(coord)))
"""
Image warping
doc/examples/applications/plot_morphology.py
+274
View File
@@ -0,0 +1,274 @@
"""
=======================
Morphological Filtering
=======================
Morphological image processing is a collection of non-linear operations related
to the shape or morphology of features in an image, such as boundaries,
skeletons, etc. In any given technique, we probe an image with a small shape or
template called a structuring element, which defines the region of interest or
neighborhood around a pixel.
In this document we outline the following basic morphological operations:
1. Erosion
2. Dilation
3. Opening
4. Closing
5. White Tophat
6. Black Tophat
7. Skeletonize
8. Convex Hull
To get started, let's load an image using ``io.imread``. Note that morphology
functions only work on gray-scale or binary images, so we set ``as_grey=True``.
"""
import matplotlib.pyplot as plt
from skimage.data import data_dir
from skimage.util import img_as_ubyte
from skimage import io
plt.gray()
phantom = img_as_ubyte(io.imread(data_dir+'/phantom.png', as_grey=True))
plt.imshow(phantom)
"""
.. image:: PLOT2RST.current_figure
Let's also define a convenience function for plotting comparisons:
"""
def plot_comparison(original, filtered, filter_name):
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4))
ax1.imshow(original)
ax1.set_title('original')
ax1.axis('off')
ax2.imshow(filtered)
ax2.set_title(filter_name)
ax2.axis('off')
"""
Erosion
=======
Morphological ``erosion`` sets a pixel at (i, j) to the *minimum over all
pixels in the neighborhood centered at (i, j)*. The structuring element,
``selem``, passed to ``erosion`` is a boolean array that describes this
neighborhood. Below, we use ``disk`` to create a circular structuring element,
which we use for most of the following examples.
"""
from skimage.morphology import erosion, dilation, opening, closing, white_tophat
from skimage.morphology import black_tophat, skeletonize, convex_hull_image
from skimage.morphology import disk
selem = disk(6)
eroded = erosion(phantom, selem)
plot_comparison(phantom, eroded, 'erosion')
"""
.. image:: PLOT2RST.current_figure
Notice how the white boundary of the image disappears or gets eroded as we
increase the size of the disk. Also notice the increase in size of the two
black ellipses in the center and the disappearance of the 3 light grey
patches in the lower part of the image.
Dilation
========
Morphological ``dilation`` sets a pixel at (i, j) to the *maximum over all
pixels in the neighborhood centered at (i, j)*. Dilation enlarges bright
regions and shrinks dark regions.
"""
dilated = dilation(phantom, selem)
plot_comparison(phantom, dilated, 'dilation')
"""
.. image:: PLOT2RST.current_figure
Notice how the white boundary of the image thickens, or gets dilated, as we
increase the size of the disk. Also notice the decrease in size of the two
black ellipses in the centre, and the thickening of the light grey circle in
the center and the 3 patches in the lower part of the image.
Opening
=======
Morphological ``opening`` on an image is defined as an *erosion followed by a
dilation*. Opening can remove small bright spots (i.e. "salt") and connect
small dark cracks.
"""
opened = opening(phantom, selem)
plot_comparison(phantom, opened, 'opening')
"""
.. image:: PLOT2RST.current_figure
Since ``opening`` an image starts with an erosion operation, light regions that
are *smaller* than the structuring element are removed. The dilation operation
that follows ensures that light regions that are *larger* than the structuring
element retain their original size. Notice how the light and dark shapes in the
center their original thickness but the 3 lighter patches in the bottom get
completely eroded. The size dependence is highlighted by the outer white ring:
The parts of the ring thinner than the structuring element were completely
erased, while the thicker region at the top retains its original thickness.
Closing
=======
Morphological ``closing`` on an image is defined as a *dilation followed by an
erosion*. Closing can remove small dark spots (i.e. "pepper") and connect
small bright cracks.
To illustrate this more clearly, let's add a small crack to the white border:
"""
phantom = img_as_ubyte(io.imread(data_dir+'/phantom.png', as_grey=True))
phantom[10:30, 200:210] = 0
closed = closing(phantom, selem)
plot_comparison(phantom, closed, 'closing')
"""
.. image:: PLOT2RST.current_figure
Since ``closing`` an image starts with an dilation operation, dark regions
that are *smaller* than the structuring element are removed. The dilation
operation that follows ensures that dark regions that are *larger* than the
structuring element retain their original size. Notice how the white ellipses
at the bottom get connected because of dilation, but other dark region retain
their original sizes. Also notice how the crack we added is mostly removed.
White tophat
============
The ``white_tophat`` of an image is defined as the *image minus its
morphological opening*. This operation returns the bright spots of the image
that are smaller than the structuring element.
To make things interesting, we'll add bright and dark spots to the image:
"""
phantom = img_as_ubyte(io.imread(data_dir+'/phantom.png', as_grey=True))
phantom[340:350, 200:210] = 255
phantom[100:110, 200:210] = 0
w_tophat = white_tophat(phantom, selem)
plot_comparison(phantom, w_tophat, 'white tophat')
"""
.. image:: PLOT2RST.current_figure
As you can see, the 10-pixel wide white square is highlighted since it is
smaller than the structuring element. Also, the thin, white edges around most
of the ellipse are retained because they're smaller than the structuring
element, but the thicker region at the top disappears.
Black tophat
============
The ``black_tophat`` of an image is defined as its morphological **closing
minus the original image**. This operation returns the *dark spots of the
image that are smaller than the structuring element*.
"""
b_tophat = black_tophat(phantom, selem)
plot_comparison(phantom, b_tophat, 'black tophat')
"""
.. image:: PLOT2RST.current_figure
As you can see, the 10-pixel wide black square is highlighted since it is
smaller than the structuring element.
Duality
-------
As you should have noticed, many of these operations are simply the reverse
of another operation. This duality can be summarized as follows:
1. Erosion <-> Dilation
2. Opening <-> Closing
3. White tophat <-> Black tophat
Skeletonize
===========
Thinning is used to reduce each connected component in a binary image to a
*single-pixel wide skeleton*. It is important to note that this is performed
on binary images only.
"""
from skimage import img_as_bool
horse = ~img_as_bool(io.imread(data_dir+'/horse.png', as_grey=True))
sk = skeletonize(horse)
plot_comparison(horse, sk, 'skeletonize')
"""
.. image:: PLOT2RST.current_figure
As the name suggests, this technique is used to thin the image to 1-pixel wide
skeleton by applying thinning successively.
Convex hull
===========
The ``convex_hull_image`` is the *set of pixels included in the smallest
convex polygon that surround all white pixels in the input image*. Again note
that this is also performed on binary images.
"""
hull1 = convex_hull_image(horse)
plot_comparison(horse, hull1, 'convex hull')
"""
.. image:: PLOT2RST.current_figure
As the figure illustrates, ``convex_hull_image`` gives the smallest polygon
which covers the white or True completely in the image.
If we add a small grain to the image, we can see how the convex hull adapts to
enclose that grain:
"""
import numpy as np
horse2 = np.copy(horse)
horse2[45:50, 75:80] = 1
hull2 = convex_hull_image(horse2)
plot_comparison(horse2, hull2, 'convex hull')
"""
.. image:: PLOT2RST.current_figure
Additional Resources
====================
1. `MathWorks tutorial on morphological processing
<http://www.mathworks.com/help/images/morphology-fundamentals-dilation-and-erosion.html>`_
2. `Auckland university's tutorial on Morphological Image Processing
<http://www.cs.auckland.ac.nz/courses/compsci773s1c/lectures/ImageProcessing-html/topic4.htm>`_
3. http://en.wikipedia.org/wiki/Mathematical_morphology
"""
plt.show()
doc/examples/applications/plot_rank_filters.py
+246 -194
View File
@@ -3,11 +3,11 @@
Rank filters
============
Rank filters are non-linear filters using the local greylevels ordering to
Rank filters are non-linear filters using the local gray-level ordering to
compute the filtered value. This ensemble of filters share a common base: the
local grey-level histogram extraction computed on the neighborhood of a pixel
(defined by a 2D structuring element). If the filtered value is taken as the
middle value of the histogram, we get the classical median filter.
local gray-level histogram is computed on the neighborhood of a pixel (defined
by a 2-D structuring element). If the filtered value is taken as the middle
value of the histogram, we get the classical median filter.
Rank filters can be used for several purposes such as:
@@ -26,11 +26,9 @@ Rank filters can be used for several purposes such as:
Some well known filters are specific cases of rank filters [1]_ e.g.
morphological dilation, morphological erosion, median filters.
The different implementation availables in `skimage` are compared.
In this example, we will see how to filter a greylevel image using some of the
linear and non-linear filters availables in skimage. We use the `camera`
image from `skimage.data`.
In this example, we will see how to filter a gray-level image using some of the
linear and non-linear filters available in skimage. We use the `camera` image
from `skimage.data` for all comparisons.
.. [1] Pierre Soille, On morphological operators based on rank filters, Pattern
Recognition 35 (2002) 527-535.
@@ -40,18 +38,19 @@ image from `skimage.data`.
import numpy as np
import matplotlib.pyplot as plt
from skimage import img_as_ubyte
from skimage import data
ima = data.camera()
hist = np.histogram(ima, bins=np.arange(0, 256))
noisy_image = img_as_ubyte(data.camera())
hist = np.histogram(noisy_image, bins=np.arange(0, 256))
plt.figure(figsize=(8, 3))
plt.subplot(1, 2, 1)
plt.imshow(ima, cmap=plt.cm.gray, interpolation='nearest')
plt.imshow(noisy_image, interpolation='nearest')
plt.axis('off')
plt.subplot(1, 2, 2)
plt.plot(hist[1][:-1], hist[0], lw=2)
plt.title('histogram of grey values')
plt.title('Histogram of grey values')
"""
@@ -65,50 +64,56 @@ randomly set to 0. The **median** filter is applied to remove the noise.
.. note::
there are different implementations of median filter :
There are different implementations of median filter:
`skimage.filter.median_filter` and `skimage.filter.rank.median`
"""
noise = np.random.random(ima.shape)
nima = data.camera()
nima[noise > 0.99] = 255
nima[noise < 0.01] = 0
from skimage.filter.rank import median
from skimage.morphology import disk
fig = plt.figure(figsize=[10, 7])
noise = np.random.random(noisy_image.shape)
noisy_image = img_as_ubyte(data.camera())
noisy_image[noise > 0.99] = 255
noisy_image[noise < 0.01] = 0
fig = plt.figure(figsize=(10, 7))
lo = median(nima, disk(1))
hi = median(nima, disk(5))
ext = median(nima, disk(20))
plt.subplot(2, 2, 1)
plt.imshow(nima, cmap=plt.cm.gray, vmin=0, vmax=255)
plt.xlabel('noised image')
plt.imshow(noisy_image, vmin=0, vmax=255)
plt.title('Noisy image')
plt.axis('off')
plt.subplot(2, 2, 2)
plt.imshow(lo, cmap=plt.cm.gray, vmin=0, vmax=255)
plt.xlabel('median $r=1$')
plt.imshow(median(noisy_image, disk(1)), vmin=0, vmax=255)
plt.title('Median $r=1$')
plt.axis('off')
plt.subplot(2, 2, 3)
plt.imshow(hi, cmap=plt.cm.gray, vmin=0, vmax=255)
plt.xlabel('median $r=5$')
plt.imshow(median(noisy_image, disk(5)), vmin=0, vmax=255)
plt.title('Median $r=5$')
plt.axis('off')
plt.subplot(2, 2, 4)
plt.imshow(ext, cmap=plt.cm.gray, vmin=0, vmax=255)
plt.xlabel('median $r=20$')
plt.imshow(median(noisy_image, disk(20)), vmin=0, vmax=255)
plt.title('Median $r=20$')
plt.axis('off')
"""
.. image:: PLOT2RST.current_figure
The added noise is efficiently removed, as the image defaults are small (1 pixel
wide), a small filter radius is sufficient. As the radius is increasing, objects
with a bigger size are filtered as well, such as the camera tripod. The median
filter is commonly used for noise removal because borders are preserved.
The added noise is efficiently removed, as the image defaults are small (1
pixel wide), a small filter radius is sufficient. As the radius is increasing,
objects with bigger sizes are filtered as well, such as the camera tripod. The
median filter is often used for noise removal because borders are preserved and
e.g. salt and pepper noise typically does not distort the gray-level.
Image smoothing
================
The example hereunder shows how a local **mean** smoothes the camera man image.
The example hereunder shows how a local **mean** filter smooths the camera man
image.
"""
@@ -116,13 +121,17 @@ from skimage.filter.rank import mean
fig = plt.figure(figsize=[10, 7])
loc_mean = mean(nima, disk(10))
loc_mean = mean(noisy_image, disk(10))
plt.subplot(1, 2, 1)
plt.imshow(ima, cmap=plt.cm.gray, vmin=0, vmax=255)
plt.xlabel('original')
plt.imshow(noisy_image, vmin=0, vmax=255)
plt.title('Original')
plt.axis('off')
plt.subplot(1, 2, 2)
plt.imshow(loc_mean, cmap=plt.cm.gray, vmin=0, vmax=255)
plt.xlabel('local mean $r=10$')
plt.imshow(loc_mean, vmin=0, vmax=255)
plt.title('Local mean $r=10$')
plt.axis('off')
"""
@@ -130,35 +139,41 @@ plt.xlabel('local mean $r=10$')
One may be interested in smoothing an image while preserving important borders
(median filters already achieved this), here we use the **bilateral** filter
that restricts the local neighborhood to pixel having a greylevel similar to
that restricts the local neighborhood to pixel having a gray-level similar to
the central one.
.. note::
a different implementation is available for color images in
A different implementation is available for color images in
`skimage.filter.denoise_bilateral`.
"""
from skimage.filter.rank import bilateral_mean
ima = data.camera()
selem = disk(10)
noisy_image = img_as_ubyte(data.camera())
bilat = bilateral_mean(ima.astype(np.uint16), disk(20), s0=10, s1=10)
bilat = bilateral_mean(noisy_image.astype(np.uint16), disk(20), s0=10, s1=10)
# display results
fig = plt.figure(figsize=[10, 7])
plt.subplot(2, 2, 1)
plt.imshow(ima, cmap=plt.cm.gray)
plt.xlabel('original')
plt.imshow(noisy_image, cmap=plt.cm.gray)
plt.title('Original')
plt.axis('off')
plt.subplot(2, 2, 3)
plt.imshow(bilat, cmap=plt.cm.gray)
plt.xlabel('bilateral mean')
plt.title('Bilateral mean')
plt.axis('off')
plt.subplot(2, 2, 2)
plt.imshow(ima[200:350, 350:450], cmap=plt.cm.gray)
plt.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)
plt.axis('off')
plt.subplot(2, 2, 4)
plt.imshow(bilat[200:350, 350:450], cmap=plt.cm.gray)
plt.axis('off')
"""
@@ -175,7 +190,7 @@ We compare here how the global histogram equalization is applied locally.
The equalized image [2]_ has a roughly linear cumulative distribution function
for each pixel neighborhood. The local version [3]_ of the histogram
equalization emphasizes every local greylevel variations.
equalization emphasizes every local gray-level variations.
.. [2] http://en.wikipedia.org/wiki/Histogram_equalization
.. [3] http://en.wikipedia.org/wiki/Adaptive_histogram_equalization
@@ -185,101 +200,112 @@ equalization emphasizes every local greylevel variations.
from skimage import exposure
from skimage.filter import rank
ima = data.camera()
noisy_image = img_as_ubyte(data.camera())
# equalize globally and locally
glob = exposure.equalize(ima) * 255
loc = rank.equalize(ima, disk(20))
glob = exposure.equalize(noisy_image) * 255
loc = rank.equalize(noisy_image, disk(20))
# extract histogram for each image
hist = np.histogram(ima, bins=np.arange(0, 256))
hist = np.histogram(noisy_image, bins=np.arange(0, 256))
glob_hist = np.histogram(glob, bins=np.arange(0, 256))
loc_hist = np.histogram(loc, bins=np.arange(0, 256))
plt.figure(figsize=(10, 10))
plt.subplot(321)
plt.imshow(ima, cmap=plt.cm.gray, interpolation='nearest')
plt.imshow(noisy_image, interpolation='nearest')
plt.axis('off')
plt.subplot(322)
plt.plot(hist[1][:-1], hist[0], lw=2)
plt.title('histogram of grey values')
plt.title('Histogram of gray values')
plt.subplot(323)
plt.imshow(glob, cmap=plt.cm.gray, interpolation='nearest')
plt.imshow(glob, interpolation='nearest')
plt.axis('off')
plt.subplot(324)
plt.plot(glob_hist[1][:-1], glob_hist[0], lw=2)
plt.title('histogram of grey values')
plt.title('Histogram of gray values')
plt.subplot(325)
plt.imshow(loc, cmap=plt.cm.gray, interpolation='nearest')
plt.imshow(loc, interpolation='nearest')
plt.axis('off')
plt.subplot(326)
plt.plot(loc_hist[1][:-1], loc_hist[0], lw=2)
plt.title('histogram of grey values')
plt.title('Histogram of gray values')
"""
.. image:: PLOT2RST.current_figure
another way to maximize the number of greylevels used for an image is to apply
a local autoleveling, i.e. here a pixel greylevel is proportionally remapped
between local minimum and local maximum.
Another way to maximize the number of gray-levels used for an image is to apply
a local auto-leveling, i.e. the gray-value of a pixel is proportionally
remapped between local minimum and local maximum.
The following example shows how local autolevel enhances the camara man picture.
The following example shows how local auto-level enhances the camara man
picture.
"""
from skimage.filter.rank import autolevel
ima = data.camera()
selem = disk(10)
noisy_image = img_as_ubyte(data.camera())
auto = autolevel(ima.astype(np.uint16), disk(20))
auto = autolevel(noisy_image.astype(np.uint16), disk(20))
# display results
fig = plt.figure(figsize=[10, 7])
plt.subplot(1, 2, 1)
plt.imshow(ima, cmap=plt.cm.gray)
plt.xlabel('original')
plt.imshow(noisy_image, cmap=plt.cm.gray)
plt.title('Original')
plt.axis('off')
plt.subplot(1, 2, 2)
plt.imshow(auto, cmap=plt.cm.gray)
plt.xlabel('local autolevel')
plt.title('Local autolevel')
plt.axis('off')
"""
.. image:: PLOT2RST.current_figure
This filter is very sensitive to local outlayers, see the little white spot in
the sky left part. This is due to a local maximum which is very high comparing
to the rest of the neighborhood. One can moderate this using the percentile
version of the autolevel filter which uses given percentiles (one inferior,
one superior) in place of local minimum and maximum. The example below
illustrates how the percentile parameters influence the local autolevel result.
This filter is very sensitive to local outliers, see the little white spot in
the left part of the sky. This is due to a local maximum which is very high
comparing to the rest of the neighborhood. One can moderate this using the
percentile version of the auto-level filter which uses given percentiles (one
inferior, one superior) in place of local minimum and maximum. The example
below illustrates how the percentile parameters influence the local auto-level
result.
"""
from skimage.filter.rank import percentile_autolevel
from skimage.filter.rank import autolevel_percentile
image = data.camera()
selem = disk(20)
loc_autolevel = autolevel(image, selem=selem)
loc_perc_autolevel0 = percentile_autolevel(image, selem=selem, p0=.00, p1=1.0)
loc_perc_autolevel1 = percentile_autolevel(image, selem=selem, p0=.01, p1=.99)
loc_perc_autolevel2 = percentile_autolevel(image, selem=selem, p0=.05, p1=.95)
loc_perc_autolevel3 = percentile_autolevel(image, selem=selem, p0=.1, p1=.9)
loc_perc_autolevel0 = autolevel_percentile(image, selem=selem, p0=.00, p1=1.0)
loc_perc_autolevel1 = autolevel_percentile(image, selem=selem, p0=.01, p1=.99)
loc_perc_autolevel2 = autolevel_percentile(image, selem=selem, p0=.05, p1=.95)
loc_perc_autolevel3 = autolevel_percentile(image, selem=selem, p0=.1, p1=.9)
fig, axes = plt.subplots(nrows=3, figsize=(7, 8))
ax0, ax1, ax2 = axes
plt.gray()
ax0.imshow(np.hstack((image, loc_autolevel)))
ax0.set_title('original / autolevel')
ax0.set_title('Original / auto-level')
ax1.imshow(
np.hstack((loc_perc_autolevel0, loc_perc_autolevel1)), vmin=0, vmax=255)
ax1.set_title('percentile autolevel 0%,1%')
ax1.set_title('Percentile auto-level 0%,1%')
ax2.imshow(
np.hstack((loc_perc_autolevel2, loc_perc_autolevel3)), vmin=0, vmax=255)
ax2.set_title('percentile autolevel 5% and 10%')
ax2.set_title('Percentile auto-level 5% and 10%')
for ax in axes:
ax.axis('off')
@@ -289,29 +315,35 @@ for ax in axes:
.. image:: PLOT2RST.current_figure
The morphological contrast enhancement filter replaces the central pixel by the
local maximum if the original pixel value is closest to local maximum, otherwise
by the minimum local.
local maximum if the original pixel value is closest to local maximum,
otherwise by the minimum local.
"""
from skimage.filter.rank import morph_contr_enh
from skimage.filter.rank import enhance_contrast
ima = data.camera()
noisy_image = img_as_ubyte(data.camera())
enh = morph_contr_enh(ima, disk(5))
enh = enhance_contrast(noisy_image, disk(5))
# display results
fig = plt.figure(figsize=[10, 7])
plt.subplot(2, 2, 1)
plt.imshow(ima, cmap=plt.cm.gray)
plt.xlabel('original')
plt.imshow(noisy_image, cmap=plt.cm.gray)
plt.title('Original')
plt.axis('off')
plt.subplot(2, 2, 3)
plt.imshow(enh, cmap=plt.cm.gray)
plt.xlabel('local morphlogical contrast enhancement')
plt.title('Local morphological contrast enhancement')
plt.axis('off')
plt.subplot(2, 2, 2)
plt.imshow(ima[200:350, 350:450], cmap=plt.cm.gray)
plt.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)
plt.axis('off')
plt.subplot(2, 2, 4)
plt.imshow(enh[200:350, 350:450], cmap=plt.cm.gray)
plt.axis('off')
"""
@@ -322,24 +354,30 @@ percentile *p0* and *p1* instead of the local minimum and maximum.
"""
from skimage.filter.rank import percentile_morph_contr_enh
from skimage.filter.rank import enhance_contrast_percentile
ima = data.camera()
noisy_image = img_as_ubyte(data.camera())
penh = percentile_morph_contr_enh(ima, disk(5), p0=.1, p1=.9)
penh = enhance_contrast_percentile(noisy_image, disk(5), p0=.1, p1=.9)
# display results
fig = plt.figure(figsize=[10, 7])
plt.subplot(2, 2, 1)
plt.imshow(ima, cmap=plt.cm.gray)
plt.xlabel('original')
plt.imshow(noisy_image, cmap=plt.cm.gray)
plt.title('Original')
plt.axis('off')
plt.subplot(2, 2, 3)
plt.imshow(penh, cmap=plt.cm.gray)
plt.xlabel('local percentile morphlogical\n contrast enhancement')
plt.title('Local percentile morphological\n contrast enhancement')
plt.axis('off')
plt.subplot(2, 2, 2)
plt.imshow(ima[200:350, 350:450], cmap=plt.cm.gray)
plt.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)
plt.axis('off')
plt.subplot(2, 2, 4)
plt.imshow(penh[200:350, 350:450], cmap=plt.cm.gray)
plt.axis('off')
"""
@@ -348,20 +386,20 @@ plt.imshow(penh[200:350, 350:450], cmap=plt.cm.gray)
Image threshold
===============
The Otsu's threshold [1]_ method can be applied locally using the local
greylevel distribution. In the example below, for each pixel, an "optimal"
threshold is determined by maximizing the variance between two classes of pixels
of the local neighborhood defined by a structuring element.
The Otsu threshold [1]_ method can be applied locally using the local gray-
level distribution. In the example below, for each pixel, an "optimal"
threshold is determined by maximizing the variance between two classes of
pixels of the local neighborhood defined by a structuring element.
The example compares the local threshold with the global threshold
`skimage.filter.threshold_otsu`.
.. note::
Local thresholding is much slower than global one. There exists a function
for global Otsu thresholding: `skimage.filter.threshold_otsu`.
Local is much slower than global thresholding. A function for global Otsu
thresholding can be found in : `skimage.filter.threshold_otsu`.
.. [1] http://en.wikipedia.org/wiki/Otsu's_method
.. [4] http://en.wikipedia.org/wiki/Otsu's_method
"""
@@ -382,27 +420,35 @@ t_glob_otsu = threshold_otsu(p8)
glob_otsu = p8 >= t_glob_otsu
plt.figure()
plt.subplot(2, 2, 1)
plt.imshow(p8, cmap=plt.cm.gray)
plt.xlabel('original')
plt.title('Original')
plt.colorbar()
plt.axis('off')
plt.subplot(2, 2, 2)
plt.imshow(t_loc_otsu, cmap=plt.cm.gray)
plt.xlabel('local Otsu ($radius=%d$)' % radius)
plt.title('Local Otsu ($r=%d$)' % radius)
plt.colorbar()
plt.axis('off')
plt.subplot(2, 2, 3)
plt.imshow(p8 >= t_loc_otsu, cmap=plt.cm.gray)
plt.xlabel('original>=local Otsu' % t_glob_otsu)
plt.title('Original >= local Otsu' % t_glob_otsu)
plt.axis('off')
plt.subplot(2, 2, 4)
plt.imshow(glob_otsu, cmap=plt.cm.gray)
plt.xlabel('global Otsu ($t=%d$)' % t_glob_otsu)
plt.title('Global Otsu ($t=%d$)' % t_glob_otsu)
plt.axis('off')
"""
.. image:: PLOT2RST.current_figure
The following example shows how local Otsu's threshold handles a global level
shift applied to a synthetic image .
The following example shows how local Otsu thresholding handles a global level
shift applied to a synthetic image.
"""
@@ -413,13 +459,18 @@ m = (np.tile(x, (n, 1)) * np.linspace(0.1, 1, n) * 128 + 128).astype(np.uint8)
radius = 10
t = rank.otsu(m, disk(radius))
plt.figure()
plt.subplot(1, 2, 1)
plt.imshow(m)
plt.xlabel('original')
plt.title('Original')
plt.axis('off')
plt.subplot(1, 2, 2)
plt.imshow(m >= t, interpolation='nearest')
plt.xlabel('local Otsu ($radius=%d$)' % radius)
plt.title('Local Otsu ($r=%d$)' % radius)
plt.axis('off')
"""
@@ -428,7 +479,7 @@ plt.xlabel('local Otsu ($radius=%d$)' % radius)
Image morphology
================
Local maximum and local minimum are the base operators for greylevel
Local maximum and local minimum are the base operators for gray-level
morphology.
.. note::
@@ -436,33 +487,41 @@ morphology.
`skimage.dilate` and `skimage.erode` are equivalent filters (see below for
comparison).
Here is an example of the classical morphological greylevel filters: opening,
Here is an example of the classical morphological gray-level filters: opening,
closing and morphological gradient.
"""
from skimage.filter.rank import maximum, minimum, gradient
ima = data.camera()
noisy_image = img_as_ubyte(data.camera())
closing = maximum(minimum(ima, disk(5)), disk(5))
opening = minimum(maximum(ima, disk(5)), disk(5))
grad = gradient(ima, disk(5))
closing = maximum(minimum(noisy_image, disk(5)), disk(5))
opening = minimum(maximum(noisy_image, disk(5)), disk(5))
grad = gradient(noisy_image, disk(5))
# display results
fig = plt.figure(figsize=[10, 7])
plt.subplot(2, 2, 1)
plt.imshow(ima, cmap=plt.cm.gray)
plt.xlabel('original')
plt.imshow(noisy_image, cmap=plt.cm.gray)
plt.title('Original')
plt.axis('off')
plt.subplot(2, 2, 2)
plt.imshow(closing, cmap=plt.cm.gray)
plt.xlabel('greylevel closing')
plt.title('Gray-level closing')
plt.axis('off')
plt.subplot(2, 2, 3)
plt.imshow(opening, cmap=plt.cm.gray)
plt.xlabel('greylevel opening')
plt.title('Gray-level opening')
plt.axis('off')
plt.subplot(2, 2, 4)
plt.imshow(grad, cmap=plt.cm.gray)
plt.xlabel('morphological gradient')
plt.title('Morphological gradient')
plt.axis('off')
"""
@@ -471,13 +530,14 @@ plt.xlabel('morphological gradient')
Feature extraction
===================
Local histogram can be exploited to compute local entropy, which is related to
Local histograms can be exploited to compute local entropy, which is related to
the local image complexity. Entropy is computed using base 2 logarithm i.e. the
filter returns the minimum number of bits needed to encode local greylevel
filter returns the minimum number of bits needed to encode local gray-level
distribution.
`skimage.rank.entropy` returns local entropy on a given structuring element.
The following example shows this filter applied on 8- and 16- bit images.
`skimage.rank.entropy` returns the local entropy on a given structuring
element. The following example shows applies this filter on 8- and 16-bit
images.
.. note::
@@ -492,47 +552,36 @@ from skimage.morphology import disk
import numpy as np
import matplotlib.pyplot as plt
# defining a 8- and a 16-bit test images
a8 = data.camera()
a16 = data.camera().astype(np.uint16) * 4
image = data.camera()
ent8 = entropy(a8, disk(5)) # pixel value contain 10x the local entropy
ent16 = entropy(a16, disk(5)) # pixel value contain 1000x the local entropy
plt.figure(figsize=(10, 4))
# display results
plt.figure(figsize=(10, 10))
plt.subplot(2, 2, 1)
plt.imshow(a8, cmap=plt.cm.gray)
plt.xlabel('8-bit image')
plt.subplot(1, 2, 1)
plt.imshow(image, cmap=plt.cm.gray)
plt.title('Image')
plt.colorbar()
plt.axis('off')
plt.subplot(2, 2, 2)
plt.imshow(ent8, cmap=plt.cm.jet)
plt.xlabel('entropy*10')
plt.colorbar()
plt.subplot(2, 2, 3)
plt.imshow(a16, cmap=plt.cm.gray)
plt.xlabel('16-bit image')
plt.colorbar()
plt.subplot(2, 2, 4)
plt.imshow(ent16, cmap=plt.cm.jet)
plt.xlabel('entropy*1000')
plt.subplot(1, 2, 2)
plt.imshow(entropy(image, disk(5)), cmap=plt.cm.jet)
plt.title('Entropy')
plt.colorbar()
plt.axis('off')
"""
.. image:: PLOT2RST.current_figure
Implementation
================
==============
The central part of the `skimage.rank` filters is build on a sliding window that
update local greylevel histogram. This approach limits the algorithm complexity
to O(n) where n is the number of image pixels. The complexity is also limited
with respect to the structuring element size.
The central part of the `skimage.rank` filters is build on a sliding window
that updates the local gray-level histogram. This approach limits the algorithm
complexity to O(n) where n is the number of image pixels. The complexity is
also limited with respect to the structuring element size.
In the following we compare the performance of different implementations
available in `skimage`.
"""
@@ -583,10 +632,10 @@ def ndi_med(image, n):
Comparison between
* `rank.maximum`
* `cmorph.dilate`
* `filter.rank.maximum`
* `morphology.dilate`
on increasing structuring element size
on increasing structuring element size:
"""
@@ -603,18 +652,18 @@ for r in e_range:
rec = np.asarray(rec)
plt.figure()
plt.title('increasing element size')
plt.ylabel('time (ms)')
plt.xlabel('element radius')
plt.title('Performance with respect to element size')
plt.ylabel('Time (ms)')
plt.title('Element radius')
plt.plot(e_range, rec)
plt.legend(['crank.maximum', 'cmorph.dilate'])
plt.legend(['filter.rank.maximum', 'morphology.dilate'])
"""
and increasing image size
.. image:: PLOT2RST.current_figure
and increasing image size:
"""
r = 9
@@ -623,7 +672,7 @@ elem = disk(r + 1)
rec = []
s_range = range(100, 1000, 100)
for s in s_range:
a = (np.random.random((s, s)) * 256).astype('uint8')
a = (np.random.random((s, s)) * 256).astype(np.uint8)
(rc, ms_rc) = cr_max(a, elem)
(rcm, ms_rcm) = cm_dil(a, elem)
rec.append((ms_rc, ms_rcm))
@@ -631,11 +680,11 @@ for s in s_range:
rec = np.asarray(rec)
plt.figure()
plt.title('increasing image size')
plt.ylabel('time (ms)')
plt.xlabel('image size')
plt.title('Performance with respect to image size')
plt.ylabel('Time (ms)')
plt.title('Image size')
plt.plot(s_range, rec)
plt.legend(['crank.maximum', 'cmorph.dilate'])
plt.legend(['filter.rank.maximum', 'morphology.dilate'])
"""
@@ -644,11 +693,11 @@ plt.legend(['crank.maximum', 'cmorph.dilate'])
Comparison between:
* `rank.median`
* `ctmf.median_filter`
* `ndimage.percentile`
* `filter.rank.median`
* `filter.median_filter`
* `scipy.ndimage.percentile`
on increasing structuring element size
on increasing structuring element size:
"""
@@ -666,27 +715,29 @@ for r in e_range:
rec = np.asarray(rec)
plt.figure()
plt.title('increasing element size')
plt.title('Performance with respect to element size')
plt.plot(e_range, rec)
plt.legend(['rank.median', 'ctmf.median_filter', 'ndimage.percentile'])
plt.ylabel('time (ms)')
plt.xlabel('element radius')
plt.legend(['filter.rank.median', 'filter.median_filter',
'scipy.ndimage.percentile'])
plt.ylabel('Time (ms)')
plt.title('Element radius')
"""
.. image:: PLOT2RST.current_figure
comparison of outcome of the three methods
Comparison of outcome of the three methods:
"""
plt.figure()
plt.imshow(np.hstack((rc, rctmf, rndi)))
plt.xlabel('rank.median vs ctmf.median_filter vs ndimage.percentile')
plt.title('filter.rank.median vs filtermedian_filter vs scipy.ndimage.percentile')
plt.axis('off')
"""
.. image:: PLOT2RST.current_figure
and increasing image size
and increasing image size:
"""
@@ -696,7 +747,7 @@ elem = disk(r + 1)
rec = []
s_range = [100, 200, 500, 1000]
for s in s_range:
a = (np.random.random((s, s)) * 256).astype('uint8')
a = (np.random.random((s, s)) * 256).astype(np.uint8)
(rc, ms_rc) = cr_med(a, elem)
rctmf, ms_rctmf = ctmf_med(a, r)
rndi, ms_ndi = ndi_med(a, r)
@@ -705,11 +756,12 @@ for s in s_range:
rec = np.asarray(rec)
plt.figure()
plt.title('increasing image size')
plt.title('Performance with respect to image size')
plt.plot(s_range, rec)
plt.legend(['rank.median', 'ctmf.median_filter', 'ndimage.percentile'])
plt.ylabel('time (ms)')
plt.xlabel('image size')
plt.legend(['filter.rank.median', 'filter.median_filter',
'scipy.ndimage.percentile'])
plt.ylabel('Time (ms)')
plt.title('Image size')
"""
.. image:: PLOT2RST.current_figure
doc/examples/plot_16bitbilateral.py
-47
View File
@@ -1,47 +0,0 @@
"""
==============================
Bilateral mean
==============================
This example compares
* local mean
* percentile mean
* bilateral mean
build on the local histogram distribution
local mean uses all pixels belonging to the structuring element to compute average gray level,
percentile mean uses only values between percentiles p0 and p1 (here 10% and 90%),
whereas bilateral mean uses only pixels of the structuring element having a gray level situated inside
g-s0 and g+s1 (here g-500 and g+500).
The filters are applied on a 16 bit image (actual bitdepth is 12bit).
Percentile and usual mean give here similar results, these filters smooth the complete image (background and details).
Bilateral mean exhibits a high filtering rate for continuous area (i.e. background) while image higher frequencies
remains untouched.
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage import data
from skimage.morphology import disk
import skimage.filter.rank as rank
a16 = (data.coins()).astype('uint16') * 16
selem = disk(20)
f1 = rank.percentile_mean(a16, selem=selem, p0=.1, p1=.9)
f2 = rank.bilateral_mean(a16, selem=selem, s0=500, s1=500)
f3 = rank.mean(a16, selem=selem)
# display results
fig, axes = plt.subplots(nrows=3, figsize=(15, 10))
ax0, ax1, ax2 = axes
ax0.imshow(np.hstack((a16, f1)))
ax0.set_title('percentile mean')
ax1.imshow(np.hstack((a16, f2)))
ax1.set_title('bilateral mean')
ax2.imshow(np.hstack((a16, f3)))
ax2.set_title('local mean')
plt.show()
doc/examples/plot_canny.py
+3 -1
View File
@@ -13,12 +13,15 @@ thresholding on the gradient magnitude.
The Canny has three adjustable parameters: the width of the Gaussian (the
noisier the image, the greater the width), and the low and high threshold for
the hysteresis thresholding.
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage
from skimage import filter
# Generate noisy image of a square
im = np.zeros((128, 128))
im[32:-32, 32:-32] = 1
@@ -52,6 +55,5 @@ plt.title('Canny filter, $\sigma=3$', fontsize=20)
plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9,
bottom=0.02, left=0.02, right=0.98)
plt.show()
doc/examples/plot_censure_keypoints.py
+54
View File
@@ -0,0 +1,54 @@
"""
=========================
CenSurE Feature Detection
=========================
In this example we detect and plot the CenSurE (Center Surround Extrema)
features at various scales using Difference of Boxes, Octagon and Star shaped
bi-level filters.
"""
from skimage.feature import keypoints_censure
from skimage.data import lena
from skimage.color import rgb2gray
import matplotlib.pyplot as plt
# Initializing the parameters for Censure keypoints
img = lena()
gray_img = rgb2gray(img)
min_scale = 2
max_scale = 6
non_max_threshold = 0.15
line_threshold = 10
_, ax = plt.subplots(nrows=(max_scale - min_scale - 1), ncols=3,
figsize=(6, 6))
plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.94,
bottom=0.02, left=0.06, right=0.98)
# Detecting Censure keypoints for the following filters
for col, mode in enumerate(['dob', 'octagon', 'star']):
ax[0, col].set_title(mode.upper(), fontsize=12)
keypoints, scales = keypoints_censure(gray_img, min_scale, max_scale,
mode, non_max_threshold,
line_threshold)
# Plotting Censure features at all the scales
for row, scale in enumerate(range(min_scale + 1, max_scale)):
mask = scales == scale
x = keypoints[mask, 1]
y = keypoints[mask, 0]
s = 0.5 * 2 ** (scale + min_scale + 1)
ax[row, col].imshow(img)
ax[row, col].scatter(x, y, s, facecolors='none', edgecolors='b')
ax[row, col].set_xticks([])
ax[row, col].set_yticks([])
ax[row, col].axis((0, img.shape[1], img.shape[0], 0))
if col == 0:
ax[row, col].set_ylabel('Scale %d' % scale, fontsize=12)
plt.show()
doc/examples/plot_circular_elliptical_hough_transform.py
Executable
+147
View File
@@ -0,0 +1,147 @@
"""
========================================
Circular and Elliptical Hough Transforms
========================================
The Hough transform in its simplest form is a `method to detect
straight lines <http://en.wikipedia.org/wiki/Hough_transform>`__
but it can also be used to detect circles or ellipses.
The algorithm assumes that the edge is detected and it is robust against
noise or missing points.
Circle detection
================
In the following example, the Hough transform is used to detect
coin positions and match their edges. We provide a range of
plausible radii. For each radius, two circles are extracted and
we finally keep the five most prominent candidates.
The result shows that coin positions are well-detected.
Algorithm overview
------------------
Given a black circle on a white background, we first guess its
radius (or a range of radii) to construct a new circle.
This circle is applied on each black pixel of the original picture
and the coordinates of this circle are voting in an accumulator.
From this geometrical construction, the original circle center
position receives the highest score.
Note that the accumulator size is built to be larger than the
original picture in order to detect centers outside the frame.
Its size is extended by two times the larger radius.
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage import data, filter, color
from skimage.transform import hough_circle
from skimage.feature import peak_local_max
from skimage.draw import circle_perimeter
from skimage.util import img_as_ubyte
# Load picture and detect edges
image = img_as_ubyte(data.coins()[0:95, 70:370])
edges = filter.canny(image, sigma=3, low_threshold=10, high_threshold=50)
fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(6, 6))
# Detect two radii
hough_radii = np.arange(15, 30, 2)
hough_res = hough_circle(edges, hough_radii)
centers = []
accums = []
radii = []
for radius, h in zip(hough_radii, hough_res):
# For each radius, extract two circles
peaks = peak_local_max(h, num_peaks=2)
centers.extend(peaks)
accums.extend(h[peaks[:, 0], peaks[:, 1]])
radii.extend([radius, radius])
# Draw the most prominent 5 circles
image = color.gray2rgb(image)
for idx in np.argsort(accums)[::-1][:5]:
center_x, center_y = centers[idx]
radius = radii[idx]
cx, cy = circle_perimeter(center_y, center_x, radius)
image[cy, cx] = (220, 20, 20)
ax.imshow(image, cmap=plt.cm.gray)
plt.show()
"""
Ellipse detection
=================
In this second example, the aim is to detect the edge of a coffee cup.
Basically, this is a projection of a circle, i.e. an ellipse.
The problem to solve is much more difficult because five parameters have to be
determined, instead of three for circles.
Algorithm overview
------------------
The algorithm takes two different points belonging to the ellipse. It assumes
that it is the main axis. A loop on all the other points determines how much
an ellipse passes to them. A good match corresponds to high accumulator values.
A full description of the algorithm can be found in reference [1]_.
References
----------
.. [1] Xie, Yonghong, and Qiang Ji. "A new efficient ellipse detection
method." Pattern Recognition, 2002. Proceedings. 16th International
Conference on. Vol. 2. IEEE, 2002
"""
import matplotlib.pyplot as plt
from skimage import data, filter, color
from skimage.transform import hough_ellipse
from skimage.draw import ellipse_perimeter
# Load picture, convert to grayscale and detect edges
image_rgb = data.load('coffee.png')[0:220, 100:450]
image_gray = color.rgb2gray(image_rgb)
edges = filter.canny(image_gray, sigma=2.0,
low_threshold=0.55, high_threshold=0.8)
# Perform a Hough Transform
# The accuracy corresponds to the bin size of a major axis.
# The value is chosen in order to get a single high accumulator.
# The threshold eliminates low accumulators
accum = hough_ellipse(edges, accuracy=10, threshold=170, min_size=50)
accum.sort(key=lambda x:x[5])
# Estimated parameters for the ellipse
center_y = int(accum[-1][0])
center_x = int(accum[-1][1])
xradius = int(accum[-1][2])
yradius = int(accum[-1][3])
angle = np.pi - accum[-1][4]
# Draw the ellipse on the original image
cx, cy = ellipse_perimeter(center_y, center_x,
yradius, xradius, orientation=angle)
image_rgb[cy, cx] = (0, 0, 1)
# Draw the edge (white) and the resulting ellipse (red)
edges = color.gray2rgb(edges)
edges[cy, cx] = (250, 0, 0)
fig = plt.subplots(figsize=(10, 6))
plt.subplot(1, 2, 1)
plt.title('Original picture')
plt.imshow(image_rgb)
plt.subplot(1, 2, 2)
plt.title('Edge (white) and result (red)')
plt.imshow(edges)
plt.show()
doc/examples/plot_circular_hough_transform.py
-72
View File
@@ -1,72 +0,0 @@
"""
========================
Circular Hough Transform
========================
The Hough transform in its simplest form is a `method to detect
straight lines <http://en.wikipedia.org/wiki/Hough_transform>`__
but it can also be used to detect circles.
In the following example, the Hough transform is used to detect
coin positions and match their edges. We provide a range of
plausible radii. For each radius, two circles are extracted and
we finally keep the five most prominent candidates.
The result shows that coin positions are well-detected.
Algorithm overview
------------------
Given a black circle on a white background, we first guess its
radius (or a range of radii) to construct a new circle.
This circle is applied on each black pixel of the original picture
and the coordinates of this circle are voting in an accumulator.
From this geometrical construction, the original circle center
position receives the highest score.
Note that the accumulator size is built to be larger than the
original picture in order to detect centers outside the frame.
Its size is extended by two times the larger radius.
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage import data, filter, color
from skimage.transform import hough_circle
from skimage.feature import peak_local_max
from skimage.draw import circle_perimeter
# Load picture and detect edges
image = data.coins()[0:95, 70:370]
edges = filter.canny(image, sigma=3, low_threshold=10, high_threshold=50)
fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(6, 6))
# Detect two radii
hough_radii = np.arange(15, 30, 2)
hough_res = hough_circle(edges, hough_radii)
centers = []
accums = []
radii = []
for radius, h in zip(hough_radii, hough_res):
# For each radius, extract two circles
peaks = peak_local_max(h, num_peaks=2)
centers.extend(peaks - hough_radii.max())
accums.extend(h[peaks[:, 0], peaks[:, 1]])
radii.extend([radius, radius])
# Draw the most prominent 5 circles
image = color.gray2rgb(image)
for idx in np.argsort(accums)[::-1][:5]:
center_x, center_y = centers[idx]
radius = radii[idx]
cx, cy = circle_perimeter(center_y, center_x, radius)
image[cy, cx] = (220, 20, 20)
ax.imshow(image, cmap=plt.cm.gray)
plt.show()
doc/examples/plot_contours.py
+3 -5
View File
@@ -15,13 +15,12 @@ Cubes: A High Resolution 3D Surface Construction Algorithm. Computer Graphics
(SIGGRAPH 87 Proceedings) 21(4) July 1987, p. 163-170).
"""
from skimage import data
from skimage import measure
import numpy as np
import matplotlib.pyplot as plt
from skimage import measure
# Construct some test data
x, y = np.ogrid[-np.pi:np.pi:100j, -np.pi:np.pi:100j]
r = np.sin(np.exp((np.sin(x)**3 + np.cos(y)**2)))
@@ -39,4 +38,3 @@ plt.axis('image')
plt.xticks([])
plt.yticks([])
plt.show()
doc/examples/plot_convex_hull.py
+19 -4
View File
@@ -13,12 +13,12 @@ A good overview of the algorithm is given on `Steve Eddin's blog
<http://blogs.mathworks.com/steve/2011/10/04/binary-image-convex-hull-algorithm-notes/>`__.
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage.morphology import convex_hull_image
image = np.array(
[[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0, 0],
@@ -27,9 +27,24 @@ image = np.array(
[0, 1, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=float)
chull = convex_hull_image(image)
image[chull] += 1.7
image -= -1.7
original_image = np.copy(image)
chull = convex_hull_image(image)
image[chull] += 1
# image is now:
#[[ 0. 0. 0. 0. 0. 0. 0. 0. 0.]
# [ 0. 0. 0. 0. 2. 0. 0. 0. 0.]
# [ 0. 0. 0. 2. 1. 2. 0. 0. 0.]
# [ 0. 0. 2. 1. 1. 1. 2. 0. 0.]
# [ 0. 2. 1. 1. 1. 1. 1. 2. 0.]
# [ 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
fig = plt.subplots(figsize=(10, 6))
plt.subplot(1, 2, 1)
plt.title('Original picture')
plt.imshow(original_image, cmap=plt.cm.gray, interpolation='nearest')
plt.subplot(1, 2, 2)
plt.title('Transformed picture')
plt.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
plt.show()
doc/examples/plot_corner.py
+1 -1
View File
@@ -10,7 +10,6 @@ position of corners.
.. [2] http://en.wikipedia.org/wiki/Interest_point_detection
"""
from matplotlib import pyplot as plt
from skimage import data
@@ -18,6 +17,7 @@ from skimage.feature import corner_harris, corner_subpix, corner_peaks
from skimage.transform import warp, AffineTransform
from skimage.draw import ellipse
tform = AffineTransform(scale=(1.3, 1.1), rotation=1, shear=0.7,
translation=(210, 50))
image = warp(data.checkerboard(), tform.inverse, output_shape=(350, 350))
doc/examples/plot_daisy.py
-1
View File
@@ -11,7 +11,6 @@ representations.
In this example a limited number of DAISY descriptors are extracted at a large
scale for illustrative purposes.
"""
from skimage.feature import daisy
from skimage import data
import matplotlib.pyplot as plt
doc/examples/plot_denoise.py
+2 -2
View File
@@ -25,13 +25,13 @@ A bilateral filter is an edge-preserving and noise reducing filter. It averages
pixels based on their spatial closeness and radiometric similarity.
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage import data, color, img_as_float
from skimage import data, img_as_float
from skimage.filter import denoise_tv_chambolle, denoise_bilateral
lena = img_as_float(data.lena())
lena = lena[220:300, 220:320]
doc/examples/plot_edge_filter.py
+31
View File
@@ -0,0 +1,31 @@
"""
==============
Edge operators
==============
Edge operators are used in image processing within edge detection algorithms.
They are discrete differentiation operators, computing an approximation of the
gradient of the image intensity function.
"""
import matplotlib.pyplot as plt
from skimage.data import camera
from skimage.filter import roberts, sobel
image = camera()
edge_roberts = roberts(image)
edge_sobel = sobel(image)
fig, (ax0, ax1) = plt.subplots(ncols=2)
ax0.imshow(edge_roberts, cmap=plt.cm.gray)
ax0.set_title('Roberts Edge Detection')
ax0.axis('off')
ax1.imshow(edge_sobel, cmap=plt.cm.gray)
ax1.set_title('Sobel Edge Detection')
ax1.axis('off')
plt.show()
doc/examples/plot_entropy.py
+17 -29
View File
@@ -1,44 +1,32 @@
"""
===================
=======
Entropy
===================
=======
Image entropy is a quantity which is used to describe the amount of information
coded in an image.
"""
import matplotlib.pyplot as plt
from skimage import data
from skimage.filter.rank import entropy
from skimage.morphology import disk
import numpy as np
import matplotlib.pyplot as plt
from skimage.util import img_as_ubyte
# defining a 8- and a 16-bit test images
a8 = data.camera()
a16 = data.camera().astype(np.uint16)*4
ent8 = entropy(a8,disk(5)) # pixel value contain 10x the local entropy
ent16 = entropy(a16,disk(5)) # pixel value contain 1000x the local entropy
image = img_as_ubyte(data.camera())
# display results
plt.figure(figsize=(10, 10))
fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 4))
plt.subplot(2,2,1)
plt.imshow(a8, cmap=plt.cm.gray)
plt.xlabel('8-bit image')
plt.colorbar()
img0 = ax0.imshow(image, cmap=plt.cm.gray)
ax0.set_title('Image')
ax0.axis('off')
plt.colorbar(img0, ax=ax0)
plt.subplot(2,2,2)
plt.imshow(ent8, cmap=plt.cm.jet)
plt.xlabel('entropy*10')
plt.colorbar()
img1 = ax1.imshow(entropy(image, disk(5)), cmap=plt.cm.jet)
ax1.set_title('Entropy')
ax1.axis('off')
plt.colorbar(img1, ax=ax1)
plt.subplot(2,2,3)
plt.imshow(a16, cmap=plt.cm.gray)
plt.xlabel('16-bit image')
plt.colorbar()
plt.subplot(2,2,4)
plt.imshow(ent16, cmap=plt.cm.jet)
plt.xlabel('entropy*1000')
plt.colorbar()
plt.show()
doc/examples/plot_equalize.py
+2 -3
View File
@@ -17,13 +17,12 @@ that fall within the 2nd and 98th percentiles [2]_.
.. [2] http://homepages.inf.ed.ac.uk/rbf/HIPR2/stretch.htm
"""
import matplotlib.pyplot as plt
import numpy as np
from skimage import data, img_as_float
from skimage import exposure
import matplotlib.pyplot as plt
import numpy as np
def plot_img_and_hist(img, axes, bins=256):
"""Plot an image along with its histogram and cumulative histogram.
doc/examples/plot_gabor.py
+135
View File
@@ -0,0 +1,135 @@
"""
=============================================
Gabor filter banks for texture classification
=============================================
In this example, we will see how to classify textures based on Gabor filter
banks. Frequency and orientation representations of the Gabor filter are similar
to those of the human visual system.
The images are filtered using the real parts of various different Gabor filter
kernels. The mean and variance of the filtered images are then used as features
for classification, which is based on the least squared error for simplicity.
"""
from __future__ import print_function
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from scipy import ndimage as nd
from skimage import data
from skimage.util import img_as_float
from skimage.filter import gabor_kernel
matplotlib.rcParams['font.size'] = 9
def compute_feats(image, kernels):
feats = np.zeros((len(kernels), 2), dtype=np.double)
for k, kernel in enumerate(kernels):
filtered = nd.convolve(image, kernel, mode='wrap')
feats[k, 0] = filtered.mean()
feats[k, 1] = filtered.var()
return feats
def match(feats, ref_feats):
min_error = np.inf
min_i = None
for i in range(ref_feats.shape[0]):
error = np.sum((feats - ref_feats[i, :])**2)
if error < min_error:
min_error = error
min_i = i
return min_i
# prepare filter bank kernels
kernels = []
for theta in range(4):
theta = theta / 4. * np.pi
for sigma in (1, 3):
for frequency in (0.05, 0.25):
kernel = np.real(gabor_kernel(frequency, theta=theta,
sigma_x=sigma, sigma_y=sigma))
kernels.append(kernel)
shrink = (slice(0, None, 3), slice(0, None, 3))
brick = img_as_float(data.load('brick.png'))[shrink]
grass = img_as_float(data.load('grass.png'))[shrink]
wall = img_as_float(data.load('rough-wall.png'))[shrink]
image_names = ('brick', 'grass', 'wall')
images = (brick, grass, wall)
# prepare reference features
ref_feats = np.zeros((3, len(kernels), 2), dtype=np.double)
ref_feats[0, :, :] = compute_feats(brick, kernels)
ref_feats[1, :, :] = compute_feats(grass, kernels)
ref_feats[2, :, :] = compute_feats(wall, kernels)
print('Rotated images matched against references using Gabor filter banks:')
print('original: brick, rotated: 30deg, match result: ', end='')
feats = compute_feats(nd.rotate(brick, angle=190, reshape=False), kernels)
print(image_names[match(feats, ref_feats)])
print('original: brick, rotated: 70deg, match result: ', end='')
feats = compute_feats(nd.rotate(brick, angle=70, reshape=False), kernels)
print(image_names[match(feats, ref_feats)])
print('original: grass, rotated: 145deg, match result: ', end='')
feats = compute_feats(nd.rotate(grass, angle=145, reshape=False), kernels)
print(image_names[match(feats, ref_feats)])
def power(image, kernel):
# Normalize images for better comparison.
image = (image - image.mean()) / image.std()
return np.sqrt(nd.convolve(image, np.real(kernel), mode='wrap')**2 +
nd.convolve(image, np.imag(kernel), mode='wrap')**2)
# Plot a selection of the filter bank kernels and their responses.
results = []
kernel_params = []
for theta in (0, 1):
theta = theta / 4. * np.pi
for frequency in (0.1, 0.4):
kernel = gabor_kernel(frequency, theta=theta)
params = 'theta=%d,\nfrequency=%.2f' % (theta * 180 / np.pi, frequency)
kernel_params.append(params)
# Save kernel and the power image for each image
results.append((kernel, [power(img, kernel) for img in images]))
fig, axes = plt.subplots(nrows=5, ncols=4, figsize=(9, 6))
plt.gray()
fig.suptitle('Image responses for Gabor filter kernels', fontsize=15)
axes[0][0].axis('off')
# Plot original images
for label, img, ax in zip(image_names, images, axes[0][1:]):
ax.imshow(img)
ax.set_title(label)
ax.axis('off')
for label, (kernel, powers), ax_row in zip(kernel_params, results, axes[1:]):
# Plot Gabor kernel
ax = ax_row[0]
ax.imshow(np.real(kernel), interpolation='nearest')
ax.set_ylabel(label)
ax.set_xticks([])
ax.set_yticks([])
# Plot Gabor responses with the contrast normalized for each filter
vmin = np.min(powers)
vmax = np.max(powers)
for patch, ax in zip(powers, ax_row[1:]):
ax.imshow(patch, vmin=vmin, vmax=vmax)
ax.axis('off')
plt.show()
doc/examples/plot_gabors_from_lena.py
-3
View File
@@ -3,8 +3,6 @@
Gabors / Primary Visual Cortex "Simple Cells" from Lena
=======================================================
(under construction)
How to build a (bio-plausible) "sparse" dictionary (or 'codebook', or
'filterbank') for e.g. image classification without any fancy math and
with just standard python scientific libraries?
@@ -37,7 +35,6 @@ is not rocket science.
Interaction, and Functional Architecture in the Cat's Visual Cortex,
J. Physiol. 160 pp. 106-154 1962
"""
import numpy as np
from scipy.cluster.vq import kmeans2
from scipy import ndimage as ndi
doc/examples/plot_glcm.py
+2 -1
View File
@@ -19,10 +19,11 @@ this example) would be to train a classifier, such as logistic
regression, to label image patches from new images.
"""
import matplotlib.pyplot as plt
from skimage.feature import greycomatrix, greycoprops
from skimage import data
import matplotlib.pyplot as plt
PATCH_SIZE = 21
doc/examples/plot_hog.py
+4 -3
View File
@@ -1,4 +1,4 @@
r'''
"""
===============================
Histogram of Oriented Gradients
===============================
@@ -77,12 +77,13 @@ References
.. [2] David G. Lowe, "Distinctive image features from scale-invariant
keypoints," International Journal of Computer Vision, 60, 2 (2004),
pp. 91-110.
'''
"""
import matplotlib.pyplot as plt
from skimage.feature import hog
from skimage import data, color, exposure
import matplotlib.pyplot as plt
image = color.rgb2gray(data.lena())
doc/examples/plot_holes_and_peaks.py
+1
View File
@@ -10,6 +10,7 @@ image represent the maximum and minimum possible values of the reconstructed
image.
We start with an image containing both peaks and holes:
"""
import matplotlib.pyplot as plt
doc/examples/plot_hough_transform.py
-132
View File
@@ -1,132 +0,0 @@
r'''
===============
Hough transform
===============
The Hough transform in its simplest form is a `method to detect
straight lines <http://en.wikipedia.org/wiki/Hough_transform>`__.
In the following example, we construct an image with a line
intersection. We then use the Hough transform to explore a parameter
space for straight lines that may run through the image.
Algorithm overview
------------------
Usually, lines are parameterised as :math:`y = mx + c`, with a
gradient :math:`m` and y-intercept `c`. However, this would mean that
:math:`m` goes to infinity for vertical lines. Instead, we therefore
construct a segment perpendicular to the line, leading to the origin.
The line is represented by the length of that segment, :math:`r`, and
the angle it makes with the x-axis, :math:`\theta`.
The Hough transform constructs a histogram array representing the
parameter space (i.e., an :math:`M \times N` matrix, for :math:`M`
different values of the radius and :math:`N` different values of
:math:`\theta`). For each parameter combination, :math:`r` and
:math:`\theta`, we then find the number of non-zero pixels in the
input image that would fall close to the corresponding line, and
increment the array at position :math:`(r, \theta)` appropriately.
We can think of each non-zero pixel "voting" for potential line
candidates. The local maxima in the resulting histogram indicates the
parameters of the most probably lines. In our example, the maxima
occur at 45 and 135 degrees, corresponding to the normal vector
angles of each line.
Another approach is the Progressive Probabilistic Hough Transform
[1]_. It is based on the assumption that using a random subset of
voting points give a good approximation to the actual result, and that
lines can be extracted during the voting process by walking along
connected components. This returns the beginning and end of each
line segment, which is useful.
The function `probabilistic_hough` has three parameters: a general
threshold that is applied to the Hough accumulator, a minimum line
length and the line gap that influences line merging. In the example
below, we find lines longer than 10 with a gap less than 3 pixels.
References
----------
.. [1] C. Galamhos, J. Matas and J. Kittler,"Progressive probabilistic
Hough transform for line detection", in IEEE Computer Society
Conference on Computer Vision and Pattern Recognition, 1999.
.. [2] Duda, R. O. and P. E. Hart, "Use of the Hough Transformation to
Detect Lines and Curves in Pictures," Comm. ACM, Vol. 15,
pp. 11-15 (January, 1972)
'''
from skimage.transform import hough, hough_peaks, probabilistic_hough
from skimage.filter import canny
from skimage import data
import numpy as np
import matplotlib.pyplot as plt
# Construct test image
image = np.zeros((100, 100))
# Classic straight-line Hough transform
idx = np.arange(25, 75)
image[idx[::-1], idx] = 255
image[idx, idx] = 255
h, theta, d = hough(image)
plt.figure(figsize=(8, 4))
plt.subplot(131)
plt.imshow(image, cmap=plt.cm.gray)
plt.title('Input image')
plt.subplot(132)
plt.imshow(np.log(1 + h),
extent=[np.rad2deg(theta[-1]), np.rad2deg(theta[0]),
d[-1], d[0]],
cmap=plt.cm.gray, aspect=1/1.5)
plt.title('Hough transform')
plt.xlabel('Angles (degrees)')
plt.ylabel('Distance (pixels)')
plt.subplot(133)
plt.imshow(image, cmap=plt.cm.gray)
rows, cols = image.shape
for _, angle, dist in zip(*hough_peaks(h, theta, d)):
y0 = (dist - 0 * np.cos(angle)) / np.sin(angle)
y1 = (dist - cols * np.cos(angle)) / np.sin(angle)
plt.plot((0, cols), (y0, y1), '-r')
plt.axis((0, cols, rows, 0))
plt.title('Detected lines')
# Line finding, using the Probabilistic Hough Transform
image = data.camera()
edges = canny(image, 2, 1, 25)
lines = probabilistic_hough(edges, threshold=10, line_length=5, line_gap=3)
plt.figure(figsize=(8, 3))
plt.subplot(131)
plt.imshow(image, cmap=plt.cm.gray)
plt.title('Input image')
plt.subplot(132)
plt.imshow(edges, cmap=plt.cm.gray)
plt.title('Canny edges')
plt.subplot(133)
plt.imshow(edges * 0)
for line in lines:
p0, p1 = line
plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
plt.title('Probabilistic Hough')
plt.axis('image')
plt.show()
doc/examples/plot_ihc_color_separation.py
+72
View File
@@ -0,0 +1,72 @@
"""
==============================================
Immunohistochemical staining colors separation
==============================================
In this example we separate the immunohistochemical (IHC) staining from the
hematoxylin counterstaining. The separation is achieved with the method
described in [1]_, known as "color deconvolution".
The IHC staining expression of the FHL2 protein is here revealed with
Diaminobenzidine (DAB) which gives a brown color.
.. [1] A. C. Ruifrok and D. A. Johnston, "Quantification of histochemical
staining by color deconvolution.," Analytical and quantitative
cytology and histology / the International Academy of Cytology [and]
American Society of Cytology, vol. 23, no. 4, pp. 291-9, Aug. 2001.
"""
import matplotlib.pyplot as plt
from skimage import data
from skimage.color import rgb2hed
ihc_rgb = data.immunohistochemistry()
ihc_hed = rgb2hed(ihc_rgb)
fig, axes = plt.subplots(2, 2, figsize=(7, 6))
ax0, ax1, ax2, ax3 = axes.ravel()
ax0.imshow(ihc_rgb)
ax0.set_title("Original image")
ax1.imshow(ihc_hed[:, :, 0], cmap=plt.cm.gray)
ax1.set_title("Hematoxylin")
ax2.imshow(ihc_hed[:, :, 1], cmap=plt.cm.gray)
ax2.set_title("Eosin")
ax3.imshow(ihc_hed[:, :, 2], cmap=plt.cm.gray)
ax3.set_title("DAB")
for ax in axes.ravel():
ax.axis('off')
fig.subplots_adjust(hspace=0.3)
"""
.. image:: PLOT2RST.current_figure
Now we can easily manipulate the hematoxylin and DAB "channels":
"""
import numpy as np
from skimage.exposure import rescale_intensity
# Rescale hematoxylin and DAB signals and give them a fluorescence look
h = rescale_intensity(ihc_hed[:, :, 0], out_range=(0, 1))
d = rescale_intensity(ihc_hed[:, :, 2], out_range=(0, 1))
zdh = np.dstack((np.zeros_like(h), d, h))
plt.figure()
plt.imshow(zdh)
plt.title("Stain separated image (rescaled)")
plt.axis('off')
plt.show()
"""
.. image:: PLOT2RST.current_figure
"""
doc/examples/plot_join_segmentations.py
+12 -15
View File
@@ -8,19 +8,19 @@ segmentations. The `skimage.segmentation.join_segmentations` function
computes the join of two segmentations, in which a pixel is placed in
the same segment if and only if it is in the same segment in _both_
segmentations.
"""
"""
import numpy as np
from scipy import ndimage as nd
import matplotlib.pyplot as plt
import matplotlib as mpl
from skimage.filter import sobel
from skimage.segmentation import slic, join_segmentations
from skimage.morphology import watershed
from skimage.color import label2rgb
from skimage import data
coins = data.coins()
# make segmentation using edge-detection and watershed
@@ -43,24 +43,21 @@ seg2 = slic(coins_colour, n_segments=30, max_iter=160, sigma=1, ratio=9,
# combine the two
segj = join_segmentations(seg1, seg2)
### Display the result ###
# make a random colormap for a set number of values
def random_cmap(im):
np.random.seed(9)
cmap_array = np.concatenate(
(np.zeros((1, 3)), np.random.rand(np.ceil(im.max()), 3)))
return mpl.colors.ListedColormap(cmap_array)
# show the segmentations
fig, axes = plt.subplots(ncols=4, figsize=(9, 2.5))
axes[0].imshow(coins, cmap=plt.cm.gray, interpolation='nearest')
axes[0].set_title('Image')
axes[1].imshow(seg1, cmap=random_cmap(seg1), interpolation='nearest')
color1 = label2rgb(seg1, image=coins, bg_label=0)
axes[1].imshow(color1, interpolation='nearest')
axes[1].set_title('Sobel+Watershed')
axes[2].imshow(seg2, cmap=random_cmap(seg2), interpolation='nearest')
color2 = label2rgb(seg2, image=coins, image_alpha=0.5)
axes[2].imshow(color2, interpolation='nearest')
axes[2].set_title('SLIC superpixels')
axes[3].imshow(segj, cmap=random_cmap(segj), interpolation='nearest')
color3 = label2rgb(segj, image=coins, image_alpha=0.5)
axes[3].imshow(color3, interpolation='nearest')
axes[3].set_title('Join')
for ax in axes:
doc/examples/plot_label.py
+3 -2
View File
@@ -12,7 +12,6 @@ steps are applied:
4. Measure image regions to filter small objects
"""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
@@ -22,6 +21,7 @@ from skimage.filter import threshold_otsu
from skimage.segmentation import clear_border
from skimage.morphology import label, closing, square
from skimage.measure import regionprops
from skimage.color import label2rgb
image = data.coins()[50:-50, 50:-50]
@@ -38,9 +38,10 @@ clear_border(cleared)
label_image = label(cleared)
borders = np.logical_xor(bw, cleared)
label_image[borders] = -1
image_label_overlay = label2rgb(label_image, image=image)
fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(6, 6))
ax.imshow(label_image, cmap='jet')
ax.imshow(image_label_overlay)
for region in regionprops(label_image, ['Area', 'BoundingBox']):
doc/examples/plot_line_hough_transform.py
+130
View File
@@ -0,0 +1,130 @@
r"""
=============================
Straight line Hough transform
=============================
The Hough transform in its simplest form is a `method to detect straight lines
<http://en.wikipedia.org/wiki/Hough_transform>`__.
In the following example, we construct an image with a line intersection. We
then use the Hough transform to explore a parameter space for straight lines
that may run through the image.
Algorithm overview
------------------
Usually, lines are parameterised as :math:`y = mx + c`, with a gradient
:math:`m` and y-intercept `c`. However, this would mean that :math:`m` goes to
infinity for vertical lines. Instead, we therefore construct a segment
perpendicular to the line, leading to the origin. The line is represented by the
length of that segment, :math:`r`, and the angle it makes with the x-axis,
:math:`\theta`.
The Hough transform constructs a histogram array representing the parameter
space (i.e., an :math:`M \times N` matrix, for :math:`M` different values of the
radius and :math:`N` different values of :math:`\theta`). For each parameter
combination, :math:`r` and :math:`\theta`, we then find the number of non-zero
pixels in the input image that would fall close to the corresponding line, and
increment the array at position :math:`(r, \theta)` appropriately.
We can think of each non-zero pixel "voting" for potential line candidates. The
local maxima in the resulting histogram indicates the parameters of the most
probably lines. In our example, the maxima occur at 45 and 135 degrees,
corresponding to the normal vector angles of each line.
Another approach is the Progressive Probabilistic Hough Transform [1]_. It is
based on the assumption that using a random subset of voting points give a good
approximation to the actual result, and that lines can be extracted during the
voting process by walking along connected components. This returns the beginning
and end of each line segment, which is useful.
The function `probabilistic_hough` has three parameters: a general threshold
that is applied to the Hough accumulator, a minimum line length and the line gap
that influences line merging. In the example below, we find lines longer than 10
with a gap less than 3 pixels.
References
----------
.. [1] C. Galamhos, J. Matas and J. Kittler,"Progressive probabilistic
Hough transform for line detection", in IEEE Computer Society
Conference on Computer Vision and Pattern Recognition, 1999.
.. [2] Duda, R. O. and P. E. Hart, "Use of the Hough Transformation to
Detect Lines and Curves in Pictures," Comm. ACM, Vol. 15,
pp. 11-15 (January, 1972)
"""
from skimage.transform import (hough_line, hough_line_peaks,
probabilistic_hough_line)
from skimage.filter import canny
from skimage import data
import numpy as np
import matplotlib.pyplot as plt
# Construct test image
image = np.zeros((100, 100))
# Classic straight-line Hough transform
idx = np.arange(25, 75)
image[idx[::-1], idx] = 255
image[idx, idx] = 255
h, theta, d = hough_line(image)
plt.figure(figsize=(8, 4))
plt.subplot(131)
plt.imshow(image, cmap=plt.cm.gray)
plt.title('Input image')
plt.subplot(132)
plt.imshow(np.log(1 + h),
extent=[np.rad2deg(theta[-1]), np.rad2deg(theta[0]),
d[-1], d[0]],
cmap=plt.cm.gray, aspect=1/1.5)
plt.title('Hough transform')
plt.xlabel('Angles (degrees)')
plt.ylabel('Distance (pixels)')
plt.subplot(133)
plt.imshow(image, cmap=plt.cm.gray)
rows, cols = image.shape
for _, angle, dist in zip(*hough_line_peaks(h, theta, d)):
y0 = (dist - 0 * np.cos(angle)) / np.sin(angle)
y1 = (dist - cols * np.cos(angle)) / np.sin(angle)
plt.plot((0, cols), (y0, y1), '-r')
plt.axis((0, cols, rows, 0))
plt.title('Detected lines')
# Line finding, using the Probabilistic Hough Transform
image = data.camera()
edges = canny(image, 2, 1, 25)
lines = probabilistic_hough_line(edges, threshold=10, line_length=5, line_gap=3)
plt.figure(figsize=(8, 3))
plt.subplot(131)
plt.imshow(image, cmap=plt.cm.gray)
plt.title('Input image')
plt.subplot(132)
plt.imshow(edges, cmap=plt.cm.gray)
plt.title('Canny edges')
plt.subplot(133)
plt.imshow(edges * 0)
for line in lines:
p0, p1 = line
plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
plt.title('Probabilistic Hough')
plt.axis('image')
plt.show()
doc/examples/plot_local_binary_pattern.py
+170 -30
View File
@@ -4,24 +4,159 @@ Local Binary Pattern for texture classification
===============================================
In this example, we will see how to classify textures based on LBP (Local
Binary Pattern). The histogram of the LBP result is a good measure to classify
textures. For simplicity the histogram distributions are then tested against
each other using the Kullback-Leibler-Divergence.
Binary Pattern). LBP looks at points surrounding a central point and tests
whether the surrounding points are greater than or less than the central point
(i.e. gives a binary result).
Before trying out LBP on an image, it helps to look at a schematic of LBPs.
The below code is just used to plot the schematic.
"""
from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
METHOD = 'uniform'
plt.rcParams['font.size'] = 9
def plot_circle(ax, center, radius, color):
circle = plt.Circle(center, radius, facecolor=color, edgecolor='0.5')
ax.add_patch(circle)
def plot_lbp_model(ax, binary_values):
"""Draw the schematic for a local binary pattern."""
# Geometry spec
theta = np.deg2rad(45)
R = 1
r = 0.15
w = 1.5
gray = '0.5'
# Draw the central pixel.
plot_circle(ax, (0, 0), radius=r, color=gray)
# Draw the surrounding pixels.
for i, facecolor in enumerate(binary_values):
x = R * np.cos(i * theta)
y = R * np.sin(i * theta)
plot_circle(ax, (x, y), radius=r, color=str(facecolor))
# Draw the pixel grid.
for x in np.linspace(-w, w, 4):
ax.axvline(x, color=gray)
ax.axhline(x, color=gray)
# Tweak the layout.
ax.axis('image')
ax.axis('off')
size = w + 0.2
ax.set_xlim(-size, size)
ax.set_ylim(-size, size)
fig, axes = plt.subplots(ncols=5, figsize=(7, 2))
titles = ['flat', 'flat', 'edge', 'corner', 'non-uniform']
binary_patterns = [np.zeros(8),
np.ones(8),
np.hstack([np.ones(4), np.zeros(4)]),
np.hstack([np.zeros(3), np.ones(5)]),
[1, 0, 0, 1, 1, 1, 0, 0]]
for ax, values, name in zip(axes, binary_patterns, titles):
plot_lbp_model(ax, values)
ax.set_title(name)
"""
.. image:: PLOT2RST.current_figure
The figure above shows example results with black (or white) representing
pixels that are less (or more) intense than the central pixel. When surrounding
pixels are all black or all white, then that image region is flat (i.e.
featureless). Groups of continuous black or white pixels are considered
"uniform" patterns that can be interpreted as corners or edges. If pixels
switch back-and-forth between black and white pixels, the pattern is considered
"non-uniform".
When using LBP to detect texture, you measure a collection of LBPs over an
image patch and look at the distribution of these LBPs. Lets apply LBP to
a brick texture.
"""
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import scipy.ndimage as nd
import skimage.feature as ft
from skimage.transform import rotate
from skimage.feature import local_binary_pattern
from skimage import data
from skimage.color import label2rgb
# settings for LBP
METHOD = 'uniform'
P = 16
R = 2
matplotlib.rcParams['font.size'] = 9
radius = 3
n_points = 8 * radius
def overlay_labels(image, lbp, labels):
mask = np.logical_or.reduce([lbp == each for each in labels])
return label2rgb(mask, image=image, bg_label=0, alpha=0.5)
def highlight_bars(bars, indexes):
for i in indexes:
bars[i].set_facecolor('r')
image = data.load('brick.png')
lbp = local_binary_pattern(image, n_points, radius, METHOD)
def hist(ax, lbp):
n_bins = lbp.max() + 1
return ax.hist(lbp.ravel(), normed=True, bins=n_bins, range=(0, n_bins),
facecolor='0.5')
# plot histograms of LBP of textures
fig, (ax_img, ax_hist) = plt.subplots(nrows=2, ncols=3, figsize=(9, 6))
plt.gray()
titles = ('edge', 'flat', 'corner')
w = width = radius - 1
edge_labels = range(n_points // 2 - w, n_points // 2 + w + 1)
flat_labels = list(range(0, w + 1)) + list(range(n_points - w, n_points + 2))
i_14 = n_points // 4 # 1/4th of the histogram
i_34 = 3 * (n_points // 4) # 3/4th of the histogram
corner_labels = (list(range(i_14 - w, i_14 + w + 1)) +
list(range(i_34 - w, i_34 + w + 1)))
label_sets = (edge_labels, flat_labels, corner_labels)
for ax, labels in zip(ax_img, label_sets):
ax.imshow(overlay_labels(image, lbp, labels))
for ax, labels, name in zip(ax_hist, label_sets, titles):
counts, _, bars = hist(ax, lbp)
highlight_bars(bars, labels)
ax.set_ylim(ymax=np.max(counts[:-1]))
ax.set_xlim(xmax=n_points + 2)
ax.set_title(name)
ax_hist[0].set_ylabel('Percentage')
for ax in ax_img:
ax.axis('off')
"""
.. image:: PLOT2RST.current_figure
The above plot highlights flat, edge-like, and corner-like regions of the
image.
The histogram of the LBP result is a good measure to classify textures. Here,
we test the histogram distributions against each other using the
Kullback-Leibler-Divergence.
"""
# settings for LBP
radius = 2
n_points = 8 * radius
def kullback_leibler_divergence(p, q):
@@ -34,11 +169,12 @@ def kullback_leibler_divergence(p, q):
def match(refs, img):
best_score = 10
best_name = None
lbp = ft.local_binary_pattern(img, P, R, METHOD)
hist, _ = np.histogram(lbp, normed=True, bins=P + 2, range=(0, P + 2))
lbp = local_binary_pattern(img, n_points, radius, METHOD)
n_bins = lbp.max() + 1
hist, _ = np.histogram(lbp, normed=True, bins=n_bins, range=(0, n_bins))
for name, ref in refs.items():
ref_hist, _ = np.histogram(ref, normed=True, bins=P + 2,
range=(0, P + 2))
ref_hist, _ = np.histogram(ref, normed=True, bins=n_bins,
range=(0, n_bins))
score = kullback_leibler_divergence(hist, ref_hist)
if score < best_score:
best_score = score
@@ -51,19 +187,19 @@ grass = data.load('grass.png')
wall = data.load('rough-wall.png')
refs = {
'brick': ft.local_binary_pattern(brick, P, R, METHOD),
'grass': ft.local_binary_pattern(grass, P, R, METHOD),
'wall': ft.local_binary_pattern(wall, P, R, METHOD)
'brick': local_binary_pattern(brick, n_points, radius, METHOD),
'grass': local_binary_pattern(grass, n_points, radius, METHOD),
'wall': local_binary_pattern(wall, n_points, radius, METHOD)
}
# classify rotated textures
print 'Rotated images matched against references using LBP:'
print 'original: brick, rotated: 30deg, match result:',
print match(refs, nd.rotate(brick, angle=30, reshape=False))
print 'original: brick, rotated: 70deg, match result:',
print match(refs, nd.rotate(brick, angle=70, reshape=False))
print 'original: grass, rotated: 145deg, match result:',
print match(refs, nd.rotate(grass, angle=145, reshape=False))
print('Rotated images matched against references using LBP:')
print('original: brick, rotated: 30deg, match result: ',
match(refs, rotate(brick, angle=30, resize=False)))
print('original: brick, rotated: 70deg, match result: ',
match(refs, rotate(brick, angle=70, resize=False)))
print('original: grass, rotated: 145deg, match result: ',
match(refs, rotate(grass, angle=145, resize=False)))
# plot histograms of LBP of textures
fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows=2, ncols=3,
@@ -72,16 +208,20 @@ plt.gray()
ax1.imshow(brick)
ax1.axis('off')
ax4.hist(refs['brick'].ravel(), normed=True, bins=P + 2, range=(0, P + 2))
hist(ax4, refs['brick'])
ax4.set_ylabel('Percentage')
ax2.imshow(grass)
ax2.axis('off')
ax5.hist(refs['grass'].ravel(), normed=True, bins=P + 2, range=(0, P + 2))
hist(ax5, refs['grass'])
ax5.set_xlabel('Uniform LBP values')
ax3.imshow(wall)
ax3.axis('off')
ax6.hist(refs['wall'].ravel(), normed=True, bins=P + 2, range=(0, P + 2))
hist(ax6, refs['wall'])
"""
.. image:: PLOT2RST.current_figure
"""
plt.show()
doc/examples/plot_local_equalize.py
+17 -12
View File
@@ -1,27 +1,32 @@
"""
===============================
============================
Local Histogram Equalization
===============================
============================
This examples enhances an image with low contrast, using a method called
*local histogram equalization*, which "spreads out the most frequent intensity
values" in an image .
The equalized image [1]_ has a roughly linear cumulative distribution function for each pixel neighborhood.
The local version [2]_ of the histogram equalization emphasized every local graylevel variations.
This examples enhances an image with low contrast, using a method called *local
histogram equalization*, which spreads out the most frequent intensity values in
an image.
The equalized image [1]_ has a roughly linear cumulative distribution function
for each pixel neighborhood.
The local version [2]_ of the histogram equalization emphasized every local
graylevel variations.
References
----------
.. [1] http://en.wikipedia.org/wiki/Histogram_equalization
.. [2] http://en.wikipedia.org/wiki/Adaptive_histogram_equalization
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage import data
from skimage.util.dtype import dtype_range
from skimage.util import img_as_ubyte
from skimage import exposure
from skimage.morphology import disk
import matplotlib.pyplot as plt
import numpy as np
from skimage.filter import rank
@@ -52,7 +57,7 @@ def plot_img_and_hist(img, axes, bins=256):
# Load an example image
img = data.moon()
img = img_as_ubyte(data.moon())
# Contrast stretching
p2 = np.percentile(img, 2)
doc/examples/plot_local_otsu.py
+14 -12
View File
@@ -1,14 +1,16 @@
"""
=====================
====================
Local Otsu Threshold
=====================
This example shows how Otsu's threshold [1]_ method can be applied locally.
For each pixel, an "optimal" threshold is determined by maximizing the variance between two classes of pixels
of the local neighborhood defined by a structuring element.
====================
This example shows how Otsu's threshold [1]_ method can be applied locally. For
each pixel, an "optimal" threshold is determined by maximizing the variance
between two classes of pixels of the local neighborhood defined by a structuring
element.
The example compares the local threshold with the global threshold.
.. note: local threshold is much slower than global one.
.. note: local is much slower than global thresholding
.. [1] http://en.wikipedia.org/wiki/Otsu's_method
@@ -16,12 +18,12 @@ The example compares the local threshold with the global threshold.
import matplotlib.pyplot as plt
from skimage import data
from skimage.morphology.selem import disk
import skimage.filter.rank as rank
from skimage.filter import threshold_otsu
from skimage.morphology import disk
from skimage.filter import threshold_otsu, rank
from skimage.util import img_as_ubyte
p8 = data.page()
p8 = img_as_ubyte(data.page())
radius = 10
selem = disk(radius)
@@ -42,8 +44,8 @@ plt.xlabel('local Otsu ($radius=%d$)' % radius)
plt.colorbar()
plt.subplot(2, 2, 3)
plt.imshow(p8 >= loc_otsu, cmap=plt.cm.gray)
plt.xlabel('original>=local Otsu' % t_glob_otsu)
plt.xlabel('original >= local Otsu' % t_glob_otsu)
plt.subplot(2, 2, 4)
plt.imshow(glob_otsu, cmap=plt.cm.gray)
plt.xlabel('global Otsu ($t=%d$)' % t_glob_otsu)
plt.xlabel('global Otsu ($t = %d$)' % t_glob_otsu)
plt.show()
doc/examples/plot_marching_cubes.py
+56
View File
@@ -0,0 +1,56 @@
"""
==============
Marching Cubes
==============
Marching cubes is an algorithm to extract a 2D surface mesh from a 3D volume.
This can be conceptualized as a 3D generalization of isolines on topographical
or weather maps. It works by iterating across the volume, looking for regions
which cross the level of interest. If such regions are found, triangulations
are generated and added to an output mesh. The final result is a set of
vertices and a set of triangular faces.
The algorithm requires a data volume and an isosurface value. For example, in
CT imaging Hounsfield units of +700 to +3000 represent bone. So, one potential
input would be a reconstructed CT set of data and the value +700, to extract
a mesh for regions of bone or bone-like density.
This implementation also works correctly on anisotropic datasets, where the
voxel spacing is not equal for every spatial dimension, through use of the
`sampling` kwarg.
"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from skimage import measure
from skimage.draw import ellipsoid
# Generate a level set about zero of two identical ellipsoids in 3D
ellip_base = ellipsoid(6, 10, 16, levelset=True)
ellip_double = np.concatenate((ellip_base[:-1, ...],
ellip_base[2:, ...]), axis=0)
# Use marching cubes to obtain the surface mesh of these ellipsoids
verts, faces = measure.marching_cubes(ellip_double, 0)
# Display resulting triangular mesh using Matplotlib. This can also be done
# with mayavi (see skimage.measure.marching_cubes docstring).
fig = plt.figure(figsize=(10, 12))
ax = fig.add_subplot(111, projection='3d')
# Fancy indexing: `verts[faces]` to generate a collection of triangles
mesh = Poly3DCollection(verts[faces])
ax.add_collection3d(mesh)
ax.set_xlabel("x-axis: a = 6 per ellipsoid")
ax.set_ylabel("y-axis: b = 10")
ax.set_zlabel("z-axis: c = 16")
ax.set_xlim(0, 24) # a = 6 (times two for 2nd ellipsoid)
ax.set_ylim(0, 20) # b = 10
ax.set_zlim(0, 32) # c = 16
plt.show()
doc/examples/plot_marked_watershed.py
+6 -5
View File
@@ -1,7 +1,7 @@
"""
================================
===============================
Markers for watershed transform
================================
===============================
The watershed is a classical algorithm used for **segmentation**, that
is, for separating different objects in an image.
@@ -16,13 +16,14 @@ See Wikipedia_ for more details on the algorithm.
from scipy import ndimage
import matplotlib.pyplot as plt
from skimage.morphology import watershed, disk
from skimage import data
# original data
from skimage.filter import rank
from skimage.util import img_as_ubyte
image = data.camera()
image = img_as_ubyte(data.camera())
# denoise image
denoised = rank.median(image, disk(2))
doc/examples/plot_matching.py
+144
View File
@@ -0,0 +1,144 @@
"""
============================
Robust matching using RANSAC
============================
In this simplified example we first generate two synthetic images as if they
were taken from different view points.
In the next step we find interest points in both images and find
correspondences based on a weighted sum of squared differences of a small
neighborhood around them. Note, that this measure is only robust towards
linear radiometric and not geometric distortions and is thus only usable with
slight view point changes.
After finding the correspondences we end up having a set of source and
destination coordinates which can be used to estimate the geometric
transformation between both images. However, many of the correspondences are
faulty and simply estimating the parameter set with all coordinates is not
sufficient. Therefore, the RANSAC algorithm is used on top of the normal model
to robustly estimate the parameter set by detecting outliers.
"""
from __future__ import print_function
import numpy as np
from matplotlib import pyplot as plt
from skimage import data
from skimage.util import img_as_float
from skimage.feature import corner_harris, corner_subpix, corner_peaks
from skimage.transform import warp, AffineTransform
from skimage.exposure import rescale_intensity
from skimage.color import rgb2gray
from skimage.measure import ransac
# generate synthetic checkerboard image and add gradient for the later matching
checkerboard = img_as_float(data.checkerboard())
img_orig = np.zeros(list(checkerboard.shape) + [3])
img_orig[..., 0] = checkerboard
gradient_r, gradient_c = np.mgrid[0:img_orig.shape[0],
0:img_orig.shape[1]] / float(img_orig.shape[0])
img_orig[..., 1] = gradient_r
img_orig[..., 2] = gradient_c
img_orig = rescale_intensity(img_orig)
img_orig_gray = rgb2gray(img_orig)
# warp synthetic image
tform = AffineTransform(scale=(0.9, 0.9), rotation=0.2, translation=(20, -10))
img_warped = warp(img_orig, tform.inverse, output_shape=(200, 200))
img_warped_gray = rgb2gray(img_warped)
# extract corners using Harris' corner measure
coords_orig = corner_peaks(corner_harris(img_orig_gray), threshold_rel=0.001,
min_distance=5)
coords_warped = corner_peaks(corner_harris(img_warped_gray),
threshold_rel=0.001, min_distance=5)
# determine sub-pixel corner position
coords_orig_subpix = corner_subpix(img_orig_gray, coords_orig, window_size=9)
coords_warped_subpix = corner_subpix(img_warped_gray, coords_warped,
window_size=9)
def gaussian_weights(window_ext, sigma=1):
y, x = np.mgrid[-window_ext:window_ext+1, -window_ext:window_ext+1]
g = np.zeros(y.shape, dtype=np.double)
g[:] = np.exp(-0.5 * (x**2 / sigma**2 + y**2 / sigma**2))
g /= 2 * np.pi * sigma * sigma
return g
def match_corner(coord, window_ext=5):
r, c = np.round(coord)
window_orig = img_orig[r-window_ext:r+window_ext+1,
c-window_ext:c+window_ext+1, :]
# weight pixels depending on distance to center pixel
weights = gaussian_weights(window_ext, 3)
weights = np.dstack((weights, weights, weights))
# compute sum of squared differences to all corners in warped image
SSDs = []
for cr, cc in coords_warped:
window_warped = img_warped[cr-window_ext:cr+window_ext+1,
cc-window_ext:cc+window_ext+1, :]
SSD = np.sum(weights * (window_orig - window_warped)**2)
SSDs.append(SSD)
# use corner with minimum SSD as correspondence
min_idx = np.argmin(SSDs)
return coords_warped_subpix[min_idx]
# find correspondences using simple weighted sum of squared differences
src = []
dst = []
for coord in coords_orig_subpix:
src.append(coord)
dst.append(match_corner(coord))
src = np.array(src)
dst = np.array(dst)
# estimate affine transform model using all coordinates
model = AffineTransform()
model.estimate(src, dst)
# robustly estimate affine transform model with RANSAC
model_robust, inliers = ransac((src, dst), AffineTransform, min_samples=3,
residual_threshold=2, max_trials=100)
outliers = inliers == False
# compare "true" and estimated transform parameters
print(tform.scale, tform.translation, tform.rotation)
print(model.scale, model.translation, model.rotation)
print(model_robust.scale, model_robust.translation, model_robust.rotation)
# visualize correspondences
img_combined = np.concatenate((img_orig_gray, img_warped_gray), axis=1)
fig, ax = plt.subplots(nrows=2, ncols=1)
plt.gray()
ax[0].imshow(img_combined, interpolation='nearest')
ax[0].axis('off')
ax[0].axis((0, 400, 200, 0))
ax[0].set_title('Correct correspondences')
ax[1].imshow(img_combined, interpolation='nearest')
ax[1].axis('off')
ax[1].axis((0, 400, 200, 0))
ax[1].set_title('Faulty correspondences')
for ax_idx, (m, color) in enumerate(((inliers, 'g'), (outliers, 'r'))):
ax[ax_idx].plot((src[m, 1], dst[m, 1] + 200), (src[m, 0], dst[m, 0]), '-',
color=color)
ax[ax_idx].plot(src[m, 1], src[m, 0], '.', markersize=10, color=color)
ax[ax_idx].plot(dst[m, 1] + 200, dst[m, 0], '.', markersize=10,
color=color)
plt.show()
doc/examples/plot_medial_transform.py
+5 -7
View File
@@ -3,7 +3,7 @@
Medial axis skeletonization
===========================
The medial axis of an object is the set of all points having more than one
The medial axis of an object is the set of all points having more than one
closest point on the object's boundary. It is often called the **topological
skeleton**, because it is a 1-pixel wide skeleton of the object, with the same
connectivity as the original object.
@@ -15,11 +15,11 @@ argument ``return_distance=True``), it is possible to compute the distance to
the background for all points of the medial axis with this function. This gives
an estimate of the local width of the objects.
For a skeleton with fewer branches, there exists another skeletonization
For a skeleton with fewer branches, there exists another skeletonization
algorithm in ``skimage``: ``skimage.morphology.skeletonize``, that computes
a skeleton by iterative morphological thinnings.
"""
"""
import numpy as np
from scipy import ndimage
from skimage.morphology import medial_axis
@@ -33,7 +33,7 @@ def microstructure(l=256):
Parameters
----------
l: int, optional
l: int, optional
linear size of the returned image
"""
@@ -64,7 +64,5 @@ plt.imshow(dist_on_skel, cmap=plt.cm.spectral, interpolation='nearest')
plt.contour(data, [0.5], colors='w')
plt.axis('off')
plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
right=1)
plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
plt.show()
doc/examples/plot_otsu.py
-2
View File
@@ -14,7 +14,6 @@ the intra-class variance.
.. [1] http://en.wikipedia.org/wiki/Otsu's_method
"""
import matplotlib.pyplot as plt
from skimage.data import camera
@@ -42,4 +41,3 @@ plt.title('Thresholded')
plt.axis('off')
plt.show()
doc/examples/plot_peak_local_max.py
+2 -3
View File
@@ -1,7 +1,7 @@
"""
===============================================================================
====================
Finding local maxima
===============================================================================
====================
The ``peak_local_max`` function returns the coordinates of local peaks (maxima)
in an image. A maximum filter is used for finding local maxima. This operation
@@ -47,4 +47,3 @@ plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9,
bottom=0.02, left=0.02, right=0.98)
plt.show()
doc/examples/plot_piecewise_affine.py
+1 -1
View File
@@ -4,8 +4,8 @@ Piecewise Affine Transformation
===============================
This example shows how to use the Piecewise Affine Transformation.
"""
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage.transform import PiecewiseAffineTransform, warp
doc/examples/plot_polygon.py
+5 -2
View File
@@ -5,10 +5,13 @@ Approximate and subdivide polygons
This example shows how to approximate (Douglas-Peucker algorithm) and subdivide
(B-Splines) polygonal chains.
"""
from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
from skimage.draw import ellipse
from skimage.measure import find_contours, approximate_polygon, \
subdivide_polygon
@@ -45,7 +48,7 @@ for _ in range(5):
# approximate subdivided polygon with Douglas-Peucker algorithm
appr_hand = approximate_polygon(new_hand, tolerance=0.02)
print "Number of coordinates:", len(hand), len(new_hand), len(appr_hand)
print("Number of coordinates:", len(hand), len(new_hand), len(appr_hand))
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(9, 4))
@@ -70,7 +73,7 @@ for contour in find_contours(img, 0):
ax2.plot(coords[:, 1], coords[:, 0], '-r', linewidth=2)
coords2 = approximate_polygon(contour, tolerance=39.5)
ax2.plot(coords2[:, 1], coords2[:, 0], '-g', linewidth=2)
print "Number of coordinates:", len(contour), len(coords), len(coords2)
print("Number of coordinates:", len(contour), len(coords), len(coords2))
ax2.axis((0, 800, 0, 800))
doc/examples/plot_pyramid.py
-1
View File
@@ -9,7 +9,6 @@ implement algorithms for denoising, texture discrimination, and scale- invariant
detection.
"""
import numpy as np
import matplotlib.pyplot as plt
doc/examples/plot_radon_transform.py
+173 -31
View File
@@ -3,55 +3,197 @@
Radon transform
===============
The radon transform is a technique widely used in tomography to
reconstruct an object from different projections. A projection is, for
example, the scattering data obtained as the output of a tomographic
scan.
In computed tomography, the tomography reconstruction problem is to obtain
a tomographic slice image from a set of projections [1]_. A projection is formed
by drawing a set of parallel rays through the 2D object of interest, assigning
the integral of the object's contrast along each ray to a single pixel in the
projection. A single projection of a 2D object is one dimensional. To
enable computed tomography reconstruction of the object, several projections
must be acquired, each of them corresponding to a different angle between the
rays with respect to the object. A collection of projections at several angles
is called a sinogram, which is a linear transform of the original image.
For more information see:
The inverse Radon transform is used in computed tomography to reconstruct
a 2D image from the measured projections (the sinogram). A practical, exact
implementation of the inverse Radon transform does not exist, but there are
several good approximate algorithms available.
- http://en.wikipedia.org/wiki/Radon_transform
- http://www.clear.rice.edu/elec431/projects96/DSP/bpanalysis.html
As the inverse Radon transform reconstructs the object from a set of
projections, the (forward) Radon transform can be used to simulate a
tomography experiment.
This script performs the radon transform, and reconstructs the
input image based on the resulting sinogram.
This script performs the Radon transform to simulate a tomography experiment
and reconstructs the input image based on the resulting sinogram formed by
the simulation. Two methods for performing the inverse Radon transform
and reconstructing the original image are compared: The Filtered Back
Projection (FBP) and the Simultaneous Algebraic Reconstruction
Technique (SART).
.. seealso::
- AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging",
http://www.slaney.org/pct/pct-toc.html
- http://en.wikipedia.org/wiki/Radon_transform
The forward transform
=====================
As our original image, we will use the Shepp-Logan phantom. When calculating
the Radon transform, we need to decide how many projection angles we wish
to use. As a rule of thumb, the number of projections should be about the
same as the number of pixels there are across the object (to see why this
is so, consider how many unknown pixel values must be determined in the
reconstruction process and compare this to the number of measurements
provided by the projections), and we follow that rule here. Below is the
original image and its Radon transform, often known as its _sinogram_:
"""
from __future__ import print_function, division
import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread
from skimage import data_dir
from skimage.transform import radon, iradon
from scipy.ndimage import zoom
from skimage.transform import radon, rescale
image = imread(data_dir + "/phantom.png", as_grey=True)
image = zoom(image, 0.4)
image = rescale(image, scale=0.4)
plt.figure(figsize=(8, 4.5))
plt.subplot(121)
plt.title("Original")
plt.imshow(image, cmap=plt.cm.Greys_r)
theta = np.linspace(0., 180., max(image.shape), endpoint=True)
sinogram = radon(image, theta=theta, circle=True)
plt.subplot(122)
plt.title("Radon transform\n(Sinogram)")
plt.xlabel("Projection angle (deg)")
plt.ylabel("Projection position (pixels)")
plt.imshow(sinogram, cmap=plt.cm.Greys_r,
extent=(0, 180, 0, sinogram.shape[0]), aspect='auto')
plt.subplots_adjust(hspace=0.4, wspace=0.5)
plt.show()
"""
.. image:: PLOT2RST.current_figure
Reconstruction with the Filtered Back Projection (FBP)
======================================================
The mathematical foundation of the filtered back projection is the Fourier
slice theorem [2]_. It uses Fourier transform of the projection and
interpolation in Fourier space to obtain the 2D Fourier transform of the image,
which is then inverted to form the reconstructed image. The filtered back
projection is among the fastest methods of performing the inverse Radon
transform. The only tunable parameter for the FBP is the filter, which is
applied to the Fourier transformed projections. It may be used to suppress
high frequency noise in the reconstruction. ``skimage`` provides a few
different options for the filter.
"""
from skimage.transform import iradon
reconstruction_fbp = iradon(sinogram, theta=theta, circle=True)
error = reconstruction_fbp - image
print('FBP rms reconstruction error: %.3g' % np.sqrt(np.mean(error**2)))
imkwargs = dict(vmin=-0.2, vmax=0.2)
plt.figure(figsize=(8, 4.5))
plt.subplot(121)
plt.title("Reconstruction\nFiltered back projection")
plt.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r)
plt.subplot(122)
plt.title("Reconstruction error\nFiltered back projection")
plt.imshow(reconstruction_fbp - image, cmap=plt.cm.Greys_r, **imkwargs)
plt.show()
"""
.. image:: PLOT2RST.current_figure
Reconstruction with the Simultaneous Algebraic Reconstruction Technique
=======================================================================
Algebraic reconstruction techniques for tomography are based on a
straightforward idea: for a pixelated image the value of a single ray in a
particular projection is simply a sum of all the pixels the ray passes through
on its way through the object. This is a way of expressing the forward Radon
transform. The inverse Radon transform can then be formulated as a (large) set
of linear equations. As each ray passes through a small fraction of the pixels
in the image, this set of equations is sparse, allowing iterative solvers for
sparse linear systems to tackle the system of equations. One iterative method
has been particularly popular, namely Kaczmarz' method [3]_, which has the
property that the solution will approach a least-squares solution of the
equation set.
The combination of the formulation of the reconstruction problem as a set
of linear equations and an iterative solver makes algebraic techniques
relatively flexible, hence some forms of prior knowledge can be incorporated
with relative ease.
``skimage`` provides one of the more popular variations of the algebraic
reconstruction techniques: the Simultaneous Algebraic Reconstruction Technique
(SART) [1]_ [4]_. It uses Kaczmarz' method [3]_ as the iterative solver. A good
reconstruction is normally obtained in a single iteration, making the method
computationally effective. Running one or more extra iterations will normally
improve the reconstruction of sharp, high frequency features and reduce the
mean squared error at the expense of increased high frequency noise (the user
will need to decide on what number of iterations is best suited to the problem
at hand. The implementation in ``skimage`` allows prior information of the
form of a lower and upper threshold on the reconstructed values to be supplied
to the reconstruction.
"""
from skimage.transform import iradon_sart
reconstruction_sart = iradon_sart(sinogram, theta=theta)
error = reconstruction_sart - image
print('SART (1 iteration) rms reconstruction error: %.3g'
% np.sqrt(np.mean(error**2)))
plt.figure(figsize=(8, 8.5))
plt.subplot(221)
plt.title("Original");
plt.imshow(image, cmap=plt.cm.Greys_r)
plt.title("Reconstruction\nSART")
plt.imshow(reconstruction_sart, cmap=plt.cm.Greys_r)
plt.subplot(222)
projections = radon(image, theta=[0, 45, 90])
plt.plot(projections);
plt.title("Projections at\n0, 45 and 90 degrees")
plt.xlabel("Projection axis");
plt.ylabel("Intensity");
plt.title("Reconstruction error\nSART")
plt.imshow(reconstruction_sart - image, cmap=plt.cm.Greys_r, **imkwargs)
# Run a second iteration of SART by supplying the reconstruction
# from the first iteration as an initial estimate
reconstruction_sart2 = iradon_sart(sinogram, theta=theta,
image=reconstruction_sart)
error = reconstruction_sart2 - image
print('SART (2 iterations) rms reconstruction error: %.3g'
% np.sqrt(np.mean(error**2)))
projections = radon(image)
plt.subplot(223)
plt.title("Radon transform\n(Sinogram)");
plt.xlabel("Projection axis");
plt.ylabel("Intensity");
plt.imshow(projections)
reconstruction = iradon(projections)
plt.title("Reconstruction\nSART, 2 iterations")
plt.imshow(reconstruction_sart2, cmap=plt.cm.Greys_r)
plt.subplot(224)
plt.title("Reconstruction\nfrom sinogram")
plt.imshow(reconstruction, cmap=plt.cm.Greys_r)
plt.subplots_adjust(hspace=0.4, wspace=0.5)
plt.title("Reconstruction error\nSART, 2 iterations")
plt.imshow(reconstruction_sart2 - image, cmap=plt.cm.Greys_r, **imkwargs)
plt.show()
"""
.. image:: PLOT2RST.current_figure
.. [1] AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging",
IEEE Press 1988. http://www.slaney.org/pct/pct-toc.html
.. [2] Wikipedia, Radon transform,
http://en.wikipedia.org/wiki/Radon_transform#Relationship_with_the_Fourier_transform
.. [3] S Kaczmarz, "Angenaeherte Aufloesung von Systemen linearer
Gleichungen", Bulletin International de l'Academie Polonaise des
Sciences et des Lettres 35 pp 355--357 (1937)
.. [4] AH Andersen, AC Kak, "Simultaneous algebraic reconstruction technique
(SART): a superior implementation of the ART algorithm", Ultrasonic
Imaging 6 pp 81--94 (1984)
"""
doc/examples/plot_random_walker_segmentation.py
+3 -2
View File
@@ -18,13 +18,14 @@ values, and use the random walker for the segmentation.
.. [1] *Random walks for image segmentation*, Leo Grady, IEEE Trans. Pattern
Anal. Mach. Intell. 2006 Nov; 28(11):1768-83
"""
"""
import numpy as np
from scipy import ndimage
from skimage.segmentation import random_walker
import matplotlib.pyplot as plt
from skimage.segmentation import random_walker
def microstructure(l=256):
"""
doc/examples/plot_rank_mean.py
+52
View File
@@ -0,0 +1,52 @@
"""
============
Mean filters
============
This example compares the following mean filters of the rank filter package:
* **local mean**: all pixels belonging to the structuring element to compute
average gray level.
* **percentile mean**: only use values between percentiles p0 and p1
(here 10% and 90%).
* **bilateral mean**: only use pixels of the structuring element having a gray
level situated inside g-s0 and g+s1 (here g-500 and g+500)
Percentile and usual mean give here similar results, these filters smooth the
complete image (background and details). Bilateral mean exhibits a high
filtering rate for continuous area (i.e. background) while higher image
frequencies remain untouched.
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage import data
from skimage.morphology import disk
from skimage.filter import rank
image = (data.coins()).astype(np.uint16) * 16
selem = disk(20)
percentile_result = rank.mean_percentile(image, selem=selem, p0=.1, p1=.9)
bilateral_result = rank.mean_bilateral(image, selem=selem, s0=500, s1=500)
normal_result = rank.mean(image, selem=selem)
fig, axes = plt.subplots(nrows=3, figsize=(8, 10))
ax0, ax1, ax2 = axes
ax0.imshow(np.hstack((image, percentile_result)))
ax0.set_title('Percentile mean')
ax0.axis('off')
ax1.imshow(np.hstack((image, bilateral_result)))
ax1.set_title('Bilateral mean')
ax1.axis('off')
ax2.imshow(np.hstack((image, normal_result)))
ax2.set_title('Local mean')
ax2.axis('off')
plt.show()
doc/examples/plot_ransac.py
+55
View File
@@ -0,0 +1,55 @@
"""
=========================================
Robust line model estimation using RANSAC
=========================================
In this example we see how to robustly fit a line model to faulty data using
the RANSAC algorithm.
"""
import numpy as np
from matplotlib import pyplot as plt
from skimage.measure import LineModel, ransac
np.random.seed(seed=1)
# generate coordinates of line
x = np.arange(-200, 200)
y = 0.2 * x + 20
data = np.column_stack([x, y])
# add faulty data
faulty = np.array(30 * [(180., -100)])
faulty += 5 * np.random.normal(size=faulty.shape)
data[:faulty.shape[0]] = faulty
# add gaussian noise to coordinates
noise = np.random.normal(size=data.shape)
data += 0.5 * noise
data[::2] += 5 * noise[::2]
data[::4] += 20 * noise[::4]
# fit line using all data
model = LineModel()
model.estimate(data)
# robustly fit line only using inlier data with RANSAC algorithm
model_robust, inliers = ransac(data, LineModel, min_samples=2,
residual_threshold=1, max_trials=1000)
outliers = inliers == False
# generate coordinates of estimated models
line_x = np.arange(-250, 250)
line_y = model.predict_y(line_x)
line_y_robust = model_robust.predict_y(line_x)
plt.plot(data[inliers, 0], data[inliers, 1], '.b', alpha=0.6,
label='Inlier data')
plt.plot(data[outliers, 0], data[outliers, 1], '.r', alpha=0.6,
label='Outlier data')
plt.plot(line_x, line_y, '-k', label='Line model from all data')
plt.plot(line_x, line_y_robust, '-b', label='Robust line model')
plt.legend(loc='lower left')
plt.show()
doc/examples/plot_regional_maxima.py
+3 -2
View File
@@ -11,14 +11,15 @@ First we try reconstruction by dilation starting at the edges of the image. We
initialize a seed image to the minimum intensity of the image, and set its
border to be the pixel values in the original image. These maximal pixels will
get dilated in order to reconstruct the background image.
"""
import numpy as np
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt
from skimage import data
from skimage import img_as_float
from skimage.morphology import reconstruction
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt
# Convert to float: Important for subtraction later which won't work with uint8
image = img_as_float(data.coins())
doc/examples/plot_regionprops.py
+13 -27
View File
@@ -4,8 +4,8 @@ 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
@@ -13,48 +13,34 @@ 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
from skimage.transform import rotate
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')
image = np.zeros((600, 600))
rr, cc = ellipse(300, 350, 100, 220)
image[rr,cc] = 1
image = geometric_transform(image, rotate)
image = rotate(image, angle=15, order=0)
label_img = label(image)
props = regionprops(label_img, [
'BoundingBox',
'Centroid',
'Orientation',
'MajorAxisLength',
'MinorAxisLength'
])
regions = regionprops(label_img)
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']
for props in regions:
y0, x0 = props.centroid
orientation = props.orientation
x1 = x0 + math.cos(orientation) * 0.5 * props.major_axis_length
y1 = y0 - math.sin(orientation) * 0.5 * props.major_axis_length
x2 = x0 - math.sin(orientation) * 0.5 * props.minor_axis_length
y2 = y0 - math.cos(orientation) * 0.5 * props.minor_axis_length
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']
minr, minc, maxr, maxc = props.bbox
bx = (minc, maxc, maxc, minc, minc)
by = (minr, minr, maxr, maxr, minr)
plt.plot(bx, by, '-b', linewidth=2.5)
doc/examples/plot_segmentations.py
+1
View File
@@ -58,6 +58,7 @@ of Quickshift, while ``n_segments`` chooses the number of centers for kmeans.
Pascal Fua, and Sabine Suesstrunk, SLIC Superpixels Compared to
State-of-the-art Superpixel Methods, TPAMI, May 2012.
"""
from __future__ import print_function
import matplotlib.pyplot as plt
import numpy as np
doc/examples/plot_shapes.py
+15 -5
View File
@@ -4,20 +4,22 @@ Fill shapes
===========
This example shows how to fill several different shapes:
* line
* polygon
* circle
* ellipse
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage.draw import line, polygon, circle, circle_perimeter, ellipse
from skimage.draw import line, polygon, circle, circle_perimeter, \
ellipse, ellipse_perimeter
import numpy as np
import math
img = np.zeros((500, 500, 3), 'uint8')
img = np.zeros((500, 500, 3), dtype=np.uint8)
# draw line
rr, cc = line(120, 123, 20, 400)
@@ -43,8 +45,16 @@ rr, cc = ellipse(300, 300, 100, 200, img.shape)
img[rr,cc,2] = 255
# circle
rr, cc = circle_perimeter(120, 400, 50)
rr, cc = circle_perimeter(120, 400, 15)
img[rr, cc, :] = (255, 0, 0)
# ellipses
rr, cc = ellipse_perimeter(120, 400, 60, 20, orientation=math.pi / 4.)
img[rr, cc, :] = (255, 0, 255)
rr, cc = ellipse_perimeter(120, 400, 60, 20, orientation=-math.pi / 4.)
img[rr, cc, :] = (0, 0, 255)
rr, cc = ellipse_perimeter(120, 400, 60, 20, orientation=math.pi / 2.)
img[rr, cc, :] = (255, 255, 255)
plt.imshow(img)
plt.show()
doc/examples/plot_skeleton.py
+6 -4
View File
@@ -29,10 +29,12 @@ image[-100:-10, 10:-10] = 1
image[10:-10, -100:-10] = 1
# foreground object 2
rs, cs = draw.bresenham(250, 150, 10, 280)
for i in range(10): image[rs+i, cs] = 1
rs, cs = draw.bresenham(10, 150, 250, 280)
for i in range(20): image[rs+i, cs] = 1
rs, cs = draw.line(250, 150, 10, 280)
for i in range(10):
image[rs + i, cs] = 1
rs, cs = draw.line(10, 150, 250, 280)
for i in range(20):
image[rs + i, cs] = 1
# foreground object 3
ir, ic = np.indices(image.shape)
doc/examples/plot_ssim.py
+6 -5
View File
@@ -1,4 +1,4 @@
'''
"""
===========================
Structural similarity index
===========================
@@ -20,12 +20,13 @@ but with very different mean structural similarity indices.
Transactions on Image Processing, vol. 13, no. 4, pp. 600-612,
Apr. 2004.
'''
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage import data, color, io, exposure, img_as_float
from skimage import data, 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
@@ -33,13 +34,13 @@ 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)
doc/examples/plot_swirl.py
+10 -9
View File
@@ -1,4 +1,4 @@
r"""
"""
=====
Swirl
=====
@@ -8,7 +8,7 @@ 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
@@ -19,7 +19,7 @@ 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
@@ -35,16 +35,16 @@ The corresponding call to warp is::
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)
\\theta = \\arctan(y/x)
\rho = \sqrt{(x - x_0)^2 + (y - y_0)^2},
\\rho = \sqrt{(x - x_0)^2 + (y - y_0)^2},
and then transforms them according to
@@ -56,19 +56,20 @@ and then transforms them according to
s = \mathtt{strength}
\theta' = \phi + s \, e^{-\rho / r + \theta}
\\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
transformation decays to :math:`\\approx 1/1000^{\mathsf{th}}` within the
specified radius.
"""
import matplotlib.pyplot as plt
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)
doc/examples/plot_template.py
+3 -2
View File
@@ -17,13 +17,15 @@ 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]
@@ -53,4 +55,3 @@ ax3.autoscale(False)
ax3.plot(x, y, 'o', markeredgecolor='r', markerfacecolor='none', markersize=10)
plt.show()
doc/examples/plot_view_as_blocks.py
+1 -1
View File
@@ -12,8 +12,8 @@ blocks. Then, on each block, we either pool the mean, the max or the
median value of that block. The results are displayed altogether, along
with a spline interpolation of order 3 rescaling of the original `lena`
image.
"""
"""
import numpy as np
from scipy import ndimage as ndi
from matplotlib import pyplot as plt
doc/examples/plot_watershed.py
+10 -5
View File
@@ -24,11 +24,13 @@ See Wikipedia_ for more details on the algorithm.
.. _Wikipedia: http://en.wikipedia.org/wiki/Watershed_(image_processing)
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage.morphology import watershed, is_local_maximum
from scipy import ndimage
from skimage.morphology import watershed
from skimage.feature import peak_local_max
# Generate an initial image with two overlapping circles
x, y = np.indices((80, 80))
@@ -40,9 +42,9 @@ image = np.logical_or(mask_circle1, mask_circle2)
# Now we want to separate the two objects in image
# Generate the markers as local maxima of the distance to the background
from scipy import ndimage
distance = ndimage.distance_transform_edt(image)
local_maxi = is_local_maximum(distance, image, np.ones((3, 3)))
local_maxi = peak_local_max(distance, indices=False, footprint=np.ones((3, 3)),
labels=image)
markers = ndimage.label(local_maxi)[0]
labels = watershed(-distance, markers, mask=image)
@@ -50,8 +52,11 @@ fig, axes = plt.subplots(ncols=3, figsize=(8, 2.7))
ax0, ax1, ax2 = axes
ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
ax0.set_title('Overlapping objects')
ax1.imshow(-distance, cmap=plt.cm.jet, interpolation='nearest')
ax1.set_title('Distances')
ax2.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest')
ax2.set_title('Separated objects')
for ax in axes:
ax.axis('off')
doc/ext/docscrape.py
+2 -1
View File
@@ -9,6 +9,7 @@ import pydoc
from StringIO import StringIO
from warnings import warn
class Reader(object):
"""A line-based string reader.
@@ -369,7 +370,7 @@ class NumpyDocString(object):
idx = self['index']
out = []
out += ['.. index:: %s' % idx.get('default','')]
for section, references in idx.iteritems():
for section, references in idx.items():
if section == 'default':
continue
out += [' :%s: %s' % (section, ', '.join(references))]
doc/ext/docscrape_sphinx.py
+2 -1
View File
@@ -2,6 +2,7 @@ import re, inspect, textwrap, pydoc
import sphinx
from docscrape import NumpyDocString, FunctionDoc, ClassDoc
class SphinxDocString(NumpyDocString):
def __init__(self, docstring, config={}):
self.use_plots = config.get('use_plots', False)
@@ -127,7 +128,7 @@ class SphinxDocString(NumpyDocString):
return out
out += ['.. index:: %s' % idx.get('default','')]
for section, references in idx.iteritems():
for section, references in idx.items():
if section == 'default':
continue
elif section == 'refguide':
doc/ext/plot2rst.py
+12 -12
View File
@@ -187,7 +187,7 @@ def generate_example_galleries(app):
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
print("No example directory found at", example_dir)
return
rst_dir.makedirs()
@@ -225,12 +225,12 @@ def write_gallery(gallery_index, src_dir, rst_dir, cfg, depth=0):
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'
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*'_'
print('Skipping this directory')
print(80*'_')
return
gallery_description = file(gallery_template).read()
@@ -252,11 +252,11 @@ def write_gallery(gallery_index, src_dir, rst_dir, cfg, depth=0):
try:
write_example(src_name, src_dir, rst_dir, cfg)
except Exception:
print "Exception raised while running:"
print "%s in %s" % (src_name, src_dir)
print '~' * 60
print("Exception raised while running:")
print("%s in %s" % (src_name, src_dir))
print('~' * 60)
traceback.print_exc()
print '~' * 60
print('~' * 60)
continue
link_name = sub_dir.pjoin(src_name)
@@ -354,8 +354,8 @@ def write_example(src_name, src_dir, rst_dir, cfg):
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."
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)
doc/ext/plot_directive.py
+2 -1
View File
@@ -132,6 +132,7 @@ except ImportError:
def format_template(template, **kw):
return jinja.from_string(template, **kw)
import matplotlib
import matplotlib.cbook as cbook
matplotlib.use('Agg')
@@ -234,7 +235,7 @@ def mark_plot_labels(app, document):
the "htmlonly" (or "latexonly") node to the actual figure node
itself.
"""
for name, explicit in document.nametypes.iteritems():
for name, explicit in document.nametypes.items():
if not explicit:
continue
labelid = document.nameids[name]
doc/gh-pages.py
+5 -5
View File
@@ -48,7 +48,7 @@ def sh2(cmd):
out = p.communicate()[0]
retcode = p.returncode
if retcode:
print out.rstrip()
print(out.rstrip())
raise CalledProcessError(retcode, cmd)
else:
return out.rstrip()
@@ -123,12 +123,12 @@ if __name__ == '__main__':
sh('git add %s' % tag)
sh2('git commit -m"Updated doc release: %s"' % tag)
print 'Most recent commit:'
print('Most recent commit:')
sys.stdout.flush()
sh('git --no-pager log --oneline HEAD~1..')
finally:
cd(startdir)
print
print 'Now verify the build in: %r' % dest
print "If everything looks good, run 'git push' inside doc/gh-pages."
print('')
print('Now verify the build in: %r' % dest)
print("If everything looks good, run 'git push' inside doc/gh-pages.")
doc/logo/scikit_image_logo.py
+76 -137
View File
@@ -1,80 +1,57 @@
"""
Script to draw skimage logo using Scipy logo as stencil. The easiest
starting point is the `plot_colorized_logo`; the "if-main" demonstrates its use.
starting point is the `plot_colorized_logo`.
Original snake image from pixabay [1]_
.. [1] http://pixabay.com/en/snake-green-toxic-close-yellow-3237/
"""
import numpy as np
import sys
if len(sys.argv) != 2 or sys.argv[1] != '--no-plot':
print("Run with '--no-plot' flag to generate logo silently.")
else:
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
import scipy.misc
import numpy as np
import skimage.io as sio
import skimage.filter as imfilt
from skimage import img_as_float
from skimage.color import gray2rgb, rgb2gray
from skimage.exposure import rescale_intensity
from skimage.filter import sobel
import scipy_logo
# Utility functions
# =================
def get_edges(img):
edge = np.empty(img.shape)
if len(img.shape) == 3:
for i in range(3):
edge[:, :, i] = imfilt.sobel(img[:, :, i])
else:
edge = imfilt.sobel(img)
edge = rescale_intensity(edge)
return edge
def colorize(image, color, whiten=False):
"""Return colorized image from gray scale image.
def rescale_intensity(img):
i_range = float(img.max() - img.min())
img = (img - img.min()) / i_range * 255
return np.uint8(img)
def colorize(img, color, whiten=False):
"""Return colorized image from gray scale image
Parameters
----------
img : N x M array
grayscale image
color : length-3 sequence of floats
RGB color spec. Float values should be between 0 and 1.
whiten : bool
If True, a color value less than 1 increases the image intensity.
The colorized image has values from ranging between black at the lowest
intensity to `color` at the highest. If `whiten=True`, then the color
ranges from `color` to white.
"""
color = np.asarray(color)[np.newaxis, np.newaxis, :]
img = img[:, :, np.newaxis]
image = image[:, :, np.newaxis]
if whiten:
# truncate and stretch intensity range to enhance contrast
img = np.clip(img, 80, 255)
img = rescale_intensity(img)
return np.uint8(color * (255 - img) + img)
image = rescale_intensity(image, in_range=(0.3, 1))
return color * (1 - image) + image
else:
return np.uint8(img * color)
return image * color
def prepare_axes(ax):
plt.sca(ax)
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
for spine in ax.spines.itervalues():
for spine in ax.spines.values():
spine.set_visible(False)
_rgb_stack = np.ones((1, 1, 3), dtype=bool)
def gray2rgb(arr):
"""Return RGB image from a grayscale image.
Expand h x w image to h x w x 3 image where color channels are simply copies
of the grayscale image.
"""
return arr[:, :, np.newaxis] * _rgb_stack
# Logo generating classes
# =======================
@@ -82,21 +59,17 @@ class LogoBase(object):
def __init__(self):
self.logo = scipy_logo.ScipyLogo(radius=self.radius)
self.mask_1 = self.logo.get_mask(self.img.shape, 'upper left')
self.mask_2 = self.logo.get_mask(self.img.shape, 'lower right')
self.edges = get_edges(self.img)
self.mask_1 = self.logo.get_mask(self.image.shape, 'upper left')
self.mask_2 = self.logo.get_mask(self.image.shape, 'lower right')
edges = np.array([sobel(img) for img in self.image.T]).T
# truncate and stretch intensity range to enhance contrast
self.edges = np.clip(self.edges, 0, 100)
self.edges = rescale_intensity(self.edges)
self.edges = rescale_intensity(edges, in_range=(0, 0.4))
def _crop_image(self, img):
def _crop_image(self, image):
w = 2 * self.radius
x, y = self.origin
return img[y:y+w, x:x+w]
def get_canvas(self):
return 255 * np.ones(self.img.shape, dtype=np.uint8)
return image[y:y + w, x:x + w]
def plot_curve(self, **kwargs):
self.logo.plot_snake_curve(**kwargs)
@@ -104,15 +77,13 @@ class LogoBase(object):
class SnakeLogo(LogoBase):
def __init__(self):
self.radius = 250
self.origin = (420, 0)
img = sio.imread('data/snake_pixabay.jpg')
img = self._crop_image(img)
radius = 250
origin = (420, 0)
img = img.astype(float) * 1.1
img[img > 255] = 255
self.img = img.astype(np.uint8)
def __init__(self):
image = sio.imread('data/snake_pixabay.jpg')
image = self._crop_image(image)
self.image = img_as_float(image)
LogoBase.__init__(self)
@@ -120,107 +91,75 @@ class SnakeLogo(LogoBase):
snake_color = SnakeLogo()
snake = SnakeLogo()
# turn RGB image into gray image
snake.img = np.mean(snake.img, axis=2)
snake.edges = np.mean(snake.edges, axis=2)
snake.image = rgb2gray(snake.image)
snake.edges = rgb2gray(snake.edges)
# Demo plotting functions
# =======================
def plot_colorized_logo(logo, color, edges='light', switch=False, whiten=False):
"""Convenience function to plot artificially colored logo.
def plot_colorized_logo(logo, color, edges='light', whiten=False):
"""Convenience function to plot artificially-colored logo.
The upper-left half of the logo is an edge filtered image, while the
lower-right half is unfiltered.
Parameters
----------
logo : subclass of LogoBase
color : length-3 sequence of floats
logo : LogoBase instance
color : length-3 sequence of floats or 2 length-3 sequences
RGB color spec. Float values should be between 0 and 1.
edges : {'light'|'dark'}
Specifies whether Sobel edges are drawn light or dark
switch : bool
If False, the image is drawn on the southeast half of the Scipy curve
and the edge image is drawn on northwest half.
whiten : bool
whiten : bool or 2 bools
If True, a color value less than 1 increases the image intensity.
"""
if not hasattr(color[0], '__iter__'):
color = [color] * 2
color = [color] * 2 # use same color for upper-left & lower-right
if not hasattr(whiten, '__iter__'):
whiten = [whiten] * 2
img = gray2rgb(logo.get_canvas())
whiten = [whiten] * 2 # use same setting for upper-left & lower-right
image = gray2rgb(np.ones_like(logo.image))
mask_img = gray2rgb(logo.mask_2)
mask_edge = gray2rgb(logo.mask_1)
if switch:
mask_img, mask_edge = mask_edge, mask_img
# Compose image with colorized image and edge-image.
if edges == 'dark':
lg_edge = colorize(255 - logo.edges, color[0], whiten=whiten[0])
logo_edge = colorize(1 - logo.edges, color[0], whiten=whiten[0])
else:
lg_edge = colorize(logo.edges, color[0], whiten=whiten[0])
lg_img = colorize(logo.img, color[1], whiten=whiten[1])
img[mask_img] = lg_img[mask_img]
img[mask_edge] = lg_edge[mask_edge]
logo.plot_curve(lw=5, color='w')
plt.imshow(img)
logo_edge = colorize(logo.edges, color[0], whiten=whiten[0])
logo_img = colorize(logo.image, color[1], whiten=whiten[1])
image[mask_img] = logo_img[mask_img]
image[mask_edge] = logo_edge[mask_edge]
def red_light_edges(logo, **kwargs):
plot_colorized_logo(logo, (1, 0, 0), edges='light', **kwargs)
def red_dark_edges(logo, **kwargs):
plot_colorized_logo(logo, (1, 0, 0), edges='dark', **kwargs)
def blue_light_edges(logo, **kwargs):
plot_colorized_logo(logo, (0.35, 0.55, 0.85), edges='light', **kwargs)
def blue_dark_edges(logo, **kwargs):
plot_colorized_logo(logo, (0.35, 0.55, 0.85), edges='dark', **kwargs)
def green_orange_light_edges(logo, **kwargs):
colors = ((0.6, 0.8, 0.3), (1, 0.5, 0.1))
plot_colorized_logo(logo, colors, edges='light', **kwargs)
def green_orange_dark_edges(logo, **kwargs):
colors = ((0.6, 0.8, 0.3), (1, 0.5, 0.1))
plot_colorized_logo(logo, colors, edges='dark', **kwargs)
logo.plot_curve(lw=5, color='w') # plot snake curve on current axes
plt.imshow(image)
if __name__ == '__main__':
import sys
plot = False
if len(sys.argv) < 2 or sys.argv[1] != '--no-plot':
plot = True
print "Run with '--no-plot' flag to generate logo silently."
# Colors to use for the logo:
red = (1, 0, 0)
blue = (0.35, 0.55, 0.85)
green_orange = ((0.6, 0.8, 0.3), (1, 0.5, 0.1))
def plot_all():
plotters = (red_light_edges, red_dark_edges,
blue_light_edges, blue_dark_edges,
green_orange_light_edges, green_orange_dark_edges)
f, axes_array = plt.subplots(nrows=2, ncols=len(plotters))
for plot, ax_col in zip(plotters, axes_array.T):
prepare_axes(ax_col[0])
plot(snake)
prepare_axes(ax_col[1])
plot(snake, whiten=True)
color_list = [red, blue, green_orange]
edge_list = ['light', 'dark']
f, axes = plt.subplots(nrows=len(edge_list), ncols=len(color_list))
for axes_row, edges in zip(axes, edge_list):
for ax, color in zip(axes_row, color_list):
prepare_axes(ax)
plot_colorized_logo(snake, color, edges=edges)
plt.tight_layout()
def plot_snake():
def plot_official_logo():
f, ax = plt.subplots()
prepare_axes(ax)
green_orange_dark_edges(snake, whiten=(False, True))
plot_colorized_logo(snake, green_orange, edges='dark',
whiten=(False, True))
plt.savefig('green_orange_snake.png', bbox_inches='tight')
if plot:
plot_all()
plot_snake()
if plot:
plt.show()
plot_all()
plot_official_logo()
plt.show()
doc/logo/scipy_logo.py
+6 -10
View File
@@ -3,10 +3,11 @@ Code used to trace Scipy logo.
"""
import numpy as np
import matplotlib.pyplot as plt
import skimage.io as imgio
from scipy.misc import lena
import matplotlib.nxutils as nx
from skimage import io
from skimage import data
class SymmetricAnchorPoint(object):
"""Anchor point in a parametric curve with symmetric handles
@@ -185,7 +186,7 @@ class ScipyLogo(object):
def plot_image(self, **kwargs):
ax = kwargs.pop('ax', plt.gca())
img = imgio.imread('data/scipy.png')
img = io.imread('data/scipy.png')
ax.imshow(img, **kwargs)
def get_mask(self, shape, region):
@@ -236,9 +237,7 @@ def plot_snake_overlay():
logo = ScipyLogo((670, 250), 250)
logo.plot_snake_curve()
logo.plot_circle()
img = imgio.imread('data/snake_pixabay.jpg')
#mask = logo.get_mask(img.shape, 'upper left')
#img[mask] = 255
img = io.imread('data/snake_pixabay.jpg')
plt.imshow(img)
@@ -247,9 +246,7 @@ def plot_lena_overlay():
logo = ScipyLogo((300, 300), 180)
logo.plot_snake_curve()
logo.plot_circle()
img = lena()
#mask = logo.get_mask(img.shape, 'upper left')
#img[mask] = 255
img = data.lena()
plt.imshow(img)
@@ -259,4 +256,3 @@ if __name__ == '__main__':
plot_lena_overlay()
plt.show()
doc/make.bat
+9
View File
@@ -27,6 +27,14 @@ if "%1" == "help" (
goto end
)
for %%x in (html htmlhelp latex qthelp) do (
if "%1" == "%%x" (
md source\api 2>NUL
python tools/build_modref_templates.py
)
)
if "%1" == "clean" (
for /d %%i in (build\*) do rmdir /q /s %%i
del /q /s build\*
@@ -34,6 +42,7 @@ if "%1" == "clean" (
)
if "%1" == "html" (
cd source && python random_gallery.py && python coverage_generator.py && cd ..
%SPHINXBUILD% -b html %ALLSPHINXOPTS% build/html
echo.
echo.Build finished. The HTML pages are in build/html.
doc/source/_static/docversions.js
+10 -6
View File
@@ -1,17 +1,21 @@
function insert_version_links() {
var labels = ['dev', '0.8.0', '0.7.0', '0.6', '0.5', '0.4', '0.3'];
var versions = ['dev', '0.8.0', '0.7.0', '0.6', '0.5', '0.4', '0.3'];
for (i = 0; i < labels.length; i++){
function insert_version_links() {
for (i = 0; i < versions.length; i++){
open_list = '<li>'
if (typeof(DOCUMENTATION_OPTIONS) !== 'undefined') {
if ((DOCUMENTATION_OPTIONS['VERSION'] == labels[i]) ||
if ((DOCUMENTATION_OPTIONS['VERSION'] == versions[i]) ||
(DOCUMENTATION_OPTIONS['VERSION'].match(/dev$/) && (i == 0))) {
open_list = '<li id="current">'
}
}
document.write(open_list);
document.write('<a href="URL">skimage VERSION</a> </li>\n'
.replace('VERSION', labels[i])
.replace('URL', 'http://scikit-image.org/docs/' + labels[i]));
.replace('VERSION', versions[i])
.replace('URL', 'http://scikit-image.org/docs/' + versions[i]));
}
}
function stable_version() {
return versions[1];
}
doc/source/api_changes.txt
+10 -3
View File
@@ -1,9 +1,16 @@
Version 0.9
-----------
- No longer wrap ``imread`` output in an ``Image`` class
- Change default value of `sigma` parameter in ``skimage.segmentation.slic``
to 0
Version 0.4
-----------
- Switch mask and radius arguments for ``median_filter``
Version 0.3
-----------
- Remove ``as_grey``, ``dtype`` keyword from ImageCollection
- Remove ``dtype`` from imread
- Generalise ImageCollection to accept a load_func
Version 0.4
-----------
- Switch mask and radius arguments for median_filter
doc/source/conf.py
+45 -9
View File
@@ -26,9 +26,26 @@ sys.path.append(os.path.join(curpath, '..', 'ext'))
# 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', 'plot_directive', 'plot2rst',
'sphinx.ext.autosummary', 'plot2rst',
'sphinx.ext.intersphinx']
# Determine if the matplotlib has a recent enough version of the
# plot_directive, otherwise use the local fork.
try:
from matplotlib.sphinxext import plot_directive
except ImportError:
use_matplotlib_plot_directive = False
else:
try:
use_matplotlib_plot_directive = (plot_directive.__version__ >= 2)
except AttributeError:
use_matplotlib_plot_directive = False
if use_matplotlib_plot_directive:
extensions.append('matplotlib.sphinxext.plot_directive')
else:
extensions.append('plot_directive')
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
@@ -42,8 +59,8 @@ source_suffix = '.txt'
master_doc = 'index'
# General information about the project.
project = u'skimage'
copyright = u'2011, the scikit-image team'
project = 'skimage'
copyright = '2013, the scikit-image team'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
@@ -185,13 +202,13 @@ htmlhelp_basename = 'scikitimagedoc'
#latex_paper_size = 'letter'
# The font size ('10pt', '11pt' or '12pt').
#latex_font_size = '10pt'
latex_font_size = '10pt'
# Grouping the document tree into LaTeX files. List of tuples
# (source start file, target name, title, author, documentclass [howto/manual]).
latex_documents = [
('contents', 'scikitimage.tex', u'The Image Scikit Documentation',
u'SciPy Developers', 'manual'),
('contents', 'scikit-image.tex', u'The scikit-image Documentation',
u'scikit-image development team', 'manual'),
]
# The name of an image file (relative to this directory) to place at the top of
@@ -203,13 +220,32 @@ latex_documents = [
#latex_use_parts = False
# Additional stuff for the LaTeX preamble.
#latex_preamble = ''
latex_preamble = r'''
\usepackage{enumitem}
\setlistdepth{100}
\usepackage{amsmath}
\DeclareUnicodeCharacter{00A0}{\nobreakspace}
% In the parameters section, place a newline after the Parameters header
\usepackage{expdlist}
\let\latexdescription=\description
\def\description{\latexdescription{}{} \breaklabel}
% Make Examples/etc section headers smaller and more compact
\makeatletter
\titleformat{\paragraph}{\normalsize\py@HeaderFamily}%
{\py@TitleColor}{0em}{\py@TitleColor}{\py@NormalColor}
\titlespacing*{\paragraph}{0pt}{1ex}{0pt}
\makeatother
'''
# Documents to append as an appendix to all manuals.
#latex_appendices = []
# If false, no module index is generated.
#latex_use_modindex = True
latex_use_modindex = False
# -----------------------------------------------------------------------------
# Numpy extensions
@@ -243,7 +279,7 @@ matplotlib.rcParams.update({
"""
plot_include_source = True
plot_formats = [('png', 100)]
plot_formats = [('png', 100), ('pdf', 100)]
plot2rst_index_name = 'README'
plot2rst_rcparams = {'image.cmap' : 'gray',
doc/source/contribute.txt
+1 -1
View File
@@ -1,2 +1,2 @@
.. include:: ../../TASKS.txt
.. include:: ../../DEVELOPMENT.txt
.. include:: ../../CONTRIBUTING.txt
doc/source/coverage_generator.py
+29 -30
View File
@@ -49,13 +49,13 @@ def calculate_coverage(reader):
partial_items,
total_items - (partial_items + done_items + na_items),
na_items)
return list(i / total_items for i in counts)
def read_table_titles(reader):
r"""Create a dictionary with keys as section names and values as a list of
table names
return (dict)
"""
section_titles = []
@@ -69,39 +69,39 @@ def read_table_titles(reader):
# Extract names of the tables
for name in row[1:]:
if len(name) > 0:
names.append(name)
names.append(name)
else:
break
section_titles.append(row[0])
table_names[row[0]] = names
except csv.Error, e:
sys.exit('line %d: %s' % (reader.line_num, e))
return section_titles,table_names
def table_seperator(stream,lengths,character="-"):
r"""Write out table row seperator
:Input:
- *stream* (io/stream) Stream where output is put
- *lengths* (list) A list of the lengths of the columns
- *character* (string) Character to be filled between +, defaults to "-".
"""
stream.write("+")
stream.write('+'.join([character*(length+2) for length in lengths]))
stream.write("+")
def table_row(stream,data,lengths,num_columns=None):
r"""Write out table row data
:Input:
- *stream* (io/stream) Stream where output is put
- *data* (list) List of strings containing data
- *lengths* (list) A list of the lengths of the columns
- *num_columns* (string) Number of columns, defaults to the length of the
- *num_columns* (string) Number of columns, defaults to the length of the
data array
"""
if num_columns is None:
num_columns = len(data)
@@ -115,11 +115,11 @@ def table_row(stream,data,lengths,num_columns=None):
else:
entry = MISSING_STRING
stream.write(" " + entry + " "*(lengths[i] - len(entry)) + " |")
def generate_table(reader,stream,table_name=None,
column_titles=["Functionality","Matlab","Scipy","Scipy"]):
r"""Generate a reST grid table based on the CSV data in reader
Reads CSV data from *reader* until an empty line is found and generates a
reST table based on the data into *stream*. A table name can be given for
a section and table label. All rows are read in and checked for maximum
@@ -127,13 +127,13 @@ def generate_table(reader,stream,table_name=None,
widths so that the table can be constructed. If a row contains less than
the maximum number of columns a string is inserted that defaults to the
string *MISSING_STRING* which is a global parameter.
:Input:
- reader (csv.reader) The CSV reader to read in from
- stream (iostream) Output target
- table_name (string) Optional name of table, defaults to *None*
- column_titles (list) List of column titles
"""
# Find number of columns and column widths, base number of columns is
# determined by the headers
@@ -141,7 +141,6 @@ def generate_table(reader,stream,table_name=None,
data = [column_titles]
try:
for row in reader:
# print row
if len(row[0]) == 0:
break
data.append([entry.expandtabs() for entry in row])
@@ -153,7 +152,7 @@ def generate_table(reader,stream,table_name=None,
for row in data:
for i in xrange(len(row)):
column_lengths[i] = max(column_lengths[i],len(row[i]))
# Output table header
stream.write(table_name + "\n")
if table_name is not None:
@@ -167,7 +166,7 @@ def generate_table(reader,stream,table_name=None,
stream.write("\n")
table_seperator(stream,column_lengths,character="=")
stream.write("\n")
# Output table data
for row in data[1:]:
table_row(stream,row,column_lengths,num_columns)
@@ -175,28 +174,28 @@ def generate_table(reader,stream,table_name=None,
table_seperator(stream,column_lengths,character='-')
stream.write("\n")
stream.write("\n\n")
def generate_page(csv_path,stream,page_title="Coverage Tables"):
r"""Generate coverage table page
Generates all reST for all tables contained in the CSV file at *csv_path*
and output it to *stream*.
:Input:
- *csv_path* (path) Path to CSV file
- *stream* (iostream) Output stream
- *page_title* (string) Optional page title, defaults to
- *page_title* (string) Optional page title, defaults to
``Coverage Tables``.
"""
# Open reader
csv_file = open(csv_path,'U')
# Sniffer does not seem to work all the time even when an Excel
# spread sheet is being used
# dialect = csv.Sniffer().sniff(csv_file.read(1024))
# csv_file.seek(0)
# reader = csv.reader(csv_file, dialect)
reader = csv.reader(csv_file)
item_counts = calculate_coverage(reader)
csv_file.seek(0)
@@ -254,21 +253,21 @@ if __name__ == "__main__":
output_path = './coverage_table.txt'
if len(sys.argv) > 1:
if sys.argv[1][:5].lower() == "help":
print "Coverage Table Generator: coverage_generator.py"
print " Usage: coverage_generator.py [csv] [output]"
print " csv - Path to csv file, defaults to ./coverage.csv"
print " output - Ouput path, defaults to ./coverage_table.txt"
print ""
print("Coverage Table Generator: coverage_generator.py")
print(" Usage: coverage_generator.py [csv] [output]")
print(" csv - Path to csv file, defaults to ./coverage.csv")
print(" output - Ouput path, defaults to ./coverage_table.txt")
print('')
sys.exit(0)
if len(sys.argv) == 2:
csv_path = os.path.abspath(sys.argv[1])
if len(sys.argv) == 3:
output_path = os.path.abspath(sys.argv[2])
output = open(output_path,'w')
generate_page(csv_path,output)
output.close()
print("Generated %s from %s." % (output_path,csv_path))
doc/source/install.txt
+37 -15
View File
@@ -6,8 +6,15 @@ Pre-built installation
are kindly provided by Christoph Gohlke.
The latest stable release is also included as part of the `Enthought Python
Distribution (EPD) <http://enthought.com/products/epd.php>`__ and `Python(x,y)
<http://code.google.com/p/pythonxy/wiki/Welcome>`__.
Distribution (EPD) <http://enthought.com/products/epd.php>`__, `Python(x,y)
<http://code.google.com/p/pythonxy/wiki/Welcome>`__ and
`Anaconda <https://store.continuum.io/cshop/anaconda/>`__.
On Debian and Ubuntu, a Debian package ``python-skimage`` can be found in
`the Neurodebian repository <http://neuro.debian.net>`__. Follow `the
instructions <http://neuro.debian.net/#how-to-use-this-repository>`__ to
add Neurodebian to your system package manager, then look for
``python-skimage`` in the package manager.
On systems that support setuptools, the package can be installed from the
`Python packaging index <http://pypi.python.org/pypi/scikit-image>`__ using
@@ -28,36 +35,51 @@ Installation from source
Obtain the source from the git-repository at
`http://github.com/scikit-image/scikit-image
<http://github.com/scikit-image/scikit-image>`_.
by running
::
<http://github.com/scikit-image/scikit-image>`_ by running::
git clone http://github.com/scikit-image/scikit-image.git
in a terminal (You will need to have git installed on your machine).
in a terminal (you will need to have git installed on your machine).
If you do not have git installed, you can also download a zipball from
`https://github.com/scikit-image/scikit-image/zipball/master
<https://github.com/scikit-image/scikit-image/zipball/master>`_.
The SciKit can be installed globally using
::
The SciKit can be installed globally using::
python setup.py install
or locally using
::
or locally using::
python setup.py install --prefix=${HOME}
If you prefer, you can use it without installing, by simply adding
this path to your PYTHONPATH variable and compiling extensions
this path to your ``PYTHONPATH`` variable and compiling extensions
in-place::
python setup.py build_ext -i
Building with bento
-------------------
``scikit-image`` can also be built using `bento
<http://cournape.github.io/Bento/>`__. Bento depends on `WAF
<https://code.google.com/p/waf/>`__ for compilation.
Follow the `Bento installation instructions
<http://cournape.github.io/Bento/html/install.html>`__ and `download the WAF
source <http://code.google.com/p/waf/downloads/list>`__.
Tell Bento where to find WAF by setting the ``WAFDIR`` environment variable::
export WAFDIR=<path/to/waf>
From the ``scikit-image`` source directory::
bentomaker configure
bentomaker build -j # (add -i for in-place build)
bentomaker install # (when not builing in-place)
Depending on file permissions, the install commands may need to be run as sudo.
.. include:: ../../DEPENDS.txt
doc/source/plots/hough_tf.py
+7 -6
View File
@@ -1,17 +1,18 @@
import numpy as np
import matplotlib.pyplot as plt
from skimage.transform import hough
from skimage.transform import hough_line
from skimage.draw import line
img = np.zeros((100, 150), dtype=bool)
img[30, :] = 1
img[:, 65] = 1
img[35:45, 35:50] = 1
for i in range(90):
img[i, i] = 1
rr, cc = line(60, 130, 80, 10)
img[rr, cc] = 1
img += np.random.random(img.shape) > 0.95
out, angles, d = hough(img)
out, angles, d = hough_line(img)
plt.subplot(1, 2, 1)
@@ -20,8 +21,8 @@ plt.title('Input image')
plt.subplot(1, 2, 2)
plt.imshow(out, cmap=plt.cm.bone,
extent=(d[0], d[-1],
np.rad2deg(angles[0]), np.rad2deg(angles[-1])))
extent=(np.rad2deg(angles[-1]), np.rad2deg(angles[0]),
d[-1], d[0]))
plt.title('Hough transform')
plt.xlabel('Angle (degree)')
plt.ylabel('Distance (pixel)')
doc/source/themes/scikit-image/static/css/custom.css
+15
View File
@@ -41,6 +41,12 @@ h6 {
font-size: 13px;
line-height: 15px;
}
blockquote {
border-left: 0;
}
dt {
font-weight: normal;
}
.logo {
float: left;
@@ -73,6 +79,10 @@ h6 {
padding-left: 15px;
}
#current {
font-weight: bold;
}
.headerlink {
margin-left: 10px;
color: #ddd;
@@ -222,3 +232,8 @@ p.admonition-title {
width: 200px;
text-align: center !important;
}
/* misc */
div.math {
text-align: center;
}
doc/source/user_guide/data_types.txt
+2 -2
View File
@@ -34,9 +34,9 @@ violates these assumptions about the dtype range::
>>> from skimage import img_as_float
>>> image = np.arange(0, 50, 10, dtype=np.uint8)
>>> print image.astype(np.float) # These float values are out of range.
>>> print(image.astype(np.float)) # These float values are out of range.
[ 0. 10. 20. 30. 40.]
>>> print img_as_float(image)
>>> print(img_as_float(image))
[ 0. 0.03921569 0.07843137 0.11764706 0.15686275]
doc/source/user_guide/getting_help.txt
+1 -1
View File
@@ -24,7 +24,7 @@ Contributing examples to the gallery can be done on github (see
Search field
------------
The ``quick search`` field located in the sidebar of the html
The ``quick search`` field located in the navigation bar of the html
documentation can be used to search for specific keywords (segmentation,
rescaling, denoising, etc.).
doc/tools/apigen.py
+3 -3
View File
@@ -176,7 +176,7 @@ class ApiDocWriter(object):
''' Parse module defined in *uri* '''
filename = self._uri2path(uri)
if filename is None:
print filename, 'erk'
print(filename, 'erk')
# nothing that we could handle here.
return ([],[])
f = open(filename, 'rt')
@@ -260,7 +260,7 @@ class ApiDocWriter(object):
# get the names of all classes and functions
functions, classes = self._parse_module_with_import(uri)
if not len(functions) and not len(classes) and DEBUG:
print 'WARNING: Empty -', uri # dbg
print('WARNING: Empty -', uri) # dbg
return ''
# Make a shorter version of the uri that omits the package name for
@@ -449,7 +449,7 @@ class ApiDocWriter(object):
relpath = (outdir + os.path.sep).replace(relative_to + os.path.sep, '')
else:
relpath = outdir
print "outdir: ", relpath
print("outdir: ", relpath)
idx = open(path,'wt')
w = idx.write
w('.. AUTO-GENERATED FILE -- DO NOT EDIT!\n\n')
doc/tools/build_modref_templates.py
+2 -2
View File
@@ -13,7 +13,7 @@ from distutils.version import LooseVersion as V
#*****************************************************************************
def abort(error):
print '*WARNING* API documentation not generated: %s'%error
print('*WARNING* API documentation not generated: %s' % error)
exit()
if __name__ == '__main__':
@@ -54,4 +54,4 @@ if __name__ == '__main__':
]
docwriter.write_api_docs(outdir)
docwriter.write_index(outdir, 'api', relative_to='source/api')
print '%d files written' % len(docwriter.written_modules)
print('%d files written' % len(docwriter.written_modules))
doc/tools/plot_pr.py
+4 -4
View File
@@ -48,7 +48,7 @@ def fetch_PRs(user='scikit-image', repo='scikit-image', state='open'):
fetch_status = 'Fetching page %(page)d (state=%(state)s)' % params + \
' from %(user)s/%(repo)s...' % config
print fetch_status
print(fetch_status)
f = urllib.urlopen(
'https://api.github.com/repos/%(user)s/%(repo)s/pulls?%(params)s' \
@@ -61,7 +61,7 @@ def fetch_PRs(user='scikit-image', repo='scikit-image', state='open'):
if 'message' in page_data and page_data['message'] == "Not Found":
page_data = []
print 'Warning: Repo not found (%(user)s/%(repo)s)' % config
print('Warning: Repo not found (%(user)s/%(repo)s)' % config)
else:
data.extend(page_data)
@@ -69,7 +69,7 @@ def fetch_PRs(user='scikit-image', repo='scikit-image', state='open'):
try:
PRs = json.loads(open(cache, 'r').read())
print 'Loaded PRs from cache...'
print('Loaded PRs from cache...')
except IOError:
PRs = fetch_PRs(user='stefanv', repo='scikits.image', state='closed')
@@ -81,7 +81,7 @@ except IOError:
cf.flush()
nrs = [pr['number'] for pr in PRs]
print 'Processing %d pull requests...' % len(nrs)
print('Processing %d pull requests...' % len(nrs))
dates = [dateutil.parser.parse(pr['created_at']) for pr in PRs]
requirements.txt
+3
View File
@@ -0,0 +1,3 @@
cython>=0.17
matplotlib>=1.0
numpy>=1.6
setup.py
+2 -5
View File
@@ -21,7 +21,7 @@ VERSION = '0.8.1'
PYTHON_VERSION = (2, 5)
DEPENDENCIES = {
'numpy': (1, 6),
'Cython': (0, 15),
'Cython': (0, 17),
}
@@ -30,10 +30,7 @@ import sys
import re
import setuptools
from numpy.distutils.core import setup
try:
from distutils.command.build_py import build_py_2to3 as build_py
except ImportError:
from distutils.command.build_py import build_py
from distutils.command.build_py import build_py
def configuration(parent_package='', top_path=None):
skimage/__init__.py
+52 -30
View File
@@ -38,8 +38,6 @@ util
Utility Functions
-----------------
get_log
Returns the ``skimage`` log. Use this to print debug output.
img_as_float
Convert an image to floating point format, with values in [0, 1].
img_as_uint
@@ -54,6 +52,7 @@ img_as_ubyte
import os.path as _osp
import imp as _imp
import functools as _functools
from skimage._shared.utils import deprecated as _deprecated
pkg_dir = _osp.abspath(_osp.dirname(__file__))
data_dir = _osp.join(pkg_dir, 'data')
@@ -62,6 +61,7 @@ try:
from .version import version as __version__
except ImportError:
__version__ = "unbuilt-dev"
del version
try:
@@ -88,6 +88,55 @@ test_verbose = _functools.partial(test, verbose=True)
test_verbose.__doc__ = test.__doc__
class _Log(Warning):
pass
class _FakeLog(object):
def __init__(self, name):
"""
Parameters
----------
name : str
Name of the log.
repeat : bool
Whether to print repeating messages more than once (False by
default).
"""
self._name = name
import warnings
warnings.simplefilter("always", _Log)
self._warnings = warnings
def _warn(self, msg, wtype):
self._warnings.warn('%s: %s' % (wtype, msg), _Log)
def debug(self, msg):
self._warn(msg, 'DEBUG')
def info(self, msg):
self._warn(msg, 'INFO')
def warning(self, msg):
self._warn(msg, 'WARNING')
warn = warning
def error(self, msg):
self._warn(msg, 'ERROR')
def critical(self, msg):
self._warn(msg, 'CRITICAL')
def addHandler(*args):
pass
def setLevel(*args):
pass
@_deprecated()
def get_log(name=None):
"""Return a console logger.
@@ -105,39 +154,12 @@ def get_log(name=None):
http://docs.python.org/library/logging.html
"""
import logging
if name is None:
name = 'skimage'
else:
name = 'skimage.' + name
log = logging.getLogger(name)
return log
return _FakeLog(name)
def _setup_log():
"""Configure root logger.
"""
import logging
import sys
formatter = logging.Formatter(
'%(name)s: %(levelname)s: %(message)s'
)
try:
handler = logging.StreamHandler(stream=sys.stdout)
except TypeError:
handler = logging.StreamHandler(strm=sys.stdout)
handler.setFormatter(formatter)
log = get_log()
log.addHandler(handler)
log.setLevel(logging.WARNING)
log.propagate = False
_setup_log()
from .util.dtype import *
skimage/_build.py
+2 -1
View File
@@ -8,7 +8,8 @@ import subprocess
try:
WindowsError
except NameError:
WindowsError = None
class WindowsError(Exception):
pass
def cython(pyx_files, working_path=''):
skimage/_shared/six.py
+423
View File
@@ -0,0 +1,423 @@
"""Utilities for writing code that runs on Python 2 and 3"""
# Copyright (c) 2010-2013 Benjamin Peterson
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
import operator
import sys
import types
__author__ = "Benjamin Peterson <benjamin@python.org>"
__version__ = "1.3.0"
# Useful for very coarse version differentiation.
PY2 = sys.version_info[0] == 2
PY3 = sys.version_info[0] == 3
if PY3:
string_types = str,
integer_types = int,
class_types = type,
text_type = str
binary_type = bytes
MAXSIZE = sys.maxsize
else:
string_types = basestring,
integer_types = (int, long)
class_types = (type, types.ClassType)
text_type = unicode
binary_type = str
if sys.platform.startswith("java"):
# Jython always uses 32 bits.
MAXSIZE = int((1 << 31) - 1)
else:
# It's possible to have sizeof(long) != sizeof(Py_ssize_t).
class X(object):
def __len__(self):
return 1 << 31
try:
len(X())
except OverflowError:
# 32-bit
MAXSIZE = int((1 << 31) - 1)
else:
# 64-bit
MAXSIZE = int((1 << 63) - 1)
del X
def _add_doc(func, doc):
"""Add documentation to a function."""
func.__doc__ = doc
def _import_module(name):
"""Import module, returning the module after the last dot."""
__import__(name)
return sys.modules[name]
class _LazyDescr(object):
def __init__(self, name):
self.name = name
def __get__(self, obj, tp):
result = self._resolve()
setattr(obj, self.name, result)
# This is a bit ugly, but it avoids running this again.
delattr(tp, self.name)
return result
class MovedModule(_LazyDescr):
def __init__(self, name, old, new=None):
super(MovedModule, self).__init__(name)
if PY3:
if new is None:
new = name
self.mod = new
else:
self.mod = old
def _resolve(self):
return _import_module(self.mod)
class MovedAttribute(_LazyDescr):
def __init__(self, name, old_mod, new_mod, old_attr=None, new_attr=None):
super(MovedAttribute, self).__init__(name)
if PY3:
if new_mod is None:
new_mod = name
self.mod = new_mod
if new_attr is None:
if old_attr is None:
new_attr = name
else:
new_attr = old_attr
self.attr = new_attr
else:
self.mod = old_mod
if old_attr is None:
old_attr = name
self.attr = old_attr
def _resolve(self):
module = _import_module(self.mod)
return getattr(module, self.attr)
class _MovedItems(types.ModuleType):
"""Lazy loading of moved objects"""
_moved_attributes = [
MovedAttribute("cStringIO", "cStringIO", "io", "StringIO"),
MovedAttribute("filter", "itertools", "builtins", "ifilter", "filter"),
MovedAttribute("input", "__builtin__", "builtins", "raw_input", "input"),
MovedAttribute("map", "itertools", "builtins", "imap", "map"),
MovedAttribute("range", "__builtin__", "builtins", "xrange", "range"),
MovedAttribute("reload_module", "__builtin__", "imp", "reload"),
MovedAttribute("reduce", "__builtin__", "functools"),
MovedAttribute("StringIO", "StringIO", "io"),
MovedAttribute("xrange", "__builtin__", "builtins", "xrange", "range"),
MovedAttribute("zip", "itertools", "builtins", "izip", "zip"),
MovedModule("builtins", "__builtin__"),
MovedModule("configparser", "ConfigParser"),
MovedModule("copyreg", "copy_reg"),
MovedModule("http_cookiejar", "cookielib", "http.cookiejar"),
MovedModule("http_cookies", "Cookie", "http.cookies"),
MovedModule("html_entities", "htmlentitydefs", "html.entities"),
MovedModule("html_parser", "HTMLParser", "html.parser"),
MovedModule("http_client", "httplib", "http.client"),
MovedModule("email_mime_multipart", "email.MIMEMultipart", "email.mime.multipart"),
MovedModule("email_mime_text", "email.MIMEText", "email.mime.text"),
MovedModule("email_mime_base", "email.MIMEBase", "email.mime.base"),
MovedModule("BaseHTTPServer", "BaseHTTPServer", "http.server"),
MovedModule("CGIHTTPServer", "CGIHTTPServer", "http.server"),
MovedModule("SimpleHTTPServer", "SimpleHTTPServer", "http.server"),
MovedModule("cPickle", "cPickle", "pickle"),
MovedModule("queue", "Queue"),
MovedModule("reprlib", "repr"),
MovedModule("socketserver", "SocketServer"),
MovedModule("tkinter", "Tkinter"),
MovedModule("tkinter_dialog", "Dialog", "tkinter.dialog"),
MovedModule("tkinter_filedialog", "FileDialog", "tkinter.filedialog"),
MovedModule("tkinter_scrolledtext", "ScrolledText", "tkinter.scrolledtext"),
MovedModule("tkinter_simpledialog", "SimpleDialog", "tkinter.simpledialog"),
MovedModule("tkinter_tix", "Tix", "tkinter.tix"),
MovedModule("tkinter_constants", "Tkconstants", "tkinter.constants"),
MovedModule("tkinter_dnd", "Tkdnd", "tkinter.dnd"),
MovedModule("tkinter_colorchooser", "tkColorChooser",
"tkinter.colorchooser"),
MovedModule("tkinter_commondialog", "tkCommonDialog",
"tkinter.commondialog"),
MovedModule("tkinter_tkfiledialog", "tkFileDialog", "tkinter.filedialog"),
MovedModule("tkinter_font", "tkFont", "tkinter.font"),
MovedModule("tkinter_messagebox", "tkMessageBox", "tkinter.messagebox"),
MovedModule("tkinter_tksimpledialog", "tkSimpleDialog",
"tkinter.simpledialog"),
MovedModule("urllib_robotparser", "robotparser", "urllib.robotparser"),
MovedModule("winreg", "_winreg"),
]
for attr in _moved_attributes:
setattr(_MovedItems, attr.name, attr)
del attr
moves = sys.modules[__name__ + ".moves"] = _MovedItems("moves")
def add_move(move):
"""Add an item to six.moves."""
setattr(_MovedItems, move.name, move)
def remove_move(name):
"""Remove item from six.moves."""
try:
delattr(_MovedItems, name)
except AttributeError:
try:
del moves.__dict__[name]
except KeyError:
raise AttributeError("no such move, %r" % (name,))
if PY3:
_meth_func = "__func__"
_meth_self = "__self__"
_func_closure = "__closure__"
_func_code = "__code__"
_func_defaults = "__defaults__"
_func_globals = "__globals__"
_iterkeys = "keys"
_itervalues = "values"
_iteritems = "items"
_iterlists = "lists"
else:
_meth_func = "im_func"
_meth_self = "im_self"
_func_closure = "func_closure"
_func_code = "func_code"
_func_defaults = "func_defaults"
_func_globals = "func_globals"
_iterkeys = "iterkeys"
_itervalues = "itervalues"
_iteritems = "iteritems"
_iterlists = "iterlists"
try:
advance_iterator = next
except NameError:
def advance_iterator(it):
return it.next()
next = advance_iterator
try:
callable = callable
except NameError:
def callable(obj):
return any("__call__" in klass.__dict__ for klass in type(obj).__mro__)
if PY3:
def get_unbound_function(unbound):
return unbound
create_bound_method = types.MethodType
Iterator = object
else:
def get_unbound_function(unbound):
return unbound.im_func
def create_bound_method(func, obj):
return types.MethodType(func, obj, obj.__class__)
class Iterator(object):
def next(self):
return type(self).__next__(self)
callable = callable
_add_doc(get_unbound_function,
"""Get the function out of a possibly unbound function""")
get_method_function = operator.attrgetter(_meth_func)
get_method_self = operator.attrgetter(_meth_self)
get_function_closure = operator.attrgetter(_func_closure)
get_function_code = operator.attrgetter(_func_code)
get_function_defaults = operator.attrgetter(_func_defaults)
get_function_globals = operator.attrgetter(_func_globals)
def iterkeys(d, **kw):
"""Return an iterator over the keys of a dictionary."""
return iter(getattr(d, _iterkeys)(**kw))
def itervalues(d, **kw):
"""Return an iterator over the values of a dictionary."""
return iter(getattr(d, _itervalues)(**kw))
def iteritems(d, **kw):
"""Return an iterator over the (key, value) pairs of a dictionary."""
return iter(getattr(d, _iteritems)(**kw))
def iterlists(d, **kw):
"""Return an iterator over the (key, [values]) pairs of a dictionary."""
return iter(getattr(d, _iterlists)(**kw))
if PY3:
def b(s):
return s.encode("latin-1")
def u(s):
return s
unichr = chr
if sys.version_info[1] <= 1:
def int2byte(i):
return bytes((i,))
else:
# This is about 2x faster than the implementation above on 3.2+
int2byte = operator.methodcaller("to_bytes", 1, "big")
byte2int = operator.itemgetter(0)
indexbytes = operator.getitem
iterbytes = iter
import io
StringIO = io.StringIO
BytesIO = io.BytesIO
else:
def b(s):
return s
def u(s):
return unicode(s, "unicode_escape")
unichr = unichr
int2byte = chr
def byte2int(bs):
return ord(bs[0])
def indexbytes(buf, i):
return ord(buf[i])
def iterbytes(buf):
return (ord(byte) for byte in buf)
import StringIO
StringIO = BytesIO = StringIO.StringIO
_add_doc(b, """Byte literal""")
_add_doc(u, """Text literal""")
if PY3:
import builtins
exec_ = getattr(builtins, "exec")
def reraise(tp, value, tb=None):
if value.__traceback__ is not tb:
raise value.with_traceback(tb)
raise value
print_ = getattr(builtins, "print")
del builtins
else:
def exec_(_code_, _globs_=None, _locs_=None):
"""Execute code in a namespace."""
if _globs_ is None:
frame = sys._getframe(1)
_globs_ = frame.f_globals
if _locs_ is None:
_locs_ = frame.f_locals
del frame
elif _locs_ is None:
_locs_ = _globs_
exec("""exec _code_ in _globs_, _locs_""")
exec_("""def reraise(tp, value, tb=None):
raise tp, value, tb
""")
def print_(*args, **kwargs):
"""The new-style print function."""
fp = kwargs.pop("file", sys.stdout)
if fp is None:
return
def write(data):
if not isinstance(data, basestring):
data = str(data)
fp.write(data)
want_unicode = False
sep = kwargs.pop("sep", None)
if sep is not None:
if isinstance(sep, unicode):
want_unicode = True
elif not isinstance(sep, str):
raise TypeError("sep must be None or a string")
end = kwargs.pop("end", None)
if end is not None:
if isinstance(end, unicode):
want_unicode = True
elif not isinstance(end, str):
raise TypeError("end must be None or a string")
if kwargs:
raise TypeError("invalid keyword arguments to print()")
if not want_unicode:
for arg in args:
if isinstance(arg, unicode):
want_unicode = True
break
if want_unicode:
newline = unicode("\n")
space = unicode(" ")
else:
newline = "\n"
space = " "
if sep is None:
sep = space
if end is None:
end = newline
for i, arg in enumerate(args):
if i:
write(sep)
write(arg)
write(end)
_add_doc(reraise, """Reraise an exception.""")
def with_metaclass(meta, *bases):
"""Create a base class with a metaclass."""
return meta("NewBase", bases, {})
skimage/_shared/transform.pxd
+2 -2
View File
@@ -1,5 +1,5 @@
cimport numpy as cnp
cdef float integrate(cnp.ndarray[float, ndim=2, mode="c"] sat,
Py_ssize_t r0, Py_ssize_t c0, Py_ssize_t r1, Py_ssize_t c1)
cdef float integrate(float[:, ::1] sat, Py_ssize_t r0, Py_ssize_t c0,
Py_ssize_t r1, Py_ssize_t c1)
skimage/_shared/transform.pyx
+2 -2
View File
@@ -5,8 +5,8 @@
cimport numpy as cnp
cdef float integrate(cnp.ndarray[float, ndim=2, mode="c"] sat,
Py_ssize_t r0, Py_ssize_t c0, Py_ssize_t r1, Py_ssize_t c1):
cdef float integrate(float[:, ::1] sat, Py_ssize_t r0, Py_ssize_t c0,
Py_ssize_t r1, Py_ssize_t c1):
"""
Using a summed area table / integral image, calculate the sum
over a given window.
skimage/_shared/utils.py
+37 -7
View File
@@ -1,8 +1,19 @@
import warnings
import functools
import sys
from . import six
__all__ = ['deprecated']
__all__ = ['deprecated', 'get_bound_method_class']
class skimage_deprecation(Warning):
"""Create our own deprecation class, since Python >= 2.7
silences deprecations by default.
"""
pass
class deprecated(object):
@@ -25,19 +36,38 @@ class deprecated(object):
def __call__(self, func):
msg = "Call to deprecated function `%s`." % func.__name__
alt_msg = ''
if self.alt_func is not None:
msg = msg + " Use `%s` instead." % self.alt_func
alt_msg = ' Use ``%s`` instead.' % self.alt_func
msg = 'Call to deprecated function ``%s``.' % func.__name__
msg += alt_msg
@functools.wraps(func)
def wrapped(*args, **kwargs):
if self.behavior == 'warn':
func_code = six.get_function_code(func)
warnings.simplefilter('always', skimage_deprecation)
warnings.warn_explicit(msg,
category=DeprecationWarning,
filename=func.func_code.co_filename,
lineno=func.func_code.co_firstlineno + 1)
category=skimage_deprecation,
filename=func_code.co_filename,
lineno=func_code.co_firstlineno + 1)
elif self.behavior == 'raise':
raise DeprecationWarning(msg)
raise skimage_deprecation(msg)
return func(*args, **kwargs)
# modify doc string to display deprecation warning
doc = '**Deprecated function**.' + alt_msg
if wrapped.__doc__ is None:
wrapped.__doc__ = doc
else:
wrapped.__doc__ = doc + '\n\n ' + wrapped.__doc__
return wrapped
def get_bound_method_class(m):
"""Return the class for a bound method.
"""
return m.im_class if sys.version < '3' else m.__self__.__class__
skimage/color/__init__.py
+107 -1
View File
@@ -1 +1,107 @@
from .colorconv import *
from .colorconv import (convert_colorspace,
guess_spatial_dimensions,
rgb2hsv,
hsv2rgb,
rgb2xyz,
xyz2rgb,
rgb2rgbcie,
rgbcie2rgb,
rgb2grey,
rgb2gray,
gray2rgb,
xyz2lab,
lab2xyz,
lab2rgb,
rgb2lab,
rgb2hed,
hed2rgb,
lab2lch,
lch2lab,
separate_stains,
combine_stains,
rgb_from_hed,
hed_from_rgb,
rgb_from_hdx,
hdx_from_rgb,
rgb_from_fgx,
fgx_from_rgb,
rgb_from_bex,
bex_from_rgb,
rgb_from_rbd,
rbd_from_rgb,
rgb_from_gdx,
gdx_from_rgb,
rgb_from_hax,
hax_from_rgb,
rgb_from_bro,
bro_from_rgb,
rgb_from_bpx,
bpx_from_rgb,
rgb_from_ahx,
ahx_from_rgb,
rgb_from_hpx,
hpx_from_rgb,
is_rgb,
is_gray)
from .colorlabel import color_dict, label2rgb
from .delta_e import (deltaE_cie76,
deltaE_ciede94,
deltaE_ciede2000,
deltaE_cmc,
)
__all__ = ['convert_colorspace',
'guess_spatial_dimensions',
'rgb2hsv',
'hsv2rgb',
'rgb2xyz',
'xyz2rgb',
'rgb2rgbcie',
'rgbcie2rgb',
'rgb2grey',
'rgb2gray',
'gray2rgb',
'xyz2lab',
'lab2xyz',
'lab2rgb',
'rgb2lab',
'rgb2hed',
'hed2rgb',
'lab2lch',
'lch2lab',
'separate_stains',
'combine_stains',
'rgb_from_hed',
'hed_from_rgb',
'rgb_from_hdx',
'hdx_from_rgb',
'rgb_from_fgx',
'fgx_from_rgb',
'rgb_from_bex',
'bex_from_rgb',
'rgb_from_rbd',
'rbd_from_rgb',
'rgb_from_gdx',
'gdx_from_rgb',
'rgb_from_hax',
'hax_from_rgb',
'rgb_from_bro',
'bro_from_rgb',
'rgb_from_bpx',
'bpx_from_rgb',
'rgb_from_ahx',
'ahx_from_rgb',
'rgb_from_hpx',
'hpx_from_rgb',
'is_rgb',
'is_gray',
'color_dict',
'label2rgb',
'deltaE_cie76',
'deltaE_ciede94',
'deltaE_ciede2000',
'deltaE_cmc',
]
skimage/color/colorconv.py
+443 -29
View File
@@ -26,10 +26,17 @@ Supported color spaces
Derived from the RGB CIE color space. Chosen such that
``x == y == z == 1/3`` at the whitepoint, and all color matching
functions are greater than zero everywhere.
* LAB CIE : Lightness, a, b
Colorspace derived from XYZ CIE that is intended to be more
perceptually uniform
* LCH CIE : Lightness, Chroma, Hue
Defined in terms of LAB CIE. C and H are the polar representation of
a and b. The polar angle C is defined to be on (0, 2*pi)
:author: Nicolas Pinto (rgb2hsv)
:author: Ralf Gommers (hsv2rgb)
:author: Travis Oliphant (XYZ and RGB CIE functions)
:author: Matt Terry (lab2lch)
:license: modified BSD
@@ -43,18 +50,44 @@ References
from __future__ import division
__all__ = ['convert_colorspace', 'rgb2hsv', 'hsv2rgb', 'rgb2xyz', 'xyz2rgb',
'rgb2rgbcie', 'rgbcie2rgb', 'rgb2grey', 'rgb2gray', 'gray2rgb',
'xyz2lab', 'lab2xyz', 'lab2rgb', 'rgb2lab', 'is_rgb', 'is_gray'
]
__docformat__ = "restructuredtext en"
import numpy as np
from scipy import linalg
from ..util import dtype
from skimage._shared.utils import deprecated
def guess_spatial_dimensions(image):
"""Make an educated guess about whether an image has a channels dimension.
Parameters
----------
image : ndarray
The input image.
Returns
-------
spatial_dims : int or None
The number of spatial dimensions of `image`. If ambiguous, the value
is `None`.
Raises
------
ValueError
If the image array has less than two or more than four dimensions.
"""
if image.ndim == 2:
return 2
if image.ndim == 3 and image.shape[-1] != 3:
return 3
if image.ndim == 3 and image.shape[-1] == 3:
return None
if image.ndim == 4 and image.shape[-1] == 3:
return 3
else:
raise ValueError("Expected 2D, 3D, or 4D array, got %iD." % image.ndim)
@deprecated()
def is_rgb(image):
"""Test whether the image is RGB or RGBA.
@@ -67,6 +100,7 @@ def is_rgb(image):
return (image.ndim == 3 and image.shape[2] in (3, 4))
@deprecated()
def is_gray(image):
"""Test whether the image is gray (i.e. has only one color band).
@@ -76,7 +110,7 @@ def is_gray(image):
Input image.
"""
return image.ndim == 2
return image.ndim in (2, 3) and not is_rgb(image)
def convert_colorspace(arr, fromspace, tospace):
@@ -133,8 +167,9 @@ def _prepare_colorarray(arr):
"""
arr = np.asanyarray(arr)
if arr.ndim != 3 or arr.shape[2] != 3:
msg = "the input array must be have a shape == (.,.,3))"
if arr.ndim not in [3, 4] or arr.shape[-1] != 3:
msg = ("the input array must be have a shape == (.., ..,[ ..,] 3)), " +
"got (" + (", ".join(map(str, arr.shape))) + ")")
raise ValueError(msg)
return dtype.img_as_float(arr)
@@ -312,6 +347,90 @@ gray_from_rgb = np.array([[0.2125, 0.7154, 0.0721],
# CIE LAB constants for Observer= 2A, Illuminant= D65
lab_ref_white = np.array([0.95047, 1., 1.08883])
# Haematoxylin-Eosin-DAB colorspace
# From original Ruifrok's paper: A. C. Ruifrok and D. A. Johnston,
# "Quantification of histochemical staining by color deconvolution.,"
# Analytical and quantitative cytology and histology / the International
# Academy of Cytology [and] American Society of Cytology, vol. 23, no. 4,
# pp. 291-9, Aug. 2001.
rgb_from_hed = np.array([[0.65, 0.70, 0.29],
[0.07, 0.99, 0.11],
[0.27, 0.57, 0.78]])
hed_from_rgb = linalg.inv(rgb_from_hed)
# Following matrices are adapted form the Java code written by G.Landini.
# The original code is available at:
# http://www.dentistry.bham.ac.uk/landinig/software/cdeconv/cdeconv.html
# Hematoxylin + DAB
rgb_from_hdx = np.array([[0.650, 0.704, 0.286],
[0.268, 0.570, 0.776],
[0.0, 0.0, 0.0]])
rgb_from_hdx[2, :] = np.cross(rgb_from_hdx[0, :], rgb_from_hdx[1, :])
hdx_from_rgb = linalg.inv(rgb_from_hdx)
# Feulgen + Light Green
rgb_from_fgx = np.array([[0.46420921, 0.83008335, 0.30827187],
[0.94705542, 0.25373821, 0.19650764],
[0.0, 0.0, 0.0]])
rgb_from_fgx[2, :] = np.cross(rgb_from_fgx[0, :], rgb_from_fgx[1, :])
fgx_from_rgb = linalg.inv(rgb_from_fgx)
# Giemsa: Methyl Blue + Eosin
rgb_from_bex = np.array([[0.834750233, 0.513556283, 0.196330403],
[0.092789, 0.954111, 0.283111],
[0.0, 0.0, 0.0]])
rgb_from_bex[2, :] = np.cross(rgb_from_bex[0, :], rgb_from_bex[1, :])
bex_from_rgb = linalg.inv(rgb_from_bex)
# FastRed + FastBlue + DAB
rgb_from_rbd = np.array([[0.21393921, 0.85112669, 0.47794022],
[0.74890292, 0.60624161, 0.26731082],
[0.268, 0.570, 0.776]])
rbd_from_rgb = linalg.inv(rgb_from_rbd)
# Methyl Green + DAB
rgb_from_gdx = np.array([[0.98003, 0.144316, 0.133146],
[0.268, 0.570, 0.776],
[0.0, 0.0, 0.0]])
rgb_from_gdx[2, :] = np.cross(rgb_from_gdx[0, :], rgb_from_gdx[1, :])
gdx_from_rgb = linalg.inv(rgb_from_gdx)
# Hematoxylin + AEC
rgb_from_hax = np.array([[0.650, 0.704, 0.286],
[0.2743, 0.6796, 0.6803],
[0.0, 0.0, 0.0]])
rgb_from_hax[2, :] = np.cross(rgb_from_hax[0, :], rgb_from_hax[1, :])
hax_from_rgb = linalg.inv(rgb_from_hax)
# Blue matrix Anilline Blue + Red matrix Azocarmine + Orange matrix Orange-G
rgb_from_bro = np.array([[0.853033, 0.508733, 0.112656],
[0.09289875, 0.8662008, 0.49098468],
[0.10732849, 0.36765403, 0.9237484]])
bro_from_rgb = linalg.inv(rgb_from_bro)
# Methyl Blue + Ponceau Fuchsin
rgb_from_bpx = np.array([[0.7995107, 0.5913521, 0.10528667],
[0.09997159, 0.73738605, 0.6680326],
[0.0, 0.0, 0.0]])
rgb_from_bpx[2, :] = np.cross(rgb_from_bpx[0, :], rgb_from_bpx[1, :])
bpx_from_rgb = linalg.inv(rgb_from_bpx)
# Alcian Blue + Hematoxylin
rgb_from_ahx = np.array([[0.874622, 0.457711, 0.158256],
[0.552556, 0.7544, 0.353744],
[0.0, 0.0, 0.0]])
rgb_from_ahx[2, :] = np.cross(rgb_from_ahx[0, :], rgb_from_ahx[1, :])
ahx_from_rgb = linalg.inv(rgb_from_ahx)
# Hematoxylin + PAS
rgb_from_hpx = np.array([[0.644211, 0.716556, 0.266844],
[0.175411, 0.972178, 0.154589],
[0.0, 0.0, 0.0]])
rgb_from_hpx[2, :] = np.cross(rgb_from_hpx[0, :], rgb_from_hpx[1, :])
hpx_from_rgb = linalg.inv(rgb_from_hpx)
#-------------------------------------------------------------
# The conversion functions that make use of the matrices above
#-------------------------------------------------------------
@@ -333,12 +452,12 @@ def _convert(matrix, arr):
The converted array.
"""
arr = _prepare_colorarray(arr)
arr = np.swapaxes(arr, 0, 2)
arr = np.swapaxes(arr, 0, -1)
oldshape = arr.shape
arr = np.reshape(arr, (3, -1))
out = np.dot(matrix, arr)
out.shape = oldshape
out = np.swapaxes(out, 2, 0)
out = np.swapaxes(out, -1, 0)
return np.ascontiguousarray(out)
@@ -393,17 +512,19 @@ def rgb2xyz(rgb):
Parameters
----------
rgb : array_like
The image in RGB format, in a 3-D array of shape (.., .., 3).
The image in RGB format, in a 3- or 4-D array of shape
(.., ..,[ ..,] 3).
Returns
-------
out : ndarray
The image in XYZ format, in a 3-D array of shape (.., .., 3).
The image in XYZ format, in a 3- or 4-D array of shape
(.., ..,[ ..,] 3).
Raises
------
ValueError
If `rgb` is not a 3-D array of shape (.., .., 3).
If `rgb` is not a 3- or 4-D array of shape (.., ..,[ ..,] 3).
Notes
-----
@@ -548,23 +669,24 @@ def gray2rgb(image):
Parameters
----------
image : array_like
Input image of shape ``(M, N)``.
Input image of shape ``(M, N [, P])``.
Returns
-------
rgb : ndarray
RGB image of shape ``(M, N, 3)``.
RGB image of shape ``(M, N, [, P], 3)``.
Raises
------
ValueError
If the input is not 2-dimensional.
If the input is not a 2- or 3-dimensional image.
"""
if is_rgb(image):
if np.squeeze(image).ndim == 3 and image.shape[2] in (3, 4):
return image
elif is_gray(image):
return np.dstack((image, image, image))
elif image.ndim != 1 and np.squeeze(image).ndim in (1, 2, 3):
image = image[..., np.newaxis]
return np.concatenate(3 * (image,), axis=-1)
else:
raise ValueError("Input image expected to be RGB, RGBA or gray.")
@@ -575,17 +697,19 @@ def xyz2lab(xyz):
Parameters
----------
xyz : array_like
The image in XYZ format, in a 3-D array of shape (.., .., 3).
The image in XYZ format, in a 3- or 4-D array of shape
(.., ..,[ ..,] 3).
Returns
-------
out : ndarray
The image in CIE-LAB format, in a 3-D array of shape (.., .., 3).
The image in CIE-LAB format, in a 3- or 4-D array of shape
(.., ..,[ ..,] 3).
Raises
------
ValueError
If `xyz` is not a 3-D array of shape (.., .., 3).
If `xyz` is not a 3-D array of shape (.., ..,[ ..,] 3).
Notes
-----
@@ -615,14 +739,14 @@ def xyz2lab(xyz):
arr[mask] = np.power(arr[mask], 1. / 3.)
arr[~mask] = 7.787 * arr[~mask] + 16. / 116.
x, y, z = arr[:, :, 0], arr[:, :, 1], arr[:, :, 2]
x, y, z = arr[..., 0], arr[..., 1], arr[..., 2]
# Vector scaling
L = (116. * y) - 16.
a = 500.0 * (x - y)
b = 200.0 * (y - z)
return np.dstack([L, a, b])
return np.concatenate([x[..., np.newaxis] for x in [L, a, b]], axis=-1)
def lab2xyz(lab):
@@ -679,17 +803,19 @@ def rgb2lab(rgb):
Parameters
----------
rgb : array_like
The image in RGB format, in a 3-D array of shape (.., .., 3).
The image in RGB format, in a 3- or 4-D array of shape
(.., ..,[ ..,] 3).
Returns
-------
out : ndarray
The image in Lab format, in a 3-D array of shape (.., .., 3).
The image in Lab format, in a 3- or 4-D array of shape
(.., ..,[ ..,] 3).
Raises
------
ValueError
If `rgb` is not a 3-D array of shape (.., .., 3).
If `rgb` is not a 3- or 4-D array of shape (.., ..,[ ..,] 3).
Notes
-----
@@ -721,3 +847,291 @@ def lab2rgb(lab):
This function uses lab2xyz and xyz2rgb.
"""
return xyz2rgb(lab2xyz(lab))
def rgb2hed(rgb):
"""RGB to Haematoxylin-Eosin-DAB (HED) color space conversion.
Parameters
----------
rgb : array_like
The image in RGB format, in a 3-D array of shape (.., .., 3).
Returns
-------
out : ndarray
The image in HED format, in a 3-D array of shape (.., .., 3).
Raises
------
ValueError
If `rgb` is not a 3-D array of shape (.., .., 3).
References
----------
.. [1] A. C. Ruifrok and D. A. Johnston, "Quantification of histochemical
staining by color deconvolution.," Analytical and quantitative
cytology and histology / the International Academy of Cytology [and]
American Society of Cytology, vol. 23, no. 4, pp. 291-9, Aug. 2001.
Examples
--------
>>> from skimage import data
>>> from skimage.color import rgb2hed
>>> ihc = data.immunohistochemistry()
>>> ihc_hed = rgb2hed(ihc)
"""
return separate_stains(rgb, hed_from_rgb)
def hed2rgb(hed):
"""Haematoxylin-Eosin-DAB (HED) to RGB color space conversion.
Parameters
----------
hed : array_like
The image in the HED color space, in a 3-D array of shape (.., .., 3).
Returns
-------
out : ndarray
The image in RGB, in a 3-D array of shape (.., .., 3).
Raises
------
ValueError
If `hed` is not a 3-D array of shape (.., .., 3).
References
----------
.. [1] A. C. Ruifrok and D. A. Johnston, "Quantification of histochemical
staining by color deconvolution.," Analytical and quantitative
cytology and histology / the International Academy of Cytology [and]
American Society of Cytology, vol. 23, no. 4, pp. 291-9, Aug. 2001.
Examples
--------
>>> from skimage import data
>>> from skimage.color import rgb2hed, hed2rgb
>>> ihc = data.immunohistochemistry()
>>> ihc_hed = rgb2hed(ihc)
>>> ihc_rgb = hed2rgb(ihc_hed)
"""
return combine_stains(hed, rgb_from_hed)
def separate_stains(rgb, conv_matrix):
"""RGB to stain color space conversion.
Parameters
----------
rgb : array_like
The image in RGB format, in a 3-D array of shape (.., .., 3).
conv_matrix: ndarray
The stain separation matrix as described by G. Landini [1]_.
Returns
-------
out : ndarray
The image in stain color space, in a 3-D array of shape (.., .., 3).
Raises
------
ValueError
If `rgb` is not a 3-D array of shape (.., .., 3).
Notes
-----
Stain separation matrices available in the ``color`` module and their
respective colorspace:
* ``hed_from_rgb``: Hematoxylin + Eosin + DAB
* ``hdx_from_rgb``: Hematoxylin + DAB
* ``fgx_from_rgb``: Feulgen + Light Green
* ``bex_from_rgb``: Giemsa stain : Methyl Blue + Eosin
* ``rbd_from_rgb``: FastRed + FastBlue + DAB
* ``gdx_from_rgb``: Methyl Green + DAB
* ``hax_from_rgb``: Hematoxylin + AEC
* ``bro_from_rgb``: Blue matrix Anilline Blue + Red matrix Azocarmine\
+ Orange matrix Orange-G
* ``bpx_from_rgb``: Methyl Blue + Ponceau Fuchsin
* ``ahx_from_rgb``: Alcian Blue + Hematoxylin
* ``hpx_from_rgb``: Hematoxylin + PAS
References
----------
.. [1] http://www.dentistry.bham.ac.uk/landinig/software/cdeconv/cdeconv.html
Examples
--------
>>> from skimage import data
>>> from skimage.color import separate_stains, hdx_from_rgb
>>> ihc = data.immunohistochemistry()
>>> ihc_hdx = separate_stains(ihc, hdx_from_rgb)
"""
rgb = dtype.img_as_float(rgb) + 2
stains = np.dot(np.reshape(-np.log(rgb), (-1, 3)), conv_matrix)
return np.reshape(stains, rgb.shape)
def combine_stains(stains, conv_matrix):
"""Stain to RGB color space conversion.
Parameters
----------
stains : array_like
The image in stain color space, in a 3-D array of shape (.., .., 3).
conv_matrix: ndarray
The stain separation matrix as described by G. Landini [1]_.
Returns
-------
out : ndarray
The image in RGB format, in a 3-D array of shape (.., .., 3).
Raises
------
ValueError
If `stains` is not a 3-D array of shape (.., .., 3).
Notes
-----
Stain combination matrices available in the ``color`` module and their
respective colorspace:
* ``rgb_from_hed``: Hematoxylin + Eosin + DAB
* ``rgb_from_hdx``: Hematoxylin + DAB
* ``rgb_from_fgx``: Feulgen + Light Green
* ``rgb_from_bex``: Giemsa stain : Methyl Blue + Eosin
* ``rgb_from_rbd``: FastRed + FastBlue + DAB
* ``rgb_from_gdx``: Methyl Green + DAB
* ``rgb_from_hax``: Hematoxylin + AEC
* ``rgb_from_bro``: Blue matrix Anilline Blue + Red matrix Azocarmine\
+ Orange matrix Orange-G
* ``rgb_from_bpx``: Methyl Blue + Ponceau Fuchsin
* ``rgb_from_ahx``: Alcian Blue + Hematoxylin
* ``rgb_from_hpx``: Hematoxylin + PAS
References
----------
.. [1] http://www.dentistry.bham.ac.uk/landinig/software/cdeconv/cdeconv.html
Examples
--------
>>> from skimage import data
>>> from skimage.color import (separate_stains, combine_stains,
... hdx_from_rgb, rgb_from_hdx)
>>> ihc = data.immunohistochemistry()
>>> ihc_hdx = separate_stains(ihc, hdx_from_rgb)
>>> ihc_rgb = combine_stains(ihc_hdx, rgb_from_hdx)
"""
from ..exposure import rescale_intensity
stains = dtype.img_as_float(stains)
logrgb2 = np.dot(-np.reshape(stains, (-1, 3)), conv_matrix)
rgb2 = np.exp(logrgb2)
return rescale_intensity(np.reshape(rgb2 - 2, stains.shape), in_range=(-1, 1))
def lab2lch(lab):
"""CIE-LAB to CIE-LCH color space conversion.
LCH is the cylindrical representation of the LAB (Cartesian) colorspace
Parameters
----------
lab : array_like
The N-D image in CIE-LAB format. The last (`N+1`th) dimension must have
at least 3 elements, corresponding to the ``L``, ``a``, and ``b`` color
channels. Subsequent elements are copied.
Returns
-------
out : ndarray
The image in LCH format, in a N-D array with same shape as input `lab`.
Raises
------
ValueError
If `lch` does not have at least 3 color channels (i.e. l, a, b).
Notes
-----
The Hue is expressed as an angle between (0, 2*pi)
Examples
--------
>>> from skimage import data
>>> from skimage.color import rgb2lab, lab2lch
>>> lena = data.lena()
>>> lena_lab = rgb2lab(lena)
>>> lena_lch = lab2lch(lena_lab)
"""
lch = _prepare_lab_array(lab)
a, b = lch[..., 1], lch[..., 2]
lch[..., 1], lch[..., 2] = _cart2polar_2pi(a, b)
return lch
def _cart2polar_2pi(x, y):
"""convert cartesian coordiantes to polar (uses non-standard theta range!)
NON-STANDARD RANGE! Maps to (0, 2*pi) rather than usual (-pi, +pi)
"""
r, t = np.hypot(x, y), np.arctan2(y, x)
t += np.where(t < 0., 2 * np.pi, 0)
return r, t
def lch2lab(lch):
"""CIE-LCH to CIE-LAB color space conversion.
LCH is the cylindrical representation of the LAB (Cartesian) colorspace
Parameters
----------
lch : array_like
The N-D image in CIE-LCH format. The last (`N+1`th) dimension must have
at least 3 elements, corresponding to the ``L``, ``a``, and ``b`` color
channels. Subsequent elements are copied.
Returns
-------
out : ndarray
The image in LAB format, with same shape as input `lch`.
Raises
------
ValueError
If `lch` does not have at least 3 color channels (i.e. l, c, h).
Examples
--------
>>> from skimage import data
>>> from skimage.color import rgb2lab, lch2lab
>>> lena = data.lena()
>>> lena_lab = rgb2lab(lena)
>>> lena_lch = lab2lch(lena_lab)
>>> lena_lab2 = lch2lab(lena_lch)
"""
lch = _prepare_lab_array(lch)
c, h = lch[..., 1], lch[..., 2]
lch[..., 1], lch[..., 2] = c * np.cos(h), c * np.sin(h)
return lch
def _prepare_lab_array(arr):
"""Ensure input for lab2lch, lch2lab are well-posed.
Arrays must be in floating point and have at least 3 elements in
last dimension. Return a new array.
"""
arr = np.asarray(arr)
shape = arr.shape
if shape[-1] < 3:
raise ValueError('Input array has less than 3 color channels')
return dtype.img_as_float(arr, force_copy=True)
skimage/color/colorlabel.py
+96
View File
@@ -0,0 +1,96 @@
import warnings
import itertools
import numpy as np
from skimage import img_as_float
from skimage._shared import six
from skimage._shared.six.moves import zip
from .colorconv import rgb2gray, gray2rgb
from . import rgb_colors
__all__ = ['color_dict', 'label2rgb', 'DEFAULT_COLORS']
DEFAULT_COLORS = ('red', 'blue', 'yellow', 'magenta', 'green',
'indigo', 'darkorange', 'cyan', 'pink', 'yellowgreen')
color_dict = rgb_colors.__dict__
def _rgb_vector(color):
"""Return RGB color as (1, 3) array.
This RGB array gets multiplied by masked regions of an RGB image, which are
partially flattened by masking (i.e. dimensions 2D + RGB -> 1D + RGB).
Parameters
----------
color : str or array
Color name in `color_dict` or RGB float values between [0, 1].
"""
if isinstance(color, six.string_types):
color = color_dict[color]
# slice to handle RGBA colors
return np.array(color[:3]).reshape(1, 3)
def label2rgb(label, image=None, colors=None, alpha=0.3,
bg_label=-1, bg_color=None, image_alpha=1):
"""Return an RGB image where color-coded labels are painted over the image.
Parameters
----------
label : array
Integer array of labels with the same shape as `image`.
image : array
Image used as underlay for labels. If the input is an RGB image, it's
converted to grayscale before coloring.
colors : list
List of colors. If the number of labels exceeds the number of colors,
then the colors are cycled.
alpha : float [0, 1]
Opacity of colorized labels. Ignored if image is `None`.
bg_label : int
Label that's treated as the background.
bg_color : str or array
Background color. Must be a name in `color_dict` or RGB float values
between [0, 1].
image_alpha : float [0, 1]
Opacity of the image.
"""
if colors is None:
colors = DEFAULT_COLORS
colors = [_rgb_vector(c) for c in colors]
if image is None:
colorized = np.zeros(label.shape + (3,), dtype=np.float64)
# Opacity doesn't make sense if no image exists.
alpha = 1
else:
if not image.shape[:2] == label.shape:
raise ValueError("`image` and `label` must be the same shape")
if image.min() < 0:
warnings.warn("Negative intensities in `image` are not supported")
image = img_as_float(rgb2gray(image))
colorized = gray2rgb(image) * image_alpha + (1 - image_alpha)
labels = list(set(label.flat))
color_cycle = itertools.cycle(colors)
if bg_label in labels:
labels.remove(bg_label)
if bg_color is not None:
labels.insert(0, bg_label)
bg_color = _rgb_vector(bg_color)
color_cycle = itertools.chain(bg_color, color_cycle)
for c, i in zip(color_cycle, labels):
mask = (label == i)
colorized[mask] = c * alpha + colorized[mask] * (1 - alpha)
return colorized
skimage/color/delta_e.py
+339
View File
@@ -0,0 +1,339 @@
"""
Functions for calculating the "distance" between colors.
Implicit in these definitions of "distance" is the notion of "Just Noticeable
Distance" (JND). This represents the distance between colors where a human can
perceive different colors. Humans are more sensitive to certain colors than
others, which different deltaE metrics correct for with varying degrees of
sophistication.
The literature often mentions 1 as the minimum distance for visual
differentiation, but more recent studies (Mahy 1994) peg JND at 2.3
The delta-E notation comes from the German word for "Sensation" (Empfindung).
Reference
---------
http://en.wikipedia.org/wiki/Color_difference
"""
from __future__ import division
import numpy as np
from skimage.color.colorconv import lab2lch, _cart2polar_2pi
def deltaE_cie76(lab1, lab2):
"""Euclidean distance between two points in Lab color space
Parameters
----------
lab1 : array_like
reference color (Lab colorspace)
lab2 : array_like
comparison color (Lab colorspace)
Returns
-------
dE : array_like
distance between colors `lab1` and `lab2`
References
----------
.. [1] http://en.wikipedia.org/wiki/Color_difference
.. [2] A. R. Robertson, "The CIE 1976 color-difference formulae,"
Color Res. Appl. 2, 7-11 (1977).
"""
lab1 = np.asarray(lab1)
lab2 = np.asarray(lab2)
L1, a1, b1 = np.rollaxis(lab1, -1)[:3]
L2, a2, b2 = np.rollaxis(lab2, -1)[:3]
return np.sqrt((L2 - L1) ** 2 + (a2 - a1) ** 2 + (b2 - b1) ** 2)
def deltaE_ciede94(lab1, lab2, kH=1, kC=1, kL=1, k1=0.045, k2=0.015):
"""Color difference according to CIEDE 94 standard
Accommodates perceptual non-uniformities through the use of application
specific scale factors (`kH`, `kC`, `kL`, `k1`, and `k2`).
Parameters
----------
lab1 : array_like
reference color (Lab colorspace)
lab2 : array_like
comparison color (Lab colorspace)
kH : float, optional
Hue scale
kC : float, optional
Chroma scale
kL : float, optional
Lightness scale
k1 : float, optional
first scale parameter
k2 : float, optional
second scale parameter
Returns
-------
dE : array_like
color difference between `lab1` and `lab2`
Notes
-----
deltaE_ciede94 is not symmetric with respect to lab1 and lab2. CIEDE94
defines the scales for the lightness, hue, and chroma in terms of the first
color. Consequently, the first color should be regarded as the "reference"
color.
`kL`, `k1`, `k2` depend on the application and default to the values
suggested for graphic arts
========== ============== ==========
Parameter Graphic Arts Textiles
========== ============== ==========
`kL` 1.000 2.000
`k1` 0.045 0.048
`k2` 0.015 0.014
========== ============== ==========
References
----------
.. [1] http://en.wikipedia.org/wiki/Color_difference
.. [2] http://www.brucelindbloom.com/index.html?Eqn_DeltaE_CIE94.html
"""
L1, C1 = np.rollaxis(lab2lch(lab1), -1)[:2]
L2, C2 = np.rollaxis(lab2lch(lab2), -1)[:2]
dL = L1 - L2
dC = C1 - C2
dH2 = get_dH2(lab1, lab2)
SL = 1
SC = 1 + k1 * C1
SH = 1 + k2 * C1
dE2 = (dL / (kL * SL)) ** 2
dE2 += (dC / (kC * SC)) ** 2
dE2 += dH2 / (kH * SH) ** 2
return np.sqrt(dE2)
def deltaE_ciede2000(lab1, lab2, kL=1, kC=1, kH=1):
"""Color difference as given by the CIEDE 2000 standard.
CIEDE 2000 is a major revision of CIDE94. The perceptual calibration is
largely based on experience with automotive paint on smooth surfaces.
Parameters
----------
lab1 : array_like
reference color (Lab colorspace)
lab2 : array_like
comparison color (Lab colorspace)
kL : float (range), optional
lightness scale factor, 1 for "acceptably close"; 2 for "imperceptible"
see deltaE_cmc
kC : float (range), optional
chroma scale factor, usually 1
kH : float (range), optional
hue scale factor, usually 1
Returns
-------
deltaE : array_like
The distance between `lab1` and `lab2`
Notes
-----
CIEDE 2000 assumes parametric weighting factors for the lightness, chroma,
and hue (`kL`, `kC`, `kH` respectively). These default to 1.
References
----------
.. [1] http://en.wikipedia.org/wiki/Color_difference
.. [2] http://www.ece.rochester.edu/~gsharma/ciede2000/ciede2000noteCRNA.pdf
(doi:10.1364/AO.33.008069)
.. [3] M. Melgosa, J. Quesada, and E. Hita, "Uniformity of some recent
color metrics tested with an accurate color-difference tolerance
dataset," Appl. Opt. 33, 8069-8077 (1994).
"""
lab1 = np.asarray(lab1)
lab2 = np.asarray(lab2)
unroll = False
if lab1.ndim == 1 and lab2.ndim == 1:
unroll = True
if lab1.ndim == 1:
lab1 = lab1[None, :]
if lab2.ndim == 1:
lab2 = lab2[None, :]
L1, a1, b1 = np.rollaxis(lab1, -1)[:3]
L2, a2, b2 = np.rollaxis(lab2, -1)[:3]
# distort `a` based on average chroma
# then convert to lch coordines from distorted `a`
# all subsequence calculations are in the new coordiantes
# (often denoted "prime" in the literature)
Cbar = 0.5 * (np.hypot(a1, b1) + np.hypot(a2, b2))
c7 = Cbar ** 7
G = 0.5 * (1 - np.sqrt(c7 / (c7 + 25 ** 7)))
scale = 1 + G
C1, h1 = _cart2polar_2pi(a1 * scale, b1)
C2, h2 = _cart2polar_2pi(a2 * scale, b2)
# recall that c, h are polar coordiantes. c==r, h==theta
# cide2000 has four terms to delta_e:
# 1) Luminance term
# 2) Hue term
# 3) Chroma term
# 4) hue Rotation term
# lightness term
Lbar = 0.5 * (L1 + L2)
tmp = (Lbar - 50) ** 2
SL = 1 + 0.015 * tmp / np.sqrt(20 + tmp)
L_term = (L2 - L1) / (kL * SL)
# chroma term
Cbar = 0.5 * (C1 + C2) # new coordiantes
SC = 1 + 0.045 * Cbar
C_term = (C2 - C1) / (kC * SC)
# hue term
h_diff = h2 - h1
h_sum = h1 + h2
CC = C1 * C2
dH = h_diff.copy()
dH[h_diff > np.pi] -= 2 * np.pi
dH[h_diff < -np.pi] += 2 * np.pi
dH[CC == 0.] = 0. # if r == 0, dtheta == 0
dH_term = 2 * np.sqrt(CC) * np.sin(dH / 2)
Hbar = h_sum.copy()
mask = np.logical_and(CC != 0., np.abs(h_diff) > np.pi)
Hbar[mask * (h_sum < 2 * np.pi)] += 2 * np.pi
Hbar[mask * (h_sum >= 2 * np.pi)] -= 2 * np.pi
Hbar[CC == 0.] *= 2
Hbar *= 0.5
T = (1 -
0.17 * np.cos(Hbar - np.deg2rad(30)) +
0.24 * np.cos(2 * Hbar) +
0.32 * np.cos(3 * Hbar + np.deg2rad(6)) -
0.20 * np.cos(4 * Hbar - np.deg2rad(63))
)
SH = 1 + 0.015 * Cbar * T
H_term = dH_term / (kH * SH)
# hue rotation
c7 = Cbar ** 7
Rc = 2 * np.sqrt(c7 / (c7 + 25 ** 7))
dtheta = np.deg2rad(30) * np.exp(-((np.rad2deg(Hbar) - 275) / 25) ** 2)
R_term = -np.sin(2 * dtheta) * Rc * C_term * H_term
# put it all together
dE2 = L_term ** 2
dE2 += C_term ** 2
dE2 += H_term ** 2
dE2 += R_term
ans = np.sqrt(dE2)
if unroll:
ans = ans[0]
return ans
def deltaE_cmc(lab1, lab2, kL=1, kC=1):
"""Color difference from the CMC l:c standard.
This color difference was developed by the Colour Measurement Committee
(CMC) of the Society of Dyers and Colourists (United Kingdom). It is
intended for use in the textile industry.
The scale factors `kL`, `kC` set the weight given to differences in
lightness and chroma relative to differences in hue. The usual values are
``kL=2``, ``kC=1`` for "acceptability" and ``kL=1``, ``kC=1`` for
"imperceptibility". Colors with ``dE > 1`` are "different" for the given
scale factors.
Parameters
----------
lab1 : array_like
reference color (Lab colorspace)
lab2 : array_like
comparison color (Lab colorspace)
Returns
-------
dE : array_like
distance between colors `lab1` and `lab2`
Notes
-----
deltaE_cmc the defines the scales for the lightness, hue, and chroma
in terms of the first color. Consequently
``deltaE_cmc(lab1, lab2) != deltaE_cmc(lab2, lab1)``
References
----------
.. [1] http://en.wikipedia.org/wiki/Color_difference
.. [2] http://www.brucelindbloom.com/index.html?Eqn_DeltaE_CIE94.html
.. [3] F. J. J. Clarke, R. McDonald, and B. Rigg, "Modification to the
JPC79 colour-difference formula," J. Soc. Dyers Colour. 100, 128-132
(1984).
"""
L1, C1, h1 = np.rollaxis(lab2lch(lab1), -1)[:3]
L2, C2, h2 = np.rollaxis(lab2lch(lab2), -1)[:3]
dC = C1 - C2
dL = L1 - L2
dH2 = get_dH2(lab1, lab2)
T = np.where(np.logical_and(np.rad2deg(h1) >= 164, np.rad2deg(h1) <= 345),
0.56 + 0.2 * np.abs(np.cos(h1 + np.deg2rad(168))),
0.36 + 0.4 * np.abs(np.cos(h1 + np.deg2rad(35)))
)
c1_4 = C1 ** 4
F = np.sqrt(c1_4 / (c1_4 + 1900))
SL = np.where(L1 < 16, 0.511, 0.040975 * L1 / (1. + 0.01765 * L1))
SC = 0.638 + 0.0638 * C1 / (1. + 0.0131 * C1)
SH = SC * (F * T + 1 - F)
dE2 = (dL / (kL * SL)) ** 2
dE2 += (dC / (kC * SC)) ** 2
dE2 += dH2 / (SH ** 2)
return np.sqrt(dE2)
def get_dH2(lab1, lab2):
"""squared hue difference term occurring in deltaE_cmc and deltaE_ciede94
Despite its name, "dH" is not a simple difference of hue values. We avoid
working directly with the hue value, since differencing angles is
troublesome. The hue term is usually written as:
c1 = sqrt(a1**2 + b1**2)
c2 = sqrt(a2**2 + b2**2)
term = (a1-a2)**2 + (b1-b2)**2 - (c1-c2)**2
dH = sqrt(term)
However, this has poor roundoff properties when a or b is dominant.
Instead, ab is a vector with elements a and b. The same dH term can be
re-written as:
|ab1-ab2|**2 - (|ab1| - |ab2|)**2
and then simplified to:
2*|ab1|*|ab2| - 2*dot(ab1, ab2)
"""
lab1 = np.asarray(lab1)
lab2 = np.asarray(lab2)
a1, b1 = np.rollaxis(lab1, -1)[1:3]
a2, b2 = np.rollaxis(lab2, -1)[1:3]
# magnitude of (a, b) is the chroma
C1 = np.hypot(a1, b1)
C2 = np.hypot(a2, b2)
term = (C1 * C2) - (a1 * a2 + b1 * b2)
return 2*term
skimage/color/rgb_colors.py
+146
View File
@@ -0,0 +1,146 @@
aliceblue = (0.941, 0.973, 1)
antiquewhite = (0.98, 0.922, 0.843)
aqua = (0, 1, 1)
aquamarine = (0.498, 1, 0.831)
azure = (0.941, 1, 1)
beige = (0.961, 0.961, 0.863)
bisque = (1, 0.894, 0.769)
black = (0, 0, 0)
blanchedalmond = (1, 0.922, 0.804)
blue = (0, 0, 1)
blueviolet = (0.541, 0.169, 0.886)
brown = (0.647, 0.165, 0.165)
burlywood = (0.871, 0.722, 0.529)
cadetblue = (0.373, 0.62, 0.627)
chartreuse = (0.498, 1, 0)
chocolate = (0.824, 0.412, 0.118)
coral = (1, 0.498, 0.314)
cornflowerblue = (0.392, 0.584, 0.929)
cornsilk = (1, 0.973, 0.863)
crimson = (0.863, 0.0784, 0.235)
cyan = (0, 1, 1)
darkblue = (0, 0, 0.545)
darkcyan = (0, 0.545, 0.545)
darkgoldenrod = (0.722, 0.525, 0.0431)
darkgray = (0.663, 0.663, 0.663)
darkgreen = (0, 0.392, 0)
darkgrey = (0.663, 0.663, 0.663)
darkkhaki = (0.741, 0.718, 0.42)
darkmagenta = (0.545, 0, 0.545)
darkolivegreen = (0.333, 0.42, 0.184)
darkorange = (1, 0.549, 0)
darkorchid = (0.6, 0.196, 0.8)
darkred = (0.545, 0, 0)
darksalmon = (0.914, 0.588, 0.478)
darkseagreen = (0.561, 0.737, 0.561)
darkslateblue = (0.282, 0.239, 0.545)
darkslategray = (0.184, 0.31, 0.31)
darkslategrey = (0.184, 0.31, 0.31)
darkturquoise = (0, 0.808, 0.82)
darkviolet = (0.58, 0, 0.827)
deeppink = (1, 0.0784, 0.576)
deepskyblue = (0, 0.749, 1)
dimgray = (0.412, 0.412, 0.412)
dimgrey = (0.412, 0.412, 0.412)
dodgerblue = (0.118, 0.565, 1)
firebrick = (0.698, 0.133, 0.133)
floralwhite = (1, 0.98, 0.941)
forestgreen = (0.133, 0.545, 0.133)
fuchsia = (1, 0, 1)
gainsboro = (0.863, 0.863, 0.863)
ghostwhite = (0.973, 0.973, 1)
gold = (1, 0.843, 0)
goldenrod = (0.855, 0.647, 0.125)
gray = (0.502, 0.502, 0.502)
green = (0, 0.502, 0)
greenyellow = (0.678, 1, 0.184)
grey = (0.502, 0.502, 0.502)
honeydew = (0.941, 1, 0.941)
hotpink = (1, 0.412, 0.706)
indianred = (0.804, 0.361, 0.361)
indigo = (0.294, 0, 0.51)
ivory = (1, 1, 0.941)
khaki = (0.941, 0.902, 0.549)
lavender = (0.902, 0.902, 0.98)
lavenderblush = (1, 0.941, 0.961)
lawngreen = (0.486, 0.988, 0)
lemonchiffon = (1, 0.98, 0.804)
lightblue = (0.678, 0.847, 0.902)
lightcoral = (0.941, 0.502, 0.502)
lightcyan = (0.878, 1, 1)
lightgoldenrodyellow = (0.98, 0.98, 0.824)
lightgray = (0.827, 0.827, 0.827)
lightgreen = (0.565, 0.933, 0.565)
lightgrey = (0.827, 0.827, 0.827)
lightpink = (1, 0.714, 0.757)
lightsalmon = (1, 0.627, 0.478)
lightseagreen = (0.125, 0.698, 0.667)
lightskyblue = (0.529, 0.808, 0.98)
lightslategray = (0.467, 0.533, 0.6)
lightslategrey = (0.467, 0.533, 0.6)
lightsteelblue = (0.69, 0.769, 0.871)
lightyellow = (1, 1, 0.878)
lime = (0, 1, 0)
limegreen = (0.196, 0.804, 0.196)
linen = (0.98, 0.941, 0.902)
magenta = (1, 0, 1)
maroon = (0.502, 0, 0)
mediumaquamarine = (0.4, 0.804, 0.667)
mediumblue = (0, 0, 0.804)
mediumorchid = (0.729, 0.333, 0.827)
mediumpurple = (0.576, 0.439, 0.859)
mediumseagreen = (0.235, 0.702, 0.443)
mediumslateblue = (0.482, 0.408, 0.933)
mediumspringgreen = (0, 0.98, 0.604)
mediumturquoise = (0.282, 0.82, 0.8)
mediumvioletred = (0.78, 0.0824, 0.522)
midnightblue = (0.098, 0.098, 0.439)
mintcream = (0.961, 1, 0.98)
mistyrose = (1, 0.894, 0.882)
moccasin = (1, 0.894, 0.71)
navajowhite = (1, 0.871, 0.678)
navy = (0, 0, 0.502)
oldlace = (0.992, 0.961, 0.902)
olive = (0.502, 0.502, 0)
olivedrab = (0.42, 0.557, 0.137)
orange = (1, 0.647, 0)
orangered = (1, 0.271, 0)
orchid = (0.855, 0.439, 0.839)
palegoldenrod = (0.933, 0.91, 0.667)
palegreen = (0.596, 0.984, 0.596)
palevioletred = (0.686, 0.933, 0.933)
papayawhip = (1, 0.937, 0.835)
peachpuff = (1, 0.855, 0.725)
peru = (0.804, 0.522, 0.247)
pink = (1, 0.753, 0.796)
plum = (0.867, 0.627, 0.867)
powderblue = (0.69, 0.878, 0.902)
purple = (0.502, 0, 0.502)
red = (1, 0, 0)
rosybrown = (0.737, 0.561, 0.561)
royalblue = (0.255, 0.412, 0.882)
saddlebrown = (0.545, 0.271, 0.0745)
salmon = (0.98, 0.502, 0.447)
sandybrown = (0.98, 0.643, 0.376)
seagreen = (0.18, 0.545, 0.341)
seashell = (1, 0.961, 0.933)
sienna = (0.627, 0.322, 0.176)
silver = (0.753, 0.753, 0.753)
skyblue = (0.529, 0.808, 0.922)
slateblue = (0.416, 0.353, 0.804)
slategray = (0.439, 0.502, 0.565)
slategrey = (0.439, 0.502, 0.565)
snow = (1, 0.98, 0.98)
springgreen = (0, 1, 0.498)
steelblue = (0.275, 0.51, 0.706)
tan = (0.824, 0.706, 0.549)
teal = (0, 0.502, 0.502)
thistle = (0.847, 0.749, 0.847)
tomato = (1, 0.388, 0.278)
turquoise = (0.251, 0.878, 0.816)
violet = (0.933, 0.51, 0.933)
wheat = (0.961, 0.871, 0.702)
white = (1, 1, 1)
whitesmoke = (0.961, 0.961, 0.961)
yellow = (1, 1, 0)
yellowgreen = (0.604, 0.804, 0.196)
skimage/color/tests/ciede2000_test_data.txt
+38
View File
@@ -0,0 +1,38 @@
# input, intermediate, and output values for CIEDE2000 dE function
# data taken from "The CIEDE2000 Color-Difference Formula: Implementation Notes, ..." http://www.ece.rochester.edu/~gsharma/ciede2000/ciede2000noteCRNA.pdf
# tab delimited data
# pair 1 L1 a1 b1 ap1 cp1 hp1 hbar1 G T SL SC SH RT dE 2 L2 a2 b2 ap2 cp2 hp2
1 1 50.0000 2.6772 -79.7751 2.6774 79.8200 271.9222 270.9611 0.0001 0.6907 1.0000 4.6578 1.8421 -1.7042 2.0425 2 50.0000 0.0000 -82.7485 0.0000 82.7485 270.0000
2 1 50.0000 3.1571 -77.2803 3.1573 77.3448 272.3395 271.1698 0.0001 0.6843 1.0000 4.6021 1.8216 -1.7070 2.8615 2 50.0000 0.0000 -82.7485 0.0000 82.7485 270.0000
3 1 50.0000 2.8361 -74.0200 2.8363 74.0743 272.1944 271.0972 0.0001 0.6865 1.0000 4.5285 1.8074 -1.7060 3.4412 2 50.0000 0.0000 -82.7485 0.0000 82.7485 270.0000
4 1 50.0000 -1.3802 -84.2814 -1.3803 84.2927 269.0618 269.5309 0.0001 0.7357 1.0000 4.7584 1.9217 -1.6809 1.0000 2 50.0000 0.0000 -82.7485 0.0000 82.7485 270.0000
5 1 50.0000 -1.1848 -84.8006 -1.1849 84.8089 269.1995 269.5997 0.0001 0.7335 1.0000 4.7700 1.9218 -1.6822 1.0000 2 50.0000 0.0000 -82.7485 0.0000 82.7485 270.0000
6 1 50.0000 -0.9009 -85.5211 -0.9009 85.5258 269.3964 269.6982 0.0001 0.7303 1.0000 4.7862 1.9217 -1.6840 1.0000 2 50.0000 0.0000 -82.7485 0.0000 82.7485 270.0000
7 1 50.0000 0.0000 0.0000 0.0000 0.0000 0.0000 126.8697 0.5000 1.2200 1.0000 1.0562 1.0229 0.0000 2.3669 2 50.0000 -1.0000 2.0000 -1.5000 2.5000 126.8697
8 1 50.0000 -1.0000 2.0000 -1.5000 2.5000 126.8697 126.8697 0.5000 1.2200 1.0000 1.0562 1.0229 0.0000 2.3669 2 50.0000 0.0000 0.0000 0.0000 0.0000 0.0000
9 1 50.0000 2.4900 -0.0010 3.7346 3.7346 359.9847 269.9854 0.4998 0.7212 1.0000 1.1681 1.0404 -0.0022 7.1792 2 50.0000 -2.4900 0.0009 -3.7346 3.7346 179.9862
10 1 50.0000 2.4900 -0.0010 3.7346 3.7346 359.9847 269.9847 0.4998 0.7212 1.0000 1.1681 1.0404 -0.0022 7.1792 2 50.0000 -2.4900 0.0010 -3.7346 3.7346 179.9847
11 1 50.0000 2.4900 -0.0010 3.7346 3.7346 359.9847 89.9839 0.4998 0.6175 1.0000 1.1681 1.0346 0.0000 7.2195 2 50.0000 -2.4900 0.0011 -3.7346 3.7346 179.9831
12 1 50.0000 2.4900 -0.0010 3.7346 3.7346 359.9847 89.9831 0.4998 0.6175 1.0000 1.1681 1.0346 0.0000 7.2195 2 50.0000 -2.4900 0.0012 -3.7346 3.7346 179.9816
13 1 50.0000 -0.0010 2.4900 -0.0015 2.4900 90.0345 180.0328 0.4998 0.9779 1.0000 1.1121 1.0365 0.0000 4.8045 2 50.0000 0.0009 -2.4900 0.0013 2.4900 270.0311
14 1 50.0000 -0.0010 2.4900 -0.0015 2.4900 90.0345 180.0345 0.4998 0.9779 1.0000 1.1121 1.0365 0.0000 4.8045 2 50.0000 0.0010 -2.4900 0.0015 2.4900 270.0345
15 1 50.0000 -0.0010 2.4900 -0.0015 2.4900 90.0345 0.0362 0.4998 1.3197 1.0000 1.1121 1.0493 0.0000 4.7461 2 50.0000 0.0011 -2.4900 0.0016 2.4900 270.0380
16 1 50.0000 2.5000 0.0000 3.7496 3.7496 0.0000 315.0000 0.4998 0.8454 1.0000 1.1406 1.0396 -0.0001 4.3065 2 50.0000 0.0000 -2.5000 0.0000 2.5000 270.0000
17 1 50.0000 2.5000 0.0000 3.4569 3.4569 0.0000 346.2470 0.3827 1.4453 1.1608 1.9547 1.4599 -0.0003 27.1492 2 73.0000 25.0000 -18.0000 34.5687 38.9743 332.4939
18 1 50.0000 2.5000 0.0000 3.4954 3.4954 0.0000 51.7766 0.3981 0.6447 1.0640 1.7498 1.1612 0.0000 22.8977 2 61.0000 -5.0000 29.0000 -6.9907 29.8307 103.5532
19 1 50.0000 2.5000 0.0000 3.5514 3.5514 0.0000 272.2362 0.4206 0.6521 1.0251 1.9455 1.2055 -0.8219 31.9030 2 56.0000 -27.0000 -3.0000 -38.3556 38.4728 184.4723
20 1 50.0000 2.5000 0.0000 3.5244 3.5244 0.0000 11.9548 0.4098 1.1031 1.0400 1.9120 1.3353 0.0000 19.4535 2 58.0000 24.0000 15.0000 33.8342 37.0102 23.9095
21 1 50.0000 2.5000 0.0000 3.7494 3.7494 0.0000 3.5056 0.4997 1.2616 1.0000 1.1923 1.0808 0.0000 1.0000 2 50.0000 3.1736 0.5854 4.7596 4.7954 7.0113
22 1 50.0000 2.5000 0.0000 3.7493 3.7493 0.0000 0.0000 0.4997 1.3202 1.0000 1.1956 1.0861 0.0000 1.0000 2 50.0000 3.2972 0.0000 4.9450 4.9450 0.0000
23 1 50.0000 2.5000 0.0000 3.7497 3.7497 0.0000 5.8190 0.4999 1.2197 1.0000 1.1486 1.0604 0.0000 1.0000 2 50.0000 1.8634 0.5757 2.7949 2.8536 11.6380
24 1 50.0000 2.5000 0.0000 3.7493 3.7493 0.0000 1.9603 0.4997 1.2883 1.0000 1.1946 1.0836 0.0000 1.0000 2 50.0000 3.2592 0.3350 4.8879 4.8994 3.9206
25 1 60.2574 -34.0099 36.2677 -34.0678 49.7590 133.2085 132.0835 0.0017 1.3010 1.1427 3.2946 1.9951 0.0000 1.2644 2 60.4626 -34.1751 39.4387 -34.2333 52.2238 130.9584
26 1 63.0109 -31.0961 -5.8663 -32.6194 33.1427 190.1951 188.8221 0.0490 0.9402 1.1831 2.4549 1.4560 0.0000 1.2630 2 62.8187 -29.7946 -4.0864 -31.2542 31.5202 187.4490
27 1 61.2901 3.7196 -5.3901 5.5668 7.7487 315.9240 310.0313 0.4966 0.6952 1.1586 1.3092 1.0717 -0.0032 1.8731 2 61.4292 2.2480 -4.9620 3.3644 5.9950 304.1385
28 1 35.0831 -44.1164 3.7933 -44.3939 44.5557 175.1161 176.4290 0.0063 1.0168 1.2148 2.9105 1.6476 0.0000 1.8645 2 35.0232 -40.0716 1.5901 -40.3237 40.3550 177.7418
29 1 22.7233 20.0904 -46.6940 20.1424 50.8532 293.3339 291.3809 0.0026 0.3636 1.4014 3.1597 1.2617 -1.2537 2.0373 2 23.0331 14.9730 -42.5619 15.0118 45.1317 289.4279
30 1 36.4612 47.8580 18.3852 47.9197 51.3256 20.9901 21.8781 0.0013 0.9239 1.1943 3.3888 1.7357 0.0000 1.4146 2 36.2715 50.5065 21.2231 50.5716 54.8444 22.7660
31 1 90.8027 -2.0831 1.4410 -3.1245 3.4408 155.2410 167.1011 0.4999 1.1546 1.6110 1.1329 1.0511 0.0000 1.4441 2 91.1528 -1.6435 0.0447 -2.4651 2.4655 178.9612
32 1 90.9257 -0.5406 -0.9208 -0.8109 1.2270 228.6315 218.4363 0.5000 1.3916 1.5930 1.0620 1.0288 0.0000 1.5381 2 88.6381 -0.8985 -0.7239 -1.3477 1.5298 208.2412
33 1 6.7747 -0.2908 -2.4247 -0.4362 2.4636 259.8025 263.0049 0.4999 0.9556 1.6517 1.1057 1.0337 -0.0004 0.6377 2 5.8714 -0.0985 -2.2286 -0.1477 2.2335 266.2073
34 1 2.0776 0.0795 -1.1350 0.1192 1.1412 275.9978 268.0910 0.5000 0.7826 1.7246 1.0383 1.0100 0.0000 0.9082 2 0.9033 -0.0636 -0.5514 -0.0954 0.5596 260.18421
skimage/color/tests/test_colorconv.py
+98 -12
View File
@@ -14,26 +14,48 @@ Authors
import os.path
import numpy as np
from numpy.testing import *
from numpy.testing import (assert_equal,
assert_almost_equal,
assert_array_almost_equal,
assert_raises,
TestCase,
)
from skimage import img_as_float
from skimage import img_as_float, img_as_ubyte
from skimage.io import imread
from skimage.color import (
rgb2hsv, hsv2rgb,
rgb2xyz, xyz2rgb,
rgb2rgbcie, rgbcie2rgb,
convert_colorspace,
rgb2grey, gray2rgb,
xyz2lab, lab2xyz,
lab2rgb, rgb2lab,
is_rgb, is_gray
)
from skimage.color import (rgb2hsv, hsv2rgb,
rgb2xyz, xyz2rgb,
rgb2hed, hed2rgb,
separate_stains,
combine_stains,
rgb2rgbcie, rgbcie2rgb,
convert_colorspace,
rgb2grey, gray2rgb,
xyz2lab, lab2xyz,
lab2rgb, rgb2lab,
is_rgb, is_gray,
lab2lch, lch2lab,
guess_spatial_dimensions
)
from skimage import data_dir, data
import colorsys
def test_guess_spatial_dimensions():
im1 = np.zeros((5, 5))
im2 = np.zeros((5, 5, 5))
im3 = np.zeros((5, 5, 3))
im4 = np.zeros((5, 5, 5, 3))
im5 = np.zeros((5,))
assert_equal(guess_spatial_dimensions(im1), 2)
assert_equal(guess_spatial_dimensions(im2), 3)
assert_equal(guess_spatial_dimensions(im3), None)
assert_equal(guess_spatial_dimensions(im4), 3)
assert_raises(ValueError, guess_spatial_dimensions, im5)
class TestColorconv(TestCase):
img_rgb = imread(os.path.join(data_dir, 'color.png'))
@@ -121,6 +143,32 @@ class TestColorconv(TestCase):
img_rgb = img_as_float(self.img_rgb)
assert_array_almost_equal(xyz2rgb(rgb2xyz(img_rgb)), img_rgb)
# RGB<->HED roundtrip with ubyte image
def test_hed_rgb_roundtrip(self):
img_rgb = img_as_ubyte(self.img_rgb)
assert_equal(img_as_ubyte(hed2rgb(rgb2hed(img_rgb))), img_rgb)
# RGB<->HED roundtrip with float image
def test_hed_rgb_float_roundtrip(self):
img_rgb = img_as_float(self.img_rgb)
assert_array_almost_equal(hed2rgb(rgb2hed(img_rgb)), img_rgb)
# RGB<->HDX roundtrip with ubyte image
def test_hdx_rgb_roundtrip(self):
from skimage.color.colorconv import hdx_from_rgb, rgb_from_hdx
img_rgb = self.img_rgb
conv = combine_stains(separate_stains(img_rgb, hdx_from_rgb),
rgb_from_hdx)
assert_equal(img_as_ubyte(conv), img_rgb)
# RGB<->HDX roundtrip with ubyte image
def test_hdx_rgb_roundtrip(self):
from skimage.color.colorconv import hdx_from_rgb, rgb_from_hdx
img_rgb = img_as_float(self.img_rgb)
conv = combine_stains(separate_stains(img_rgb, hdx_from_rgb),
rgb_from_hdx)
assert_array_almost_equal(conv, img_rgb)
# RGB to RGB CIE
def test_rgb2rgbcie_conversion(self):
gt = np.array([[[ 0.1488856 , 0.18288098, 0.19277574],
@@ -202,6 +250,43 @@ class TestColorconv(TestCase):
img_rgb = img_as_float(self.img_rgb)
assert_array_almost_equal(lab2rgb(rgb2lab(img_rgb)), img_rgb)
def test_lab_lch_roundtrip(self):
rgb = img_as_float(self.img_rgb)
lab = rgb2lab(rgb)
lab2 = lch2lab(lab2lch(lab))
assert_array_almost_equal(lab2, lab)
def test_rgb_lch_roundtrip(self):
rgb = img_as_float(self.img_rgb)
lab = rgb2lab(rgb)
lch = lab2lch(lab)
lab2 = lch2lab(lch)
rgb2 = lab2rgb(lab2)
assert_array_almost_equal(rgb, rgb2)
def test_lab_lch_0d(self):
lab0 = self._get_lab0()
lch0 = lab2lch(lab0)
lch2 = lab2lch(lab0[None, None, :])
assert_array_almost_equal(lch0, lch2[0, 0, :])
def test_lab_lch_1d(self):
lab0 = self._get_lab0()
lch0 = lab2lch(lab0)
lch1 = lab2lch(lab0[None, :])
assert_array_almost_equal(lch0, lch1[0, :])
def test_lab_lch_3d(self):
lab0 = self._get_lab0()
lch0 = lab2lch(lab0)
lch3 = lab2lch(lab0[None, None, None, :])
assert_array_almost_equal(lch0, lch3[0, 0, 0, :])
def _get_lab0(self):
rgb = img_as_float(self.img_rgb[:1, :1, :])
return rgb2lab(rgb)[0, 0, :]
def test_gray2rgb():
x = np.array([0, 0.5, 1])
assert_raises(ValueError, gray2rgb, x)
@@ -238,4 +323,5 @@ def test_is_rgb():
if __name__ == "__main__":
from numpy.testing import run_module_suite
run_module_suite()

Some files were not shown because too many files have changed in this diff Show More