mirror of
https://github.com/wassname/scikit-image.git
synced 2026-08-14 12:50:39 +08:00
Merge branch 'rank_filters' of https://github.com/odebeir/scikits-image into rank-filters
This commit is contained in:
Johannes Schönberger
+2
-1
@@ -22,4 +22,5 @@ doc/source/auto_examples/images/plot_*.png
|
||||
doc/source/auto_examples/images/thumb
|
||||
doc/source/auto_examples/applications/
|
||||
doc/source/_static/random.js
|
||||
|
||||
.idea/
|
||||
*.log
|
||||
|
||||
@@ -117,3 +117,6 @@
|
||||
|
||||
- Petter Strandmark
|
||||
Perimeter calculation in regionprops.
|
||||
|
||||
- Olivier Debeir
|
||||
Rank filters (8- and 16-bits) using sliding window.
|
||||
@@ -0,0 +1,585 @@
|
||||
"""
|
||||
===============================================================
|
||||
Rank filters
|
||||
===============================================================
|
||||
|
||||
Rank filters are non-linear filters using the local grey levels 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.
|
||||
|
||||
Rank filters can be used for several purposes such as:
|
||||
|
||||
* image quality enhancement
|
||||
e.g. image smoothing, sharpening
|
||||
|
||||
* image pre-processing
|
||||
e.g. noise reduction, contrast enhancement
|
||||
|
||||
* feature extraction
|
||||
e.g. border detection, isolated point detection
|
||||
|
||||
* post-processing
|
||||
e.g. small object removal, object grouping, contour smoothing
|
||||
|
||||
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 grey level image using some of the linear and non-linear filters
|
||||
availables in skimage. We use the ``camera`` image from ``skimage.data``.
|
||||
|
||||
.. [1] Pierre Soille, On morphological operators based on rank filters, Pattern Recognition 35 (2002) 527-535.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from skimage import data
|
||||
|
||||
ima = data.camera()
|
||||
hist = np.histogram(ima, bins=np.arange(0, 256))
|
||||
|
||||
plt.figure(figsize=(8, 3))
|
||||
plt.subplot(121)
|
||||
plt.imshow(ima, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
plt.subplot(122)
|
||||
plt.plot(hist[1][:-1], hist[0], lw=2)
|
||||
plt.title('histogram of grey values')
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Noise removal
|
||||
==============
|
||||
|
||||
some noise is added to the image, 1% of pixels are randomly set to 255, %1% are randomly set to 0.
|
||||
The **median** filter is applied to remove the noise.
|
||||
|
||||
.. note:: there is 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>.99] = 255
|
||||
nima[noise<.01] = 0
|
||||
|
||||
from skimage.filter.rank import median
|
||||
from skimage.morphology import disk
|
||||
|
||||
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.subplot(2,2,2)
|
||||
plt.imshow(lo,cmap=plt.cm.gray,vmin=0,vmax=255)
|
||||
plt.xlabel('median $r=1$')
|
||||
plt.subplot(2,2,3)
|
||||
plt.imshow(hi,cmap=plt.cm.gray,vmin=0,vmax=255)
|
||||
plt.xlabel('median $r=5$')
|
||||
plt.subplot(2,2,4)
|
||||
plt.imshow(ext,cmap=plt.cm.gray,vmin=0,vmax=255)
|
||||
plt.xlabel('median $r=20$')
|
||||
|
||||
"""
|
||||
.. 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 too such as the camera tripod.
|
||||
Median filter is commonly used for noise removal because borders are preserved.
|
||||
|
||||
Image smoothing
|
||||
================
|
||||
|
||||
The example hereunder shows how a local **mean** smooth the cameraman image.
|
||||
|
||||
"""
|
||||
from skimage.filter.rank import mean
|
||||
|
||||
fig = plt.figure(figsize=[10,7])
|
||||
|
||||
loc_mean = mean(nima,disk(10))
|
||||
plt.subplot(1,2,1)
|
||||
plt.imshow(ima,cmap=plt.cm.gray,vmin=0,vmax=255)
|
||||
plt.xlabel('original')
|
||||
plt.subplot(1,2,2)
|
||||
plt.imshow(loc_mean,cmap=plt.cm.gray,vmin=0,vmax=255)
|
||||
plt.xlabel('local mean $r=10$')
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
One may be interested in smoothing an image while preserving important borders (median filters already achieved this),
|
||||
here we use the **bilateral** filter that restrict the local neighborhood to pixel having a grey level similar to the
|
||||
central one.
|
||||
|
||||
.. note:: 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)
|
||||
|
||||
bilat = bilateral_mean(ima.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.subplot(2,2,3)
|
||||
plt.imshow(bilat,cmap=plt.cm.gray)
|
||||
plt.xlabel('bilateral mean')
|
||||
plt.subplot(2,2,2)
|
||||
plt.imshow(ima[200:350,350:450],cmap=plt.cm.gray)
|
||||
plt.subplot(2,2,4)
|
||||
plt.imshow(bilat[200:350,350:450],cmap=plt.cm.gray)
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
One can see that the large continuous part of the image (e.g.sky) are smoothed whereas other details are preserved.
|
||||
|
||||
|
||||
Contrast enhancement
|
||||
====================
|
||||
|
||||
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 emphasized every local graylevel variations.
|
||||
|
||||
.. [2] http://en.wikipedia.org/wiki/Histogram_equalization
|
||||
.. [3] http://en.wikipedia.org/wiki/Adaptive_histogram_equalization
|
||||
|
||||
"""
|
||||
|
||||
from skimage import exposure
|
||||
from skimage.filter import rank
|
||||
|
||||
ima = data.camera()
|
||||
# equalize globally and locally
|
||||
glob = exposure.equalize(ima)*255
|
||||
loc = rank.equalize(ima,disk(20))
|
||||
|
||||
# extract histogram for each image
|
||||
hist = np.histogram(ima, 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.axis('off')
|
||||
plt.subplot(322)
|
||||
plt.plot(hist[1][:-1], hist[0], lw=2)
|
||||
plt.title('histogram of grey values')
|
||||
plt.subplot(323)
|
||||
plt.imshow(glob, cmap=plt.cm.gray, 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.subplot(325)
|
||||
plt.imshow(loc, cmap=plt.cm.gray, 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')
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
an other way to maximize the number of grey level used for an image is to apply a local autoleveling,
|
||||
i.e. here a pixel grey level is proportionally remapped between local minimum and local maximum.
|
||||
|
||||
The following example show how local autolevel enhance the camaraman picture.
|
||||
"""
|
||||
from skimage.filter.rank import autolevel
|
||||
|
||||
ima = data.camera()
|
||||
selem = disk(10)
|
||||
|
||||
auto = autolevel(ima.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.subplot(1,2,2)
|
||||
plt.imshow(auto,cmap=plt.cm.gray)
|
||||
plt.xlabel('local autolevel')
|
||||
"""
|
||||
.. 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 to given percentiles (one inferior, one superior)
|
||||
in place of local minimum and maximum. The example bellow illustrate how the percentile parameters influence the
|
||||
local autolevel result.
|
||||
|
||||
"""
|
||||
from skimage.filter.rank import percentile_autolevel
|
||||
|
||||
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)
|
||||
|
||||
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')
|
||||
|
||||
ax1.imshow(np.hstack((loc_perc_autolevel0,loc_perc_autolevel1)),vmin=0,vmax=255)
|
||||
ax1.set_title('percentile autolevel 0%,1%')
|
||||
ax2.imshow(np.hstack((loc_perc_autolevel2,loc_perc_autolevel3)),vmin=0,vmax=255)
|
||||
ax2.set_title('percentile autolevel 5% and 10%')
|
||||
|
||||
for ax in axes:
|
||||
ax.axis('off')
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Morphological contrast enhancement filter replaces the central pixel by local maximum
|
||||
if the original grey level value if closest to local maximum, by the minimum local otherwise.
|
||||
|
||||
"""
|
||||
|
||||
from skimage.filter.rank import morph_contr_enh
|
||||
|
||||
ima = data.camera()
|
||||
|
||||
enh = morph_contr_enh(ima,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.subplot(2,2,3)
|
||||
plt.imshow(enh,cmap=plt.cm.gray)
|
||||
plt.xlabel('local morphlogical contrast enhancement')
|
||||
plt.subplot(2,2,2)
|
||||
plt.imshow(ima[200:350,350:450],cmap=plt.cm.gray)
|
||||
plt.subplot(2,2,4)
|
||||
plt.imshow(enh[200:350,350:450],cmap=plt.cm.gray)
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
The percentile version of the local morphological contrast enhancement, uses percentile *p0* and *p1* instead of local
|
||||
minimum and local maximum.
|
||||
|
||||
"""
|
||||
|
||||
from skimage.filter.rank import percentile_morph_contr_enh
|
||||
|
||||
ima = data.camera()
|
||||
|
||||
penh = percentile_morph_contr_enh(ima,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.subplot(2,2,3)
|
||||
plt.imshow(penh,cmap=plt.cm.gray)
|
||||
plt.xlabel('local percentile morphlogical\n contrast enhancement')
|
||||
plt.subplot(2,2,2)
|
||||
plt.imshow(ima[200:350,350:450],cmap=plt.cm.gray)
|
||||
plt.subplot(2,2,4)
|
||||
plt.imshow(penh[200:350,350:450],cmap=plt.cm.gray)
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Image morphology
|
||||
================
|
||||
|
||||
Local maximum and local minimum are the base operators for grey level morphology.
|
||||
|
||||
.. note:: ``skimage.dilate`` and ``skimage.erode`` are equivalent filters (see below for comparison).
|
||||
|
||||
Here is an example of classical morphological grey level filters : opening, closing and morphological gradient.
|
||||
|
||||
"""
|
||||
|
||||
from skimage.filter.rank import maximum,minimum,gradient
|
||||
|
||||
ima = data.camera()
|
||||
|
||||
closing = maximum(minimum(ima,disk(5)),disk(5))
|
||||
opening = minimum(maximum(ima,disk(5)),disk(5))
|
||||
grad = gradient(ima,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.subplot(2,2,2)
|
||||
plt.imshow(closing,cmap=plt.cm.gray)
|
||||
plt.xlabel('grey level closing')
|
||||
plt.subplot(2,2,3)
|
||||
plt.imshow(opening,cmap=plt.cm.gray)
|
||||
plt.xlabel('grey level opening')
|
||||
plt.subplot(2,2,4)
|
||||
plt.imshow(grad,cmap=plt.cm.gray)
|
||||
plt.xlabel('morphological gradient')
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Feature extraction
|
||||
===================
|
||||
|
||||
Local histogram 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 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.
|
||||
|
||||
.. note:: to better use the available image bit, the function returns 10x entropy for 8-bit images and 1000x entropy
|
||||
for 16-bit images.
|
||||
|
||||
"""
|
||||
|
||||
|
||||
|
||||
from skimage import data
|
||||
from skimage.filter.rank import entropy
|
||||
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
|
||||
|
||||
ent8 = entropy(a8,disk(5)) # pixel value contain 10x the local entropy
|
||||
ent16 = entropy(a16,disk(5)) # pixel value contain 1000x the local entropy
|
||||
|
||||
# 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.colorbar()
|
||||
|
||||
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.colorbar()
|
||||
plt.show()
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Implementation
|
||||
================
|
||||
|
||||
The central part of the ``skimage.rank`` filters is build on a sliding window that update local grey 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.
|
||||
|
||||
"""
|
||||
from time import time
|
||||
|
||||
from scipy.ndimage.filters import percentile_filter
|
||||
from skimage.morphology import dilation
|
||||
from skimage.filter import median_filter
|
||||
from skimage.filter.rank import median,maximum
|
||||
|
||||
def exec_and_timeit(func):
|
||||
""" Decorator that returns both function results and execution time
|
||||
(result, ms)
|
||||
"""
|
||||
def wrapper(*arg):
|
||||
t1 = time()
|
||||
res = func(*arg)
|
||||
t2 = time()
|
||||
ms = (t2-t1)*1000.0
|
||||
return (res,ms)
|
||||
return wrapper
|
||||
|
||||
|
||||
@exec_and_timeit
|
||||
def cr_med(image,selem):
|
||||
return median(image=image,selem = selem)
|
||||
|
||||
@exec_and_timeit
|
||||
def cr_max(image,selem):
|
||||
return maximum(image=image,selem = selem)
|
||||
|
||||
@exec_and_timeit
|
||||
def cm_dil(image,selem):
|
||||
return dilation(image=image,selem = selem)
|
||||
|
||||
@exec_and_timeit
|
||||
def ctmf_med(image,radius):
|
||||
return median_filter(image=image,radius=radius)
|
||||
|
||||
@exec_and_timeit
|
||||
def ndi_med(image,n):
|
||||
return percentile_filter(image,50,size=n*2-1)
|
||||
|
||||
"""
|
||||
Comparison between
|
||||
|
||||
* ``rank.maximum``
|
||||
* ``cmorph.dilate``
|
||||
|
||||
on increasing structuring element size
|
||||
"""
|
||||
|
||||
a = data.camera()
|
||||
|
||||
rec = []
|
||||
e_range = range(1,20,2)
|
||||
for r in e_range:
|
||||
elem = disk(r+1)
|
||||
rc,ms_rc = cr_max(a,elem)
|
||||
rcm,ms_rcm = cm_dil(a,elem)
|
||||
rec.append((ms_rc,ms_rcm))
|
||||
|
||||
rec = np.asarray(rec)
|
||||
|
||||
plt.figure()
|
||||
plt.title('increasing element size')
|
||||
plt.ylabel('time (ms)')
|
||||
plt.xlabel('element radius')
|
||||
plt.plot(e_range,rec)
|
||||
plt.legend(['crank.maximum','cmorph.dilate'])
|
||||
|
||||
"""
|
||||
and increasing image size
|
||||
|
||||
.. image:: PLOT2RST.current_figure
|
||||
"""
|
||||
|
||||
r = 9
|
||||
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')
|
||||
(rc,ms_rc) = cr_max(a,elem)
|
||||
(rcm,ms_rcm) = cm_dil(a,elem)
|
||||
rec.append((ms_rc,ms_rcm))
|
||||
|
||||
rec = np.asarray(rec)
|
||||
|
||||
plt.figure()
|
||||
plt.title('increasing image size')
|
||||
plt.ylabel('time (ms)')
|
||||
plt.xlabel('image size')
|
||||
plt.plot(s_range,rec)
|
||||
plt.legend(['crank.maximum','cmorph.dilate'])
|
||||
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
Comparison between:
|
||||
|
||||
* ``rank.median``
|
||||
* ``ctmf.median_filter``
|
||||
* ``ndimage.percentile``
|
||||
|
||||
on increasing structuring element size
|
||||
"""
|
||||
|
||||
|
||||
a = data.camera()
|
||||
|
||||
rec = []
|
||||
e_range = range(2,30,4)
|
||||
for r in e_range:
|
||||
elem = disk(r+1)
|
||||
rc,ms_rc = cr_med(a,elem)
|
||||
rctmf,ms_rctmf = ctmf_med(a,r)
|
||||
rndi,ms_ndi = ndi_med(a,r)
|
||||
rec.append((ms_rc,ms_rctmf,ms_ndi))
|
||||
|
||||
rec = np.asarray(rec)
|
||||
|
||||
plt.figure()
|
||||
plt.title('increasing element size')
|
||||
plt.plot(e_range,rec)
|
||||
plt.legend(['rank.median','ctmf.median_filter','ndimage.percentile'])
|
||||
plt.ylabel('time (ms)')
|
||||
plt.xlabel('element radius')
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
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')
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
and increasing image size
|
||||
|
||||
"""
|
||||
|
||||
r = 9
|
||||
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')
|
||||
(rc,ms_rc) = cr_med(a,elem)
|
||||
rctmf,ms_rctmf = ctmf_med(a,r)
|
||||
rndi,ms_ndi = ndi_med(a,r)
|
||||
rec.append((ms_rc,ms_rctmf,ms_ndi))
|
||||
|
||||
rec = np.asarray(rec)
|
||||
|
||||
plt.figure()
|
||||
plt.title('increasing image size')
|
||||
plt.plot(s_range,rec)
|
||||
plt.legend(['rank.median','ctmf.median_filter','ndimage.percentile'])
|
||||
plt.ylabel('time (ms)')
|
||||
plt.xlabel('image size')
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
"""
|
||||
|
||||
plt.show()
|
||||
@@ -0,0 +1,47 @@
|
||||
"""
|
||||
==============================
|
||||
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()
|
||||
@@ -0,0 +1,44 @@
|
||||
"""
|
||||
===================
|
||||
Entropy
|
||||
===================
|
||||
|
||||
|
||||
"""
|
||||
from skimage import data
|
||||
from skimage.filter.rank import entropy
|
||||
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
|
||||
|
||||
ent8 = entropy(a8,disk(5)) # pixel value contain 10x the local entropy
|
||||
ent16 = entropy(a16,disk(5)) # pixel value contain 1000x the local entropy
|
||||
|
||||
# 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.colorbar()
|
||||
|
||||
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.colorbar()
|
||||
plt.show()
|
||||
|
||||
@@ -0,0 +1,84 @@
|
||||
"""
|
||||
===============================
|
||||
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.
|
||||
|
||||
.. [1] http://en.wikipedia.org/wiki/Histogram_equalization
|
||||
.. [2] http://en.wikipedia.org/wiki/Adaptive_histogram_equalization
|
||||
|
||||
"""
|
||||
|
||||
from skimage import data
|
||||
from skimage.util.dtype import dtype_range
|
||||
from skimage import exposure
|
||||
from skimage.morphology import disk
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
import numpy as np
|
||||
from skimage.filter import rank
|
||||
|
||||
def plot_img_and_hist(img, axes, bins=256):
|
||||
"""Plot an image along with its histogram and cumulative histogram.
|
||||
|
||||
"""
|
||||
ax_img, ax_hist = axes
|
||||
ax_cdf = ax_hist.twinx()
|
||||
|
||||
# Display image
|
||||
ax_img.imshow(img, cmap=plt.cm.gray)
|
||||
ax_img.set_axis_off()
|
||||
|
||||
# Display histogram
|
||||
ax_hist.hist(img.ravel(), bins=bins)
|
||||
ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
|
||||
ax_hist.set_xlabel('Pixel intensity')
|
||||
|
||||
xmin, xmax = dtype_range[img.dtype.type]
|
||||
ax_hist.set_xlim(xmin, xmax)
|
||||
|
||||
# Display cumulative distribution
|
||||
img_cdf, bins = exposure.cumulative_distribution(img, bins)
|
||||
ax_cdf.plot(bins, img_cdf, 'r')
|
||||
|
||||
return ax_img, ax_hist, ax_cdf
|
||||
|
||||
|
||||
# Load an example image
|
||||
img = data.moon()
|
||||
|
||||
# Contrast stretching
|
||||
p2 = np.percentile(img, 2)
|
||||
p98 = np.percentile(img, 98)
|
||||
img_rescale = exposure.equalize(img)
|
||||
|
||||
# Equalization
|
||||
selem = disk(30)
|
||||
img_eq = rank.equalize(img,selem=selem)
|
||||
|
||||
|
||||
# Display results
|
||||
f, axes = plt.subplots(2, 3, figsize=(8, 4))
|
||||
|
||||
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0])
|
||||
ax_img.set_title('Low contrast image')
|
||||
ax_hist.set_ylabel('Number of pixels')
|
||||
|
||||
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_rescale, axes[:, 1])
|
||||
ax_img.set_title('Global equalise')
|
||||
|
||||
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_eq, axes[:, 2])
|
||||
ax_img.set_title('Local equalize')
|
||||
ax_cdf.set_ylabel('Fraction of total intensity')
|
||||
|
||||
|
||||
# prevent overlap of y-axis labels
|
||||
plt.subplots_adjust(wspace=0.4)
|
||||
plt.show()
|
||||
|
||||
@@ -0,0 +1,54 @@
|
||||
"""
|
||||
================================
|
||||
Markers for watershed transform
|
||||
================================
|
||||
|
||||
The watershed is a classical algorithm used for **segmentation**, that
|
||||
is, for separating different objects in an image.
|
||||
|
||||
Here a marker image is build from the region of low gradient inside the image.
|
||||
|
||||
See Wikipedia_ for more details on the algorithm.
|
||||
|
||||
.. _Wikipedia: http://en.wikipedia.org/wiki/Watershed_(image_processing)
|
||||
|
||||
"""
|
||||
|
||||
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
|
||||
|
||||
image = data.camera()
|
||||
|
||||
# denoise image
|
||||
denoised = rank.median(image,disk(2))
|
||||
|
||||
# find continuous region (low gradient) --> markers
|
||||
markers = rank.gradient(denoised,disk(5))<10
|
||||
markers = ndimage.label(markers)[0]
|
||||
|
||||
#local gradient
|
||||
gradient = rank.gradient(denoised,disk(2))
|
||||
|
||||
# process the watershed
|
||||
labels = watershed(gradient, markers)
|
||||
|
||||
# display results
|
||||
fig, axes = plt.subplots(ncols=4, figsize=(8, 2.7))
|
||||
ax0, ax1, ax2, ax3 = axes
|
||||
|
||||
ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
|
||||
ax1.imshow(gradient, cmap=plt.cm.spectral, interpolation='nearest')
|
||||
ax2.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest')
|
||||
ax3.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
|
||||
ax3.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest',alpha=.7)
|
||||
|
||||
for ax in axes:
|
||||
ax.axis('off')
|
||||
|
||||
plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
|
||||
plt.show()
|
||||
@@ -26,7 +26,7 @@ See Wikipedia_ for more details on the algorithm.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
from scipy import ndimage
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
from skimage.morphology import watershed, is_local_maximum
|
||||
|
||||
|
||||
@@ -0,0 +1,35 @@
|
||||
To use this to build your Cython file use the commandline options:
|
||||
|
||||
.. sourcecode:: text
|
||||
|
||||
$ python setup.py build_ext --inplace
|
||||
|
||||
|
||||
**To do**
|
||||
|
||||
* add simple examples, adapt documentation on existing examples
|
||||
|
||||
* add/check existing doc
|
||||
|
||||
* adapting tests for each type of filter
|
||||
|
||||
**General remarks**
|
||||
|
||||
Basically these filters compute local histogram for each pixel. Histogram is build using a moving window in
|
||||
order to limit redundant computation. The path followed by the moving window is given hereunder
|
||||
|
||||
...-----------------------\
|
||||
/--------------------------/
|
||||
\-------------------------- ...
|
||||
|
||||
A comparison is proposed with cmorph.dilate algorithm to show how computation costs evolve with respect to image size or
|
||||
structuring element size. This implementation gives better results for large structuring elements.
|
||||
|
||||
A local histogram is update at each pixel by introducing pixel entering the structuring element border and
|
||||
by removing those leaving it. The histogram size is 8bit (256 bins) for 8 bit images and 2 to 12 bit (up to 4096 bins)
|
||||
for 16bit image depending on the image maximum value. Image with pixels higher than 4095 raise a ValueError.
|
||||
|
||||
The filter is applied up to the image border, the neighboorhood used is adjusted accordingly. The user may provide
|
||||
a mask image (same size as input image) where non zero value are the part of the image participating the the
|
||||
histogram computation. By default all the image is filtered.
|
||||
|
||||
@@ -0,0 +1,3 @@
|
||||
from .rank import *
|
||||
from .percentile_rank import *
|
||||
from .bilateral_rank import *
|
||||
@@ -0,0 +1,18 @@
|
||||
cimport numpy as np
|
||||
|
||||
#---------------------------------------------------------------------------
|
||||
# 16 bit core kernel receives extra information about data bitdepth
|
||||
#---------------------------------------------------------------------------
|
||||
|
||||
# generic cdef functions
|
||||
cdef int int_max(int a, int b)
|
||||
cdef int int_min(int a, int b)
|
||||
|
||||
cdef _core16(
|
||||
np.uint16_t kernel(Py_ssize_t * , float, np.uint16_t, Py_ssize_t, Py_ssize_t, Py_ssize_t, float, float, Py_ssize_t, Py_ssize_t),
|
||||
np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask,
|
||||
np.ndarray[np.uint16_t, ndim=2] out,
|
||||
char shift_x, char shift_y, Py_ssize_t bitdepth,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1)
|
||||
@@ -0,0 +1,283 @@
|
||||
#cython: cdivision=True
|
||||
#cython: boundscheck=False
|
||||
#cython: nonecheck=False
|
||||
#cython: wraparound=False
|
||||
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
from libc.stdlib cimport malloc, free
|
||||
|
||||
#---------------------------------------------------------------------------
|
||||
# 16 bit core kernel receives extra information about data bitdepth
|
||||
#---------------------------------------------------------------------------
|
||||
|
||||
# generic cdef functions
|
||||
cdef inline int int_max(int a, int b):
|
||||
return a if a >= b else b
|
||||
cdef inline int int_min(int a, int b):
|
||||
return a if a <= b else b
|
||||
|
||||
cdef inline void histogram_increment(Py_ssize_t * histo, float * pop, np.uint16_t value):
|
||||
histo[value] += 1
|
||||
pop[0] += 1.
|
||||
|
||||
cdef inline void histogram_decrement(Py_ssize_t * histo, float * pop, np.uint16_t value):
|
||||
histo[value] -= 1
|
||||
pop[0] -= 1.
|
||||
|
||||
cdef inline np.uint8_t is_in_mask(Py_ssize_t rows, Py_ssize_t cols, Py_ssize_t r, Py_ssize_t c, np.uint8_t * mask):
|
||||
""" returns 1 if given(r,c) coordinate are within the image frame ([0-rows],[0-cols]) and
|
||||
inside the given mask
|
||||
returns 0 otherwise
|
||||
"""
|
||||
if r < 0 or r > rows - 1 or c < 0 or c > cols - 1:
|
||||
return 0
|
||||
else:
|
||||
if mask[r * cols + c]:
|
||||
return 1
|
||||
else:
|
||||
return 0
|
||||
|
||||
|
||||
cdef inline _core16(
|
||||
np.uint16_t kernel(Py_ssize_t *, float, np.uint16_t, Py_ssize_t, Py_ssize_t, Py_ssize_t, float, float, Py_ssize_t, Py_ssize_t),
|
||||
np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask,
|
||||
np.ndarray[np.uint16_t, ndim=2] out,
|
||||
char shift_x, char shift_y, Py_ssize_t bitdepth,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
""" Main loop, this function computes the histogram for each image point
|
||||
- data is uint8
|
||||
- result is uint8 casted
|
||||
"""
|
||||
|
||||
cdef Py_ssize_t rows = image.shape[0]
|
||||
cdef Py_ssize_t cols = image.shape[1]
|
||||
cdef Py_ssize_t srows = selem.shape[0]
|
||||
cdef Py_ssize_t scols = selem.shape[1]
|
||||
|
||||
cdef Py_ssize_t centre_r = int(selem.shape[0] / 2) + shift_y
|
||||
cdef Py_ssize_t centre_c = int(selem.shape[1] / 2) + shift_x
|
||||
|
||||
# check that structuring element center is inside the element bounding box
|
||||
assert centre_r >= 0
|
||||
assert centre_c >= 0
|
||||
assert centre_r < srows
|
||||
assert centre_c < scols
|
||||
assert bitdepth in range(2, 13)
|
||||
|
||||
maxbin_list = [0, 0, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
|
||||
midbin_list = [0, 0, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048]
|
||||
|
||||
#set maxbin and midbin
|
||||
cdef Py_ssize_t maxbin = maxbin_list[bitdepth], midbin = midbin_list[bitdepth]
|
||||
|
||||
assert (image < maxbin).all()
|
||||
|
||||
image = np.ascontiguousarray(image)
|
||||
|
||||
if mask is None:
|
||||
mask = np.ones((rows, cols), dtype=np.uint8)
|
||||
else:
|
||||
mask = np.ascontiguousarray(mask)
|
||||
|
||||
if image is out:
|
||||
raise NotImplementedError("Cannot perform rank operation in place.")
|
||||
|
||||
if out is None:
|
||||
out = np.zeros((rows, cols), dtype=np.uint16)
|
||||
else:
|
||||
out = np.ascontiguousarray(out)
|
||||
|
||||
mask = np.ascontiguousarray(mask)
|
||||
|
||||
# define pointers to the data
|
||||
cdef np.uint16_t * out_data = <np.uint16_t * >out.data
|
||||
cdef np.uint16_t * image_data = <np.uint16_t * >image.data
|
||||
cdef np.uint8_t * mask_data = <np.uint8_t * >mask.data
|
||||
|
||||
# define local variable types
|
||||
cdef Py_ssize_t r, c, rr, cc, s, value, local_max, i, even_row
|
||||
cdef float pop # number of pixels actually inside the neighborhood (float)
|
||||
|
||||
# allocate memory with malloc
|
||||
cdef Py_ssize_t max_se = srows * scols
|
||||
|
||||
# number of element in each attack border
|
||||
cdef Py_ssize_t num_se_n, num_se_s, num_se_e, num_se_w
|
||||
|
||||
# the current local histogram distribution
|
||||
cdef Py_ssize_t * histo = <Py_ssize_t * >malloc(maxbin * sizeof(Py_ssize_t))
|
||||
|
||||
# these lists contain the relative pixel row and column for each of the 4 attack borders
|
||||
# east, west, north and south
|
||||
# e.g. se_e_r lists the rows of the east structuring element border
|
||||
|
||||
cdef Py_ssize_t * se_e_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_e_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_w_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_w_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_n_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_n_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_s_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_s_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
|
||||
# build attack and release borders
|
||||
# by using difference along axis
|
||||
|
||||
t = np.hstack((selem, np.zeros((selem.shape[0], 1))))
|
||||
t_e = np.diff(t, axis=1) == -1
|
||||
|
||||
t = np.hstack((np.zeros((selem.shape[0], 1)), selem))
|
||||
t_w = np.diff(t, axis=1) == 1
|
||||
|
||||
t = np.vstack((selem, np.zeros((1, selem.shape[1]))))
|
||||
t_s = np.diff(t, axis=0) == -1
|
||||
|
||||
t = np.vstack((np.zeros((1, selem.shape[1])), selem))
|
||||
t_n = np.diff(t, axis=0) == 1
|
||||
|
||||
num_se_n = num_se_s = num_se_e = num_se_w = 0
|
||||
|
||||
for r in range(srows):
|
||||
for c in range(scols):
|
||||
if t_e[r, c]:
|
||||
se_e_r[num_se_e] = r - centre_r
|
||||
se_e_c[num_se_e] = c - centre_c
|
||||
num_se_e += 1
|
||||
if t_w[r, c]:
|
||||
se_w_r[num_se_w] = r - centre_r
|
||||
se_w_c[num_se_w] = c - centre_c
|
||||
num_se_w += 1
|
||||
if t_n[r, c]:
|
||||
se_n_r[num_se_n] = r - centre_r
|
||||
se_n_c[num_se_n] = c - centre_c
|
||||
num_se_n += 1
|
||||
if t_s[r, c]:
|
||||
se_s_r[num_se_s] = r - centre_r
|
||||
se_s_c[num_se_s] = c - centre_c
|
||||
num_se_s += 1
|
||||
|
||||
# initial population and histogram
|
||||
for i in range(maxbin):
|
||||
histo[i] = 0
|
||||
|
||||
pop = 0
|
||||
|
||||
for r in range(srows):
|
||||
for c in range(scols):
|
||||
rr = r - centre_r
|
||||
cc = c - centre_c
|
||||
if selem[r, c]:
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_increment(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
r = 0
|
||||
c = 0
|
||||
# kernel -------------------------------------------
|
||||
out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c],
|
||||
bitdepth, maxbin, midbin, p0, p1, s0, s1)
|
||||
# kernel -------------------------------------------
|
||||
|
||||
# main loop
|
||||
r = 0
|
||||
for even_row in range(0, rows, 2):
|
||||
# ---> west to east
|
||||
for c in range(1, cols):
|
||||
for s in range(num_se_e):
|
||||
rr = r + se_e_r[s]
|
||||
cc = c + se_e_c[s]
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_increment(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
for s in range(num_se_w):
|
||||
rr = r + se_w_r[s]
|
||||
cc = c + se_w_c[s] - 1
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_decrement(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
# kernel -------------------------------------------
|
||||
out_data[r * cols + c] = kernel(
|
||||
histo, pop, image_data[r * cols + c],
|
||||
bitdepth, maxbin, midbin, p0, p1, s0, s1)
|
||||
# kernel -------------------------------------------
|
||||
|
||||
r += 1 # pass to the next row
|
||||
if r >= rows:
|
||||
break
|
||||
|
||||
# ---> north to south
|
||||
for s in range(num_se_s):
|
||||
rr = r + se_s_r[s]
|
||||
cc = c + se_s_c[s]
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_increment(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
for s in range(num_se_n):
|
||||
rr = r + se_n_r[s] - 1
|
||||
cc = c + se_n_c[s]
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_decrement(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
# kernel -------------------------------------------
|
||||
out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c],
|
||||
bitdepth, maxbin, midbin, p0, p1, s0, s1)
|
||||
# kernel -------------------------------------------
|
||||
|
||||
# ---> east to west
|
||||
for c in range(cols - 2, -1, -1):
|
||||
for s in range(num_se_w):
|
||||
rr = r + se_w_r[s]
|
||||
cc = c + se_w_c[s]
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_increment(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
for s in range(num_se_e):
|
||||
rr = r + se_e_r[s]
|
||||
cc = c + se_e_c[s] + 1
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_decrement(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
# kernel -------------------------------------------
|
||||
out_data[r * cols + c] = kernel(
|
||||
histo, pop, image_data[r * cols + c],
|
||||
bitdepth, maxbin, midbin, p0, p1, s0, s1)
|
||||
# kernel -------------------------------------------
|
||||
|
||||
r += 1 # pass to the next row
|
||||
if r >= rows:
|
||||
break
|
||||
|
||||
# ---> north to south
|
||||
for s in range(num_se_s):
|
||||
rr = r + se_s_r[s]
|
||||
cc = c + se_s_c[s]
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_increment(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
for s in range(num_se_n):
|
||||
rr = r + se_n_r[s] - 1
|
||||
cc = c + se_n_c[s]
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_decrement(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
# kernel -------------------------------------------
|
||||
out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c],
|
||||
bitdepth, maxbin, midbin, p0, p1, s0, s1)
|
||||
# kernel -------------------------------------------
|
||||
|
||||
# release memory allocated by malloc
|
||||
|
||||
free(se_e_r)
|
||||
free(se_e_c)
|
||||
free(se_w_r)
|
||||
free(se_w_c)
|
||||
free(se_n_r)
|
||||
free(se_n_c)
|
||||
free(se_s_r)
|
||||
free(se_s_c)
|
||||
|
||||
free(histo)
|
||||
|
||||
return out
|
||||
@@ -0,0 +1,17 @@
|
||||
cimport numpy as np
|
||||
|
||||
# generic cdef functions
|
||||
cdef np.uint8_t uint8_max(np.uint8_t a, np.uint8_t b)
|
||||
cdef np.uint8_t uint8_min(np.uint8_t a, np.uint8_t b)
|
||||
|
||||
#---------------------------------------------------------------------------
|
||||
# 8 bit core kernel receives extra information about data inferior and superior percentiles
|
||||
#---------------------------------------------------------------------------
|
||||
|
||||
cdef _core8(
|
||||
np.uint8_t kernel(Py_ssize_t *, float, np.uint8_t, float, float, Py_ssize_t, Py_ssize_t),
|
||||
np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask,
|
||||
np.ndarray[np.uint8_t, ndim=2] out,
|
||||
char shift_x, char shift_y, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1)
|
||||
@@ -0,0 +1,274 @@
|
||||
#cython: cdivision=True
|
||||
#cython: boundscheck=False
|
||||
#cython: nonecheck=False
|
||||
#cython: wraparound=False
|
||||
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
from libc.stdlib cimport malloc, free
|
||||
|
||||
# generic cdef functions
|
||||
cdef inline np.uint8_t uint8_max(np.uint8_t a, np.uint8_t b):
|
||||
return a if a >= b else b
|
||||
cdef inline np.uint8_t uint8_min(np.uint8_t a, np.uint8_t b):
|
||||
return a if a <= b else b
|
||||
|
||||
|
||||
#---------------------------------------------------------------------------
|
||||
# 8 bit core kernel
|
||||
#---------------------------------------------------------------------------
|
||||
|
||||
cdef inline void histogram_increment(Py_ssize_t * histo, float * pop, np.uint8_t value):
|
||||
histo[value] += 1
|
||||
pop[0] += 1.
|
||||
|
||||
cdef inline void histogram_decrement(Py_ssize_t * histo, float * pop, np.uint8_t value):
|
||||
histo[value] -= 1
|
||||
pop[0] -= 1.
|
||||
|
||||
cdef inline np.uint8_t is_in_mask(Py_ssize_t rows, Py_ssize_t cols, Py_ssize_t r, Py_ssize_t c, np.uint8_t * mask):
|
||||
""" returns 1 if given(r,c) coordinate are within the image frame ([0-rows],[0-cols]) and
|
||||
inside the given mask
|
||||
returns 0 otherwise
|
||||
"""
|
||||
if r < 0 or r > rows - 1 or c < 0 or c > cols - 1:
|
||||
return 0
|
||||
else:
|
||||
if mask[r * cols + c]:
|
||||
return 1
|
||||
else:
|
||||
return 0
|
||||
|
||||
cdef inline _core8(
|
||||
np.uint8_t kernel(Py_ssize_t * , float, np.uint8_t, float, float, Py_ssize_t, Py_ssize_t),
|
||||
np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask,
|
||||
np.ndarray[np.uint8_t, ndim=2] out,
|
||||
char shift_x, char shift_y, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
""" Main loop, this function computes the histogram for each image point
|
||||
- data is uint8
|
||||
- result is uint8 casted
|
||||
"""
|
||||
|
||||
cdef Py_ssize_t rows = image.shape[0]
|
||||
cdef Py_ssize_t cols = image.shape[1]
|
||||
cdef Py_ssize_t srows = selem.shape[0]
|
||||
cdef Py_ssize_t scols = selem.shape[1]
|
||||
|
||||
cdef Py_ssize_t centre_r = int(selem.shape[0] / 2) + shift_y
|
||||
cdef Py_ssize_t centre_c = int(selem.shape[1] / 2) + shift_x
|
||||
|
||||
# check that structuring element center is inside the element bounding box
|
||||
assert centre_r >= 0
|
||||
assert centre_c >= 0
|
||||
assert centre_r < srows
|
||||
assert centre_c < scols
|
||||
|
||||
image = np.ascontiguousarray(image)
|
||||
|
||||
if mask is None:
|
||||
mask = np.ones((rows, cols), dtype=np.uint8)
|
||||
else:
|
||||
mask = np.ascontiguousarray(mask)
|
||||
|
||||
if image is out:
|
||||
raise NotImplementedError("Cannot perform rank operation in place.")
|
||||
|
||||
if out is None:
|
||||
out = np.zeros((rows, cols), dtype=np.uint8)
|
||||
else:
|
||||
out = np.ascontiguousarray(out)
|
||||
|
||||
mask = np.ascontiguousarray(mask)
|
||||
|
||||
# define pointers to the data
|
||||
|
||||
cdef np.uint8_t * out_data = <np.uint8_t * >out.data
|
||||
cdef np.uint8_t * image_data = <np.uint8_t * >image.data
|
||||
cdef np.uint8_t * mask_data = <np.uint8_t * >mask.data
|
||||
|
||||
# define local variable types
|
||||
cdef Py_ssize_t r, c, rr, cc, s, value, local_max, i, even_row
|
||||
|
||||
# number of pixels actually inside the neighborhood (float)
|
||||
cdef float pop
|
||||
|
||||
# allocate memory with malloc
|
||||
cdef Py_ssize_t max_se = srows * scols
|
||||
|
||||
# number of element in each attack border
|
||||
cdef Py_ssize_t num_se_n, num_se_s, num_se_e, num_se_w
|
||||
|
||||
# the current local histogram distribution
|
||||
cdef Py_ssize_t * histo = <Py_ssize_t * >malloc(256 * sizeof(Py_ssize_t))
|
||||
|
||||
# these lists contain the relative pixel row and column for each of the 4 attack borders
|
||||
# east, west, north and south
|
||||
# e.g. se_e_r lists the rows of the east structuring element border
|
||||
|
||||
cdef Py_ssize_t * se_e_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_e_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_w_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_w_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_n_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_n_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_s_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
cdef Py_ssize_t * se_s_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
|
||||
|
||||
# build attack and release borders
|
||||
# by using difference along axis
|
||||
|
||||
t = np.hstack((selem, np.zeros((selem.shape[0], 1))))
|
||||
t_e = np.diff(t, axis=1) == -1
|
||||
|
||||
t = np.hstack((np.zeros((selem.shape[0], 1)), selem))
|
||||
t_w = np.diff(t, axis=1) == 1
|
||||
|
||||
t = np.vstack((selem, np.zeros((1, selem.shape[1]))))
|
||||
t_s = np.diff(t, axis=0) == -1
|
||||
|
||||
t = np.vstack((np.zeros((1, selem.shape[1])), selem))
|
||||
t_n = np.diff(t, axis=0) == 1
|
||||
|
||||
num_se_n = num_se_s = num_se_e = num_se_w = 0
|
||||
|
||||
for r in range(srows):
|
||||
for c in range(scols):
|
||||
if t_e[r, c]:
|
||||
se_e_r[num_se_e] = r - centre_r
|
||||
se_e_c[num_se_e] = c - centre_c
|
||||
num_se_e += 1
|
||||
if t_w[r, c]:
|
||||
se_w_r[num_se_w] = r - centre_r
|
||||
se_w_c[num_se_w] = c - centre_c
|
||||
num_se_w += 1
|
||||
if t_n[r, c]:
|
||||
se_n_r[num_se_n] = r - centre_r
|
||||
se_n_c[num_se_n] = c - centre_c
|
||||
num_se_n += 1
|
||||
if t_s[r, c]:
|
||||
se_s_r[num_se_s] = r - centre_r
|
||||
se_s_c[num_se_s] = c - centre_c
|
||||
num_se_s += 1
|
||||
|
||||
# initial population and histogram (kernel is centered on the first row and column)
|
||||
for i in range(256):
|
||||
histo[i] = 0
|
||||
|
||||
pop = 0
|
||||
|
||||
for r in range(srows):
|
||||
for c in range(scols):
|
||||
rr = r - centre_r
|
||||
cc = c - centre_c
|
||||
if selem[r, c]:
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_increment(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
r = 0
|
||||
c = 0
|
||||
# kernel --------------------------------------------------------------------
|
||||
out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols +
|
||||
c], p0, p1, s0, s1)
|
||||
# kernel --------------------------------------------------------------------
|
||||
|
||||
# main loop
|
||||
r = 0
|
||||
for even_row in range(0, rows, 2):
|
||||
# ---> west to east
|
||||
for c in range(1, cols):
|
||||
for s in range(num_se_e):
|
||||
rr = r + se_e_r[s]
|
||||
cc = c + se_e_c[s]
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_increment(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
for s in range(num_se_w):
|
||||
rr = r + se_w_r[s]
|
||||
cc = c + se_w_c[s] - 1
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_decrement(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
# kernel --------------------------------------------------------------------
|
||||
out_data[r * cols + c] = kernel(
|
||||
histo, pop, image_data[r * cols + c], p0, p1, s0, s1)
|
||||
# kernel --------------------------------------------------------------------
|
||||
|
||||
r += 1 # pass to the next row
|
||||
if r >= rows:
|
||||
break
|
||||
|
||||
# ---> north to south
|
||||
for s in range(num_se_s):
|
||||
rr = r + se_s_r[s]
|
||||
cc = c + se_s_c[s]
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_increment(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
for s in range(num_se_n):
|
||||
rr = r + se_n_r[s] - 1
|
||||
cc = c + se_n_c[s]
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_decrement(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
# kernel --------------------------------------------------------------------
|
||||
out_data[r * cols + c] = kernel(histo, pop, image_data[r *
|
||||
cols + c], p0, p1, s0, s1)
|
||||
# kernel --------------------------------------------------------------------
|
||||
|
||||
# ---> east to west
|
||||
for c in range(cols - 2, -1, -1):
|
||||
for s in range(num_se_w):
|
||||
rr = r + se_w_r[s]
|
||||
cc = c + se_w_c[s]
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_increment(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
for s in range(num_se_e):
|
||||
rr = r + se_e_r[s]
|
||||
cc = c + se_e_c[s] + 1
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_decrement(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
# kernel --------------------------------------------------------------------
|
||||
out_data[r * cols + c] = kernel(
|
||||
histo, pop, image_data[r * cols + c], p0, p1, s0, s1)
|
||||
# kernel --------------------------------------------------------------------
|
||||
|
||||
r += 1 # pass to the next row
|
||||
if r >= rows:
|
||||
break
|
||||
|
||||
# ---> north to south
|
||||
for s in range(num_se_s):
|
||||
rr = r + se_s_r[s]
|
||||
cc = c + se_s_c[s]
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_increment(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
for s in range(num_se_n):
|
||||
rr = r + se_n_r[s] - 1
|
||||
cc = c + se_n_c[s]
|
||||
if is_in_mask(rows, cols, rr, cc, mask_data):
|
||||
histogram_decrement(histo, & pop, image_data[rr * cols + cc])
|
||||
|
||||
# kernel --------------------------------------------------------------------
|
||||
out_data[r * cols + c] = kernel(histo, pop, image_data[r *
|
||||
cols + c], p0, p1, s0, s1)
|
||||
# kernel --------------------------------------------------------------------
|
||||
|
||||
# release memory allocated by malloc
|
||||
|
||||
free(se_e_r)
|
||||
free(se_e_c)
|
||||
free(se_w_r)
|
||||
free(se_w_c)
|
||||
free(se_n_r)
|
||||
free(se_n_c)
|
||||
free(se_s_r)
|
||||
free(se_s_c)
|
||||
|
||||
free(histo)
|
||||
|
||||
return out
|
||||
@@ -0,0 +1,364 @@
|
||||
#cython: cdivision=True
|
||||
#cython: boundscheck=False
|
||||
#cython: nonecheck=False
|
||||
#cython: wraparound=False
|
||||
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
from libc.math cimport log2
|
||||
|
||||
# import main loop
|
||||
from skimage.filter.rank._core16 cimport _core16
|
||||
|
||||
# -----------------------------------------------------------------
|
||||
# kernels uint16 take extra parameter for defining the bitdepth
|
||||
# -----------------------------------------------------------------
|
||||
|
||||
cdef inline np.uint16_t kernel_autolevel(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t i, imin, imax, delta
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin - 1, -1, -1):
|
||||
if histo[i]:
|
||||
imax = i
|
||||
break
|
||||
for i in range(maxbin):
|
||||
if histo[i]:
|
||||
imin = i
|
||||
break
|
||||
delta = imax - imin
|
||||
if delta > 0:
|
||||
return < np.uint16_t > (1. * (maxbin - 1) * (g - imin) / delta)
|
||||
else:
|
||||
return < np.uint16_t > (imax - imin)
|
||||
|
||||
cdef inline np.uint16_t kernel_bottomhat(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t i
|
||||
|
||||
for i in range(maxbin):
|
||||
if histo[i]:
|
||||
break
|
||||
|
||||
return < np.uint16_t > (g - i)
|
||||
|
||||
|
||||
cdef inline np.uint16_t kernel_equalize(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t i
|
||||
cdef float sum = 0.
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin):
|
||||
sum += histo[i]
|
||||
if i >= g:
|
||||
break
|
||||
|
||||
return < np.uint16_t > (((maxbin - 1) * sum) / pop)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_gradient(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t i, imin, imax
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin - 1, -1, -1):
|
||||
if histo[i]:
|
||||
imax = i
|
||||
break
|
||||
for i in range(maxbin):
|
||||
if histo[i]:
|
||||
imin = i
|
||||
break
|
||||
return < np.uint16_t > (imax - imin)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_maximum(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t i
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin - 1, -1, -1):
|
||||
if histo[i]:
|
||||
return < np.uint16_t > (i)
|
||||
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_mean(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t i
|
||||
cdef float mean = 0.
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin):
|
||||
mean += histo[i] * i
|
||||
return < np.uint16_t > (mean / pop)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_meansubstraction(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t i
|
||||
cdef float mean = 0.
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin):
|
||||
mean += histo[i] * i
|
||||
return < np.uint16_t > ((g - mean / pop) / 2. + (midbin - 1))
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_median(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t i
|
||||
cdef float sum = pop / 2.0
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin):
|
||||
if histo[i]:
|
||||
sum -= histo[i]
|
||||
if sum < 0:
|
||||
return < np.uint16_t > (i)
|
||||
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_minimum(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t i
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin):
|
||||
if histo[i]:
|
||||
return < np.uint16_t > (i)
|
||||
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_modal(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t hmax = 0, imax = 0
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin):
|
||||
if histo[i] > hmax:
|
||||
hmax = histo[i]
|
||||
imax = i
|
||||
return < np.uint16_t > (imax)
|
||||
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_morph_contr_enh(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t i, imin, imax
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin - 1, -1, -1):
|
||||
if histo[i]:
|
||||
imax = i
|
||||
break
|
||||
for i in range(maxbin):
|
||||
if histo[i]:
|
||||
imin = i
|
||||
break
|
||||
if imax - g < g - imin:
|
||||
return < np.uint16_t > (imax)
|
||||
else:
|
||||
return < np.uint16_t > (imin)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_pop(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
return < np.uint16_t > (pop)
|
||||
|
||||
cdef inline np.uint16_t kernel_threshold(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t i
|
||||
cdef float mean = 0.
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin):
|
||||
mean += histo[i] * i
|
||||
return < np.uint16_t > (g > (mean / pop))
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_tophat(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t i
|
||||
|
||||
for i in range(maxbin - 1, -1, -1):
|
||||
if histo[i]:
|
||||
break
|
||||
|
||||
return < np.uint16_t > (i - g)
|
||||
|
||||
|
||||
cdef inline np.uint16_t kernel_entropy(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g,
|
||||
Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
|
||||
float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef Py_ssize_t i
|
||||
cdef float e,p
|
||||
|
||||
e = 0.
|
||||
|
||||
for i in range(maxbin):
|
||||
p = histo[i]/pop
|
||||
if p>0:
|
||||
e -= p*log2(p)
|
||||
|
||||
return < np.uint16_t > e*1000
|
||||
|
||||
# -----------------------------------------------------------------
|
||||
# python wrappers
|
||||
# -----------------------------------------------------------------
|
||||
|
||||
|
||||
def autolevel(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_autolevel, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def bottomhat(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_bottomhat, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def equalize(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_equalize, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def gradient(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_gradient, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def maximum(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_maximum, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def mean(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_mean, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def meansubstraction(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_meansubstraction, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def median(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_median, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def minimum(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_minimum, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def morph_contr_enh(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def modal(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_modal, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def pop(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_pop, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def threshold(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_threshold, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def tophat(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_tophat, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
def entropy(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
|
||||
return _core16(kernel_entropy, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
@@ -0,0 +1,72 @@
|
||||
#cython: cdivision=True
|
||||
#cython: boundscheck=False
|
||||
#cython: nonecheck=False
|
||||
#cython: wraparound=False
|
||||
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
|
||||
# import main loop
|
||||
from skimage.filter.rank._core16 cimport _core16
|
||||
|
||||
# -----------------------------------------------------------------
|
||||
# kernels uint16 take extra parameter for defining the bitdepth
|
||||
# -----------------------------------------------------------------
|
||||
|
||||
|
||||
cdef inline np.uint16_t kernel_mean(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
|
||||
Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef int i, bilat_pop = 0
|
||||
cdef float mean = 0.
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin):
|
||||
if (g > (i - s0)) and (g < (i + s1)):
|
||||
bilat_pop += histo[i]
|
||||
mean += histo[i] * i
|
||||
if bilat_pop:
|
||||
return < np.uint16_t > (mean / bilat_pop)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
|
||||
cdef inline np.uint16_t kernel_pop(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth, Py_ssize_t maxbin,
|
||||
Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef int i, bilat_pop = 0
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin):
|
||||
if (g > (i - s0)) and (g < (i + s1)):
|
||||
bilat_pop += histo[i]
|
||||
return < np.uint16_t > (bilat_pop)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
|
||||
# -----------------------------------------------------------------
|
||||
# python wrappers
|
||||
# -----------------------------------------------------------------
|
||||
def mean(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
|
||||
"""average gray level (clipped on uint8)
|
||||
"""
|
||||
return _core16(kernel_mean, image, selem, mask, out, shift_x, shift_y, bitdepth, 0., 0., s0, s1)
|
||||
|
||||
|
||||
def pop(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
|
||||
"""returns the number of actual pixels of the structuring element inside the mask
|
||||
"""
|
||||
return _core16(kernel_pop, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, s0, s1)
|
||||
@@ -0,0 +1,297 @@
|
||||
#cython: cdivision=True
|
||||
#cython: boundscheck=False
|
||||
#cython: nonecheck=False
|
||||
#cython: wraparound=False
|
||||
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
|
||||
# import main loop
|
||||
from skimage.filter.rank._core16 cimport _core16, int_min, int_max
|
||||
|
||||
# -----------------------------------------------------------------
|
||||
# kernels uint16 (SOFT version using percentiles)
|
||||
# -----------------------------------------------------------------
|
||||
|
||||
cdef inline np.uint16_t kernel_autolevel(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
|
||||
Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef int i, imin, imax, sum, delta
|
||||
|
||||
if pop:
|
||||
sum = 0
|
||||
p1 = 1.0 - p1
|
||||
for i in range(maxbin):
|
||||
sum += histo[i]
|
||||
if sum > p0 * pop:
|
||||
imin = i
|
||||
break
|
||||
sum = 0
|
||||
for i in range(maxbin - 1, -1, -1):
|
||||
sum += histo[i]
|
||||
if sum > p1 * pop:
|
||||
imax = i
|
||||
break
|
||||
|
||||
delta = imax - imin
|
||||
if delta > 0:
|
||||
return < np.uint16_t > (1.0 * (maxbin - 1) * (int_min(int_max(imin, g), imax) - imin) / delta)
|
||||
else:
|
||||
return < np.uint16_t > (imax - imin)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
|
||||
cdef inline np.uint16_t kernel_gradient(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
|
||||
Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef int i, imin, imax, sum, delta
|
||||
|
||||
if pop:
|
||||
sum = 0
|
||||
p1 = 1.0 - p1
|
||||
for i in range(maxbin):
|
||||
sum += histo[i]
|
||||
if sum >= p0 * pop:
|
||||
imin = i
|
||||
break
|
||||
sum = 0
|
||||
for i in range((maxbin - 1), -1, -1):
|
||||
sum += histo[i]
|
||||
if sum >= p1 * pop:
|
||||
imax = i
|
||||
break
|
||||
|
||||
return < np.uint16_t > (imax - imin)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
|
||||
cdef inline np.uint16_t kernel_mean(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
|
||||
Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef int i, sum, mean, n
|
||||
|
||||
if pop:
|
||||
sum = 0
|
||||
mean = 0
|
||||
n = 0
|
||||
for i in range(maxbin):
|
||||
sum += histo[i]
|
||||
if (sum >= p0 * pop) and (sum <= p1 * pop):
|
||||
n += histo[i]
|
||||
mean += histo[i] * i
|
||||
|
||||
if n > 0:
|
||||
return < np.uint16_t > (1.0 * mean / n)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_mean_substraction(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
|
||||
Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef int i, sum, mean, n
|
||||
|
||||
if pop:
|
||||
sum = 0
|
||||
mean = 0
|
||||
n = 0
|
||||
for i in range(maxbin):
|
||||
sum += histo[i]
|
||||
if (sum >= p0 * pop) and (sum <= p1 * pop):
|
||||
n += histo[i]
|
||||
mean += histo[i] * i
|
||||
if n > 0:
|
||||
return < np.uint16_t > ((g - (mean / n)) * .5 + midbin)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_morph_contr_enh(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
|
||||
Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef int i, imin, imax, sum, delta
|
||||
|
||||
if pop:
|
||||
sum = 0
|
||||
p1 = 1.0 - p1
|
||||
for i in range(maxbin):
|
||||
sum += histo[i]
|
||||
if sum > p0 * pop:
|
||||
imin = i
|
||||
break
|
||||
sum = 0
|
||||
for i in range((maxbin - 1), -1, -1):
|
||||
sum += histo[i]
|
||||
if sum > p1 * pop:
|
||||
imax = i
|
||||
break
|
||||
if g > imax:
|
||||
return < np.uint16_t > imax
|
||||
if g < imin:
|
||||
return < np.uint16_t > imin
|
||||
if imax - g < g - imin:
|
||||
return < np.uint16_t > imax
|
||||
else:
|
||||
return < np.uint16_t > imin
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_percentile(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
|
||||
Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef int i
|
||||
cdef float sum = 0.
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin):
|
||||
sum += histo[i]
|
||||
if sum >= p0 * pop:
|
||||
break
|
||||
|
||||
return < np.uint16_t > (i)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_pop(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
|
||||
Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef int i, sum, n
|
||||
|
||||
if pop:
|
||||
sum = 0
|
||||
n = 0
|
||||
for i in range(maxbin):
|
||||
sum += histo[i]
|
||||
if (sum >= p0 * pop) and (sum <= p1 * pop):
|
||||
n += histo[i]
|
||||
return < np.uint16_t > (n)
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
cdef inline np.uint16_t kernel_threshold(
|
||||
Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
|
||||
Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef int i
|
||||
cdef float sum = 0.
|
||||
|
||||
if pop:
|
||||
for i in range(maxbin):
|
||||
sum += histo[i]
|
||||
if sum >= p0 * pop:
|
||||
break
|
||||
|
||||
return < np.uint16_t > ((maxbin - 1) * (g >= i))
|
||||
else:
|
||||
return < np.uint16_t > (0)
|
||||
|
||||
# -----------------------------------------------------------------
|
||||
# python wrappers
|
||||
# -----------------------------------------------------------------
|
||||
|
||||
|
||||
def autolevel(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
|
||||
"""bottom hat
|
||||
"""
|
||||
return _core16(
|
||||
kernel_autolevel, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1,
|
||||
< Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def gradient(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
|
||||
"""return p0,p1 percentile gradient
|
||||
"""
|
||||
return _core16(
|
||||
kernel_gradient, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def mean(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
|
||||
"""return mean between [p0 and p1] percentiles
|
||||
"""
|
||||
return _core16(
|
||||
kernel_mean, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def mean_substraction(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
|
||||
"""return original - mean between [p0 and p1] percentiles *.5 +127
|
||||
"""
|
||||
return _core16(
|
||||
kernel_mean_substraction, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1,
|
||||
< Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def morph_contr_enh(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
|
||||
"""reforce contrast using percentiles
|
||||
"""
|
||||
return _core16(
|
||||
kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1,
|
||||
< Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def percentile(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
|
||||
"""return p0 percentile
|
||||
"""
|
||||
return _core16(
|
||||
kernel_percentile, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1,
|
||||
< Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def pop(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
|
||||
"""return nb of pixels between [p0 and p1]
|
||||
"""
|
||||
return _core16(
|
||||
kernel_pop, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1,
|
||||
< Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def threshold(np.ndarray[np.uint16_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint16_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
|
||||
"""return (maxbin-1) if g > percentile p0
|
||||
"""
|
||||
return _core16(
|
||||
kernel_threshold, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1,
|
||||
< Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
@@ -0,0 +1,458 @@
|
||||
#cython: cdivision=True
|
||||
#cython: boundscheck=False
|
||||
#cython: nonecheck=False
|
||||
#cython: wraparound=False
|
||||
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
|
||||
from libc.math cimport log2
|
||||
|
||||
# import main loop
|
||||
from skimage.filter.rank._core8 cimport _core8
|
||||
|
||||
# -----------------------------------------------------------------
|
||||
# kernels uint8
|
||||
# -----------------------------------------------------------------
|
||||
|
||||
cdef inline np.uint8_t kernel_autolevel(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1,
|
||||
Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i, imin, imax, delta
|
||||
|
||||
if pop:
|
||||
for i in range(255, -1, -1):
|
||||
if histo[i]:
|
||||
imax = i
|
||||
break
|
||||
for i in range(256):
|
||||
if histo[i]:
|
||||
imin = i
|
||||
break
|
||||
delta = imax - imin
|
||||
if delta > 0:
|
||||
return < np.uint8_t > (255. * (g - imin) / delta)
|
||||
else:
|
||||
return < np.uint8_t > (imax - imin)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_bottomhat(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1,
|
||||
Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i
|
||||
|
||||
for i in range(256):
|
||||
if histo[i]:
|
||||
break
|
||||
|
||||
return < np.uint8_t > (g - i)
|
||||
|
||||
|
||||
cdef inline np.uint8_t kernel_equalize(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1,
|
||||
Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i
|
||||
cdef float sum = 0.
|
||||
|
||||
if pop:
|
||||
for i in range(256):
|
||||
sum += histo[i]
|
||||
if i >= g:
|
||||
break
|
||||
|
||||
return < np.uint8_t > ((255 * sum) / pop)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_gradient(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i, imin, imax
|
||||
|
||||
if pop:
|
||||
for i in range(255, -1, -1):
|
||||
if histo[i]:
|
||||
imax = i
|
||||
break
|
||||
for i in range(256):
|
||||
if histo[i]:
|
||||
imin = i
|
||||
break
|
||||
return < np.uint8_t > (imax - imin)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_maximum(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i
|
||||
|
||||
if pop:
|
||||
for i in range(255, -1, -1):
|
||||
if histo[i]:
|
||||
return < np.uint8_t > (i)
|
||||
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_mean(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i
|
||||
cdef float mean = 0.
|
||||
|
||||
if pop:
|
||||
for i in range(256):
|
||||
mean += histo[i] * i
|
||||
return < np.uint8_t > (mean / pop)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_meansubstraction(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1,
|
||||
Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i
|
||||
cdef float mean = 0.
|
||||
|
||||
if pop:
|
||||
for i in range(256):
|
||||
mean += histo[i] * i
|
||||
return < np.uint8_t > ((g - mean / pop) / 2. + 127)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_median(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i
|
||||
cdef float sum = pop / 2.0
|
||||
|
||||
if pop:
|
||||
for i in range(256):
|
||||
if histo[i]:
|
||||
sum -= histo[i]
|
||||
if sum < 0:
|
||||
return < np.uint8_t > (i)
|
||||
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_minimum(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i
|
||||
|
||||
if pop:
|
||||
for i in range(256):
|
||||
if histo[i]:
|
||||
return < np.uint8_t > (i)
|
||||
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_modal(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t hmax = 0, imax = 0
|
||||
|
||||
if pop:
|
||||
for i in range(256):
|
||||
if histo[i] > hmax:
|
||||
hmax = histo[i]
|
||||
imax = i
|
||||
return < np.uint8_t > (imax)
|
||||
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_morph_contr_enh(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1,
|
||||
Py_ssize_t s0, Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i, imin, imax
|
||||
|
||||
if pop:
|
||||
for i in range(255, -1, -1):
|
||||
if histo[i]:
|
||||
imax = i
|
||||
break
|
||||
for i in range(256):
|
||||
if histo[i]:
|
||||
imin = i
|
||||
break
|
||||
if imax - g < g - imin:
|
||||
return < np.uint8_t > (imax)
|
||||
else:
|
||||
return < np.uint8_t > (imin)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_pop(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
|
||||
return < np.uint8_t > (pop)
|
||||
|
||||
cdef inline np.uint8_t kernel_threshold(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i
|
||||
cdef float mean = 0.
|
||||
|
||||
if pop:
|
||||
for i in range(256):
|
||||
mean += histo[i] * i
|
||||
return < np.uint8_t > (g > (mean / pop))
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_tophat(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i
|
||||
|
||||
for i in range(255, -1, -1):
|
||||
if histo[i]:
|
||||
break
|
||||
|
||||
return < np.uint8_t > (i - g)
|
||||
|
||||
cdef inline np.uint8_t kernel_noise_filter(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i
|
||||
cdef Py_ssize_t min_i
|
||||
|
||||
# early stop if at least one pixel of the neighborhood has the same g
|
||||
if histo[g]>0:
|
||||
return < np.uint8_t > 0
|
||||
|
||||
for i in range(g, -1, -1):
|
||||
if histo[i]:
|
||||
break
|
||||
min_i = g-i
|
||||
for i in range(g, 256):
|
||||
if histo[i]:
|
||||
break
|
||||
if i-g < min_i:
|
||||
return < np.uint8_t > (i-g)
|
||||
else:
|
||||
return < np.uint8_t > min_i
|
||||
|
||||
cdef inline np.uint8_t kernel_entropy(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i
|
||||
cdef float e,p
|
||||
|
||||
e = 0.
|
||||
|
||||
for i in range(256):
|
||||
p = histo[i]/pop
|
||||
if p>0:
|
||||
e -= p*log2(p)
|
||||
|
||||
return < np.uint8_t > e*10
|
||||
|
||||
cdef inline np.uint8_t kernel_otsu(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
|
||||
cdef Py_ssize_t i
|
||||
cdef Py_ssize_t max_i
|
||||
cdef float P, mu1, mu2, q1,new_q1, sigma_b, max_sigma_b
|
||||
cdef float mu = 0.
|
||||
|
||||
# compute local mean
|
||||
|
||||
if pop:
|
||||
for i in range(256):
|
||||
mu += histo[i] * i
|
||||
mu = (mu / pop)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
# maximizing the between class variance
|
||||
max_i = 0
|
||||
q1 = histo[0]/pop
|
||||
m1 = 0.
|
||||
max_sigma_b = 0.
|
||||
|
||||
for i in range(1,256):
|
||||
P = histo[i]/pop
|
||||
new_q1 = q1 + P
|
||||
if new_q1>0:
|
||||
mu1 = (q1*mu1 + i*P)/new_q1
|
||||
mu2 = (mu-new_q1*mu1)/(1.-new_q1)
|
||||
sigma_b = new_q1*(1.-new_q1)*(mu1-mu2)**2
|
||||
if sigma_b>max_sigma_b:
|
||||
max_sigma_b = sigma_b
|
||||
max_i = i
|
||||
q1 = new_q1
|
||||
|
||||
if g>max_i:
|
||||
return < np.uint8_t > 255
|
||||
else:
|
||||
return < np.uint8_t > 0
|
||||
# -----------------------------------------------------------------
|
||||
# python wrappers
|
||||
# used only internally
|
||||
# -----------------------------------------------------------------
|
||||
|
||||
|
||||
def autolevel(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(
|
||||
kernel_autolevel, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def bottomhat(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(
|
||||
kernel_bottomhat, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def equalize(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(
|
||||
kernel_equalize, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def gradient(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(
|
||||
kernel_gradient, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def maximum(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(kernel_maximum, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def mean(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(kernel_mean, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def meansubstraction(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(
|
||||
kernel_meansubstraction, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def median(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(kernel_median, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def minimum(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(kernel_minimum, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def morph_contr_enh(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(
|
||||
kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def modal(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(kernel_modal, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def pop(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(kernel_pop, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def threshold(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(
|
||||
kernel_threshold, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def tophat(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(kernel_tophat, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def noise_filter(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(kernel_noise_filter, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
def entropy(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(kernel_entropy, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
def otsu(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0):
|
||||
return _core8(kernel_otsu, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
@@ -0,0 +1,286 @@
|
||||
#cython: cdivision=True
|
||||
#cython: boundscheck=False
|
||||
#cython: nonecheck=False
|
||||
#cython: wraparound=False
|
||||
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
|
||||
# import main loop
|
||||
from skimage.filter.rank._core8 cimport _core8, uint8_max, uint8_min
|
||||
|
||||
# -----------------------------------------------------------------
|
||||
# kernels uint8 (SOFT version using percentiles)
|
||||
# -----------------------------------------------------------------
|
||||
|
||||
cdef inline np.uint8_t kernel_autolevel(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
cdef int i, imin, imax, sum, delta
|
||||
|
||||
if pop:
|
||||
sum = 0
|
||||
p1 = 1.0 - p1
|
||||
imin = 0
|
||||
imax = 255
|
||||
|
||||
for i in range(256):
|
||||
sum += histo[i]
|
||||
if sum > (p0 * pop):
|
||||
imin = i
|
||||
break
|
||||
sum = 0
|
||||
for i in range(255, -1, -1):
|
||||
sum += histo[i]
|
||||
if sum > (p1 * pop):
|
||||
imax = i
|
||||
break
|
||||
delta = imax - imin
|
||||
if delta > 0:
|
||||
return < np.uint8_t > (255 * (uint8_min(uint8_max(imin, g), imax) - imin) / delta)
|
||||
else:
|
||||
return < np.uint8_t > (imax - imin)
|
||||
else:
|
||||
return < np.uint8_t > (128)
|
||||
|
||||
|
||||
cdef inline np.uint8_t kernel_gradient(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
cdef int i, imin, imax, sum, delta
|
||||
|
||||
if pop:
|
||||
sum = 0
|
||||
p1 = 1.0 - p1
|
||||
for i in range(256):
|
||||
sum += histo[i]
|
||||
if sum >= p0 * pop:
|
||||
imin = i
|
||||
break
|
||||
sum = 0
|
||||
for i in range(255, -1, -1):
|
||||
sum += histo[i]
|
||||
if sum >= p1 * pop:
|
||||
imax = i
|
||||
break
|
||||
|
||||
return < np.uint8_t > (imax - imin)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
|
||||
cdef inline np.uint8_t kernel_mean(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
cdef int i, sum, mean, n
|
||||
|
||||
if pop:
|
||||
sum = 0
|
||||
mean = 0
|
||||
n = 0
|
||||
for i in range(256):
|
||||
sum += histo[i]
|
||||
if (sum >= p0 * pop) and (sum <= p1 * pop):
|
||||
n += histo[i]
|
||||
mean += histo[i] * i
|
||||
if n > 0:
|
||||
return < np.uint8_t > (1.0 * mean / n)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_mean_substraction(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1,
|
||||
Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef int i, sum, mean, n
|
||||
|
||||
if pop:
|
||||
sum = 0
|
||||
mean = 0
|
||||
n = 0
|
||||
for i in range(256):
|
||||
sum += histo[i]
|
||||
if (sum >= p0 * pop) and (sum <= p1 * pop):
|
||||
n += histo[i]
|
||||
mean += histo[i] * i
|
||||
if n > 0:
|
||||
return < np.uint8_t > ((g - (mean / n)) * .5 + 127)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_morph_contr_enh(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1,
|
||||
Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef int i, imin, imax, sum, delta
|
||||
|
||||
if pop:
|
||||
sum = 0
|
||||
p1 = 1.0 - p1
|
||||
for i in range(256):
|
||||
sum += histo[i]
|
||||
if sum >= p0 * pop:
|
||||
imin = i
|
||||
break
|
||||
sum = 0
|
||||
for i in range(255, -1, -1):
|
||||
sum += histo[i]
|
||||
if sum >= p1 * pop:
|
||||
imax = i
|
||||
break
|
||||
if g > imax:
|
||||
return < np.uint8_t > imax
|
||||
if g < imin:
|
||||
return < np.uint8_t > imin
|
||||
if imax - g < g - imin:
|
||||
return < np.uint8_t > imax
|
||||
else:
|
||||
return < np.uint8_t > imin
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_percentile(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1,
|
||||
Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef int i
|
||||
cdef float sum = 0.
|
||||
|
||||
if pop:
|
||||
for i in range(256):
|
||||
sum += histo[i]
|
||||
if sum >= p0 * pop:
|
||||
break
|
||||
|
||||
return < np.uint8_t > (i)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_pop(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1,
|
||||
Py_ssize_t s0, Py_ssize_t s1):
|
||||
cdef int i, sum, n
|
||||
|
||||
if pop:
|
||||
sum = 0
|
||||
n = 0
|
||||
for i in range(256):
|
||||
sum += histo[i]
|
||||
if (sum >= p0 * pop) and (sum <= p1 * pop):
|
||||
n += histo[i]
|
||||
return < np.uint8_t > (n)
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
cdef inline np.uint8_t kernel_threshold(
|
||||
Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
|
||||
Py_ssize_t s1):
|
||||
cdef int i
|
||||
cdef float sum = 0.
|
||||
|
||||
if pop:
|
||||
for i in range(256):
|
||||
sum += histo[i]
|
||||
if sum >= p0 * pop:
|
||||
break
|
||||
|
||||
return < np.uint8_t > (255 * (g >= i))
|
||||
else:
|
||||
return < np.uint8_t > (0)
|
||||
|
||||
# -----------------------------------------------------------------
|
||||
# python wrappers
|
||||
# -----------------------------------------------------------------
|
||||
|
||||
|
||||
def autolevel(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
|
||||
"""autolevel
|
||||
"""
|
||||
return _core8(
|
||||
kernel_autolevel, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def gradient(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
|
||||
"""return p0,p1 percentile gradient
|
||||
"""
|
||||
return _core8(
|
||||
kernel_gradient, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def mean(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
|
||||
"""return mean between [p0 and p1] percentiles
|
||||
"""
|
||||
return _core8(kernel_mean, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def mean_substraction(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
|
||||
"""return original - mean between [p0 and p1] percentiles *.5 +127
|
||||
"""
|
||||
return _core8(
|
||||
kernel_mean_substraction, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def morph_contr_enh(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
|
||||
"""reforce contrast using percentiles
|
||||
"""
|
||||
return _core8(
|
||||
kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def percentile(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
|
||||
"""return p0 percentile
|
||||
"""
|
||||
return _core8(
|
||||
kernel_percentile, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
|
||||
|
||||
def pop(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
|
||||
"""return nb of pixels between [p0 and p1]
|
||||
"""
|
||||
return _core8(kernel_pop, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0, < Py_ssize_t > 0)
|
||||
|
||||
|
||||
def threshold(np.ndarray[np.uint8_t, ndim=2] image,
|
||||
np.ndarray[np.uint8_t, ndim=2] selem,
|
||||
np.ndarray[np.uint8_t, ndim=2] mask=None,
|
||||
np.ndarray[np.uint8_t, ndim=2] out=None,
|
||||
char shift_x=0, char shift_y=0, float p0=0., float p1=0.):
|
||||
"""return 255 if g > percentile p0
|
||||
"""
|
||||
return _core8(
|
||||
kernel_threshold, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0,
|
||||
< Py_ssize_t > 0)
|
||||
@@ -0,0 +1,165 @@
|
||||
"""bilateral_rank.py - approximate bilateral rankfilter for local (custom kernel) mean
|
||||
|
||||
The local histogram is computed using a sliding window similar to the method described in
|
||||
|
||||
Reference: Huang, T. ,Yang, G. ; Tang, G.. "A fast two-dimensional median filtering algorithm",
|
||||
IEEE Transactions on Acoustics, Speech and Signal Processing, Feb 1979. Volume: 27 , Issue: 1, Page(s): 13 - 18.
|
||||
|
||||
input image can be 8 bit or 16 bit with a value < 4096 (i.e. 12 bit),
|
||||
8 bit images are casted in 16 bit
|
||||
the number of histogram bins is determined from the maximum value present in the image
|
||||
|
||||
The pixel neighborhood is defined by:
|
||||
|
||||
* the given structuring element
|
||||
|
||||
* an interval [g-s0,g+s1] in gray level around g the processed pixel gray level
|
||||
|
||||
The kernel is flat (i.e. each pixel belonging to the neighborhood contributes equally)
|
||||
|
||||
result image is 16 bit with respect to the input image
|
||||
|
||||
"""
|
||||
|
||||
from skimage import img_as_ubyte
|
||||
|
||||
import numpy as np
|
||||
from skimage.filter.rank import _crank16_bilateral
|
||||
|
||||
from skimage.filter.rank.generic import find_bitdepth
|
||||
|
||||
|
||||
__all__ = ['bilateral_mean', 'bilateral_pop']
|
||||
|
||||
|
||||
def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y, s0, s1):
|
||||
selem = img_as_ubyte(selem)
|
||||
if mask is not None:
|
||||
mask = img_as_ubyte(mask)
|
||||
if image.dtype == np.uint8:
|
||||
image = image.astype(np.uint16)
|
||||
elif image.dtype == np.uint16:
|
||||
pass
|
||||
else:
|
||||
raise TypeError("only uint8 and uint16 image supported!")
|
||||
bitdepth = find_bitdepth(image)
|
||||
if bitdepth > 11:
|
||||
raise ValueError("only uint16 <4096 image (12bit) supported!")
|
||||
return func16(
|
||||
image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, bitdepth=bitdepth + 1, out=out,
|
||||
s0=s0, s1=s1)
|
||||
|
||||
|
||||
def bilateral_mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False, s0=10, s1=10):
|
||||
"""Apply a flat kernel bilateral filter.
|
||||
|
||||
This is an edge-preserving and noise reducing denoising filter. It averages
|
||||
pixels based on their spatial closeness and radiometric similarity.
|
||||
|
||||
Spatial closeness is measured by considering only the local pixel neighborhood given by a
|
||||
structuring element (selem).
|
||||
|
||||
Radiometric similarity is defined by the gray level interval [g-s0,g+s1] where g is the current pixel gray level.
|
||||
Only pixels belonging to the structuring element AND having a gray level inside this interval are averaged.
|
||||
Return greyscale local bilateral_mean of an image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, as the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
s0, s1 : int
|
||||
define the [s0,s1] interval to be considered for computing the value.
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint16 array (uint8 image are casted to uint16)
|
||||
The result of the local bilateral mean.
|
||||
|
||||
See also
|
||||
--------
|
||||
skimage.filter.denoise_bilateral() for a gaussian bilateral filter.
|
||||
|
||||
Notes
|
||||
-----
|
||||
|
||||
* input image can be 8 bit or 16 bit with a value < 4096 (i.e. 12 bit)
|
||||
|
||||
* 8 bit images are casted in 16 bit
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from skimage import data
|
||||
>>> from skimage.morphology import disk
|
||||
>>> from skimage.filter.rank import bilateral_mean
|
||||
>>> # Load test image
|
||||
>>> ima = data.camera()
|
||||
>>> # bilateral filtering of cameraman image using a flat kernel
|
||||
>>> bilat_ima = bilateral_mean(ima, disk(20), s0=10,s1=10)
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
None, _crank16_bilateral.mean, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y,
|
||||
s0=s0, s1=s1)
|
||||
|
||||
|
||||
def bilateral_pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False, s0=10, s1=10):
|
||||
"""Return the number (population) of pixels actually inside the bilateral neighborhood,
|
||||
i.e. being inside the structuring element AND having a gray level inside the interval [g-s0,g+s1].
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, as the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
s0, s1 : int
|
||||
define the [s0,s1] interval to be considered for computing the value.
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint16 array (uint8 image are casted to uint16)
|
||||
the local number of pixels inside the bilateral neighborhood
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> # Local mean
|
||||
>>> from skimage.morphology import square
|
||||
>>> import skimage.filter.rank as rank
|
||||
>>> ima8 = 255*np.array([[0, 0, 0, 0, 0],
|
||||
... [0, 1, 1, 1, 0],
|
||||
... [0, 1, 1, 1, 0],
|
||||
... [0, 1, 1, 1, 0],
|
||||
... [0, 0, 0, 0, 0]], dtype=np.uint8)
|
||||
>>> rank.bilateral_pop(ima8, square(3), s0=10,s1=10)
|
||||
array([[3, 4, 3, 4, 3],
|
||||
[4, 4, 6, 4, 4],
|
||||
[3, 6, 9, 6, 3],
|
||||
[4, 4, 6, 4, 4],
|
||||
[3, 4, 3, 4, 3]], dtype=uint16)
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
None, _crank16_bilateral.pop, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y,
|
||||
s0=s0, s1=s1)
|
||||
|
||||
@@ -0,0 +1,77 @@
|
||||
import matplotlib.pyplot as plt
|
||||
from pprint import pprint
|
||||
|
||||
from skimage import data
|
||||
from skimage.morphology.selem import disk
|
||||
import skimage.filter.rank as rank
|
||||
|
||||
def plot_all():
|
||||
a8 = data.camera()
|
||||
a16 = a8.astype('uint16')*16
|
||||
selem = disk(5)
|
||||
|
||||
name_list = sorted([n for n in dir(rank) if n[0] is not '_'])
|
||||
print name_list
|
||||
|
||||
for n in name_list:
|
||||
if n.rfind('bilateral')==0:
|
||||
print n
|
||||
method = eval('rank.%s'%n)
|
||||
if type(method) == type(rank.maximum):
|
||||
print method
|
||||
f8 = method(a8,selem = selem,s0=10,s1=10)
|
||||
f16 = method(a16,selem = selem,s0=10,s1=10)
|
||||
plt.figure()
|
||||
plt.subplot(2,2,1)
|
||||
plt.imshow(a8)
|
||||
plt.colorbar()
|
||||
plt.subplot(2,2,2)
|
||||
plt.imshow(f8)
|
||||
plt.colorbar()
|
||||
plt.subplot(2,2,3)
|
||||
plt.imshow(f16)
|
||||
plt.colorbar()
|
||||
plt.title(method)
|
||||
for n in name_list:
|
||||
if n.rfind('percentile')==0:
|
||||
print n
|
||||
method = eval('rank.%s'%n)
|
||||
if type(method) == type(rank.maximum):
|
||||
print method
|
||||
f8 = method(a8,selem = selem,p0=.1,p1=.9)
|
||||
f16 = method(a16,selem = selem,p0=.1,p1=.9)
|
||||
plt.figure()
|
||||
plt.subplot(2,2,1)
|
||||
plt.imshow(a8)
|
||||
plt.colorbar()
|
||||
plt.subplot(2,2,2)
|
||||
plt.imshow(f8)
|
||||
plt.colorbar()
|
||||
plt.subplot(2,2,3)
|
||||
plt.imshow(f16)
|
||||
plt.colorbar()
|
||||
plt.title(method)
|
||||
for n in name_list:
|
||||
if n.find('percentile')==-1 and n.find('bilateral')==-1:
|
||||
print n
|
||||
method = eval('rank.%s'%n)
|
||||
if type(method) == type(rank.maximum):
|
||||
print method
|
||||
f8 = method(a8,selem = selem)
|
||||
f16 = method(a16,selem = selem)
|
||||
plt.figure()
|
||||
plt.subplot(2,2,1)
|
||||
plt.imshow(a8)
|
||||
plt.colorbar()
|
||||
plt.subplot(2,2,2)
|
||||
plt.imshow(f8)
|
||||
plt.colorbar()
|
||||
plt.subplot(2,2,3)
|
||||
plt.imshow(f16)
|
||||
plt.colorbar()
|
||||
plt.title(method)
|
||||
plt.show()
|
||||
|
||||
if __name__ == '__main__':
|
||||
plot_all()
|
||||
pprint(dir(rank))
|
||||
@@ -0,0 +1,32 @@
|
||||
import numpy as np
|
||||
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 denoise_bilateral
|
||||
|
||||
if __name__ == '__main__':
|
||||
a8 = data.camera()
|
||||
a16 = data.camera().astype(np.uint16)*4
|
||||
|
||||
p8 = data.page()
|
||||
|
||||
selem = disk(20)
|
||||
|
||||
otsu = rank.otsu(p8,selem)
|
||||
|
||||
|
||||
plt.figure()
|
||||
plt.subplot(1,2,1)
|
||||
plt.imshow(p8)
|
||||
plt.colorbar()
|
||||
plt.subplot(1,2,2)
|
||||
plt.imshow(otsu)
|
||||
plt.colorbar()
|
||||
plt.show()
|
||||
|
||||
|
||||
|
||||
|
||||
Binary file not shown.
@@ -0,0 +1,28 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import gdal
|
||||
|
||||
from skimage.morphology import disk
|
||||
import skimage.filter.rank as rank
|
||||
|
||||
filename = 'iko_pan_Ja1.tif'
|
||||
im16 = gdal.Open(filename).ReadAsArray().astype(np.uint16)
|
||||
|
||||
plt.figure()
|
||||
plt.imshow(im16,cmap=plt.cm.gray)
|
||||
plt.colorbar()
|
||||
|
||||
f0 = rank.median(im16,disk(1))
|
||||
f1 = rank.bilateral_mean(im16,disk(20),s0=200,s1=200)
|
||||
f2 = rank.equalize(f1,disk(10))
|
||||
f3 = rank.bottomhat(f1,disk(1))
|
||||
|
||||
plt.figure()
|
||||
plt.imshow(f2,cmap=plt.cm.gray,interpolation='nearest')
|
||||
plt.colorbar()
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,42 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from skimage import data
|
||||
from skimage.morphology.selem import disk
|
||||
import skimage.filter.rank as rank
|
||||
|
||||
print dir(rank)
|
||||
|
||||
print rank.mean
|
||||
print rank.percentile_mean
|
||||
print rank.bilateral_mean
|
||||
|
||||
a8 = data.camera()
|
||||
a16 = a8.astype('uint16')*16
|
||||
selem = disk(10)
|
||||
|
||||
f8 = rank.mean(a8,selem)
|
||||
f16 = rank.mean(a16,selem)
|
||||
|
||||
plt.figure()
|
||||
plt.imshow(np.hstack((a8,f8)))
|
||||
plt.colorbar()
|
||||
plt.figure()
|
||||
plt.imshow(np.hstack((a16,f16)))
|
||||
plt.colorbar()
|
||||
|
||||
f8 = rank.percentile_mean(a8,selem,p0=.1,p1=.9)
|
||||
f16 = rank.percentile_mean(a16,selem,p0=.1,p1=.9)
|
||||
|
||||
plt.figure()
|
||||
plt.imshow(np.hstack((a8,f8)))
|
||||
plt.colorbar()
|
||||
plt.figure()
|
||||
plt.imshow(np.hstack((a16,f16)))
|
||||
plt.colorbar()
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,11 @@
|
||||
import numpy as np
|
||||
|
||||
|
||||
def find_bitdepth(image):
|
||||
"""returns the max bith depth of a uint16 image
|
||||
"""
|
||||
umax = np.max(image)
|
||||
if umax > 2:
|
||||
return int(np.log2(umax))
|
||||
else:
|
||||
return 1
|
||||
@@ -0,0 +1,337 @@
|
||||
"""percentile_rank.py - inferior and superior ranks, provided by the user, are passed to the kernel function
|
||||
to provide a softer version of the rank filters. E.g. percentile_autolevel will stretch image levels between
|
||||
percentile [p0,p1] instead of using [min,max]. It means that isolate bright or dark pixels will not produce halos.
|
||||
|
||||
The local histogram is computed using a sliding window similar to the method described in
|
||||
|
||||
Reference: Huang, T. ,Yang, G. ; Tang, G.. "A fast two-dimensional median filtering algorithm",
|
||||
IEEE Transactions on Acoustics, Speech and Signal Processing, Feb 1979. Volume: 27 , Issue: 1, Page(s): 13 - 18.
|
||||
|
||||
input image can be 8 bit or 16 bit with a value < 4096 (i.e. 12 bit),
|
||||
for 16 bit input images, the number of histogram bins is determined from the maximum value present in the image
|
||||
|
||||
result image is 8 or 16 bit with respect to the input image
|
||||
|
||||
"""
|
||||
|
||||
from skimage import img_as_ubyte
|
||||
import numpy as np
|
||||
|
||||
from skimage.filter.rank.generic import find_bitdepth
|
||||
from skimage.filter.rank import _crank16_percentiles, _crank8_percentiles
|
||||
|
||||
__all__ = ['percentile_autolevel', 'percentile_gradient',
|
||||
'percentile_mean', 'percentile_mean_substraction',
|
||||
'percentile_morph_contr_enh', 'percentile', 'percentile_pop', 'percentile_threshold']
|
||||
|
||||
|
||||
def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y, p0, p1):
|
||||
selem = img_as_ubyte(selem)
|
||||
if mask is not None:
|
||||
mask = img_as_ubyte(mask)
|
||||
if image.dtype == np.uint8:
|
||||
return func8(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, out=out, p0=p0, p1=p1)
|
||||
elif image.dtype == np.uint16:
|
||||
bitdepth = find_bitdepth(image)
|
||||
if bitdepth > 11:
|
||||
raise ValueError("only uint16 <4096 image (12bit) supported!")
|
||||
return func16(
|
||||
image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, bitdepth=bitdepth + 1, out=out,
|
||||
p0=p0, p1=p1)
|
||||
else:
|
||||
raise TypeError("only uint8 and uint16 image supported!")
|
||||
|
||||
|
||||
def percentile_autolevel(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
|
||||
"""Return greyscale local autolevel of an image.
|
||||
|
||||
Autolevel is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, as the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
p0, p1 : float in [0.,...,1.]
|
||||
define the [p0,p1] percentile interval to be considered for computing the value.
|
||||
|
||||
Returns
|
||||
-------
|
||||
local autolevel : uint8 array or uint16 array depending on input image
|
||||
The result of the local autolevel.
|
||||
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8_percentiles.autolevel, _crank16_percentiles.autolevel, image, selem, out=out, mask=mask,
|
||||
shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
|
||||
|
||||
|
||||
def percentile_gradient(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
|
||||
"""Return greyscale local percentile_gradient of an image.
|
||||
|
||||
percentile_gradient is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, as the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
p0, p1 : float in [0.,...,1.]
|
||||
define the [p0,p1] percentile interval to be considered for computing the value.
|
||||
|
||||
Returns
|
||||
-------
|
||||
local percentile_gradient : uint8 array or uint16 array depending on input image
|
||||
The result of the local percentile_gradient.
|
||||
|
||||
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8_percentiles.gradient, _crank16_percentiles.gradient, image, selem, out=out, mask=mask,
|
||||
shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
|
||||
|
||||
|
||||
def percentile_mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
|
||||
"""Return greyscale local mean of an image.
|
||||
|
||||
Mean is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, as the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
p0, p1 : float in [0.,...,1.]
|
||||
define the [p0,p1] percentile interval to be considered for computing the value.
|
||||
|
||||
Returns
|
||||
-------
|
||||
local mean : uint8 array or uint16 array depending on input image
|
||||
The result of the local mean.
|
||||
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8_percentiles.mean, _crank16_percentiles.mean, image, selem, out=out, mask=mask,
|
||||
shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
|
||||
|
||||
|
||||
def percentile_mean_substraction(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
|
||||
"""Return greyscale local mean_substraction of an image.
|
||||
|
||||
mean_substraction is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, as the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
p0, p1 : float in [0.,...,1.]
|
||||
define the [p0,p1] percentile interval to be considered for computing the value.
|
||||
|
||||
Returns
|
||||
-------
|
||||
local mean_substraction : uint8 array or uint16 array depending on input image
|
||||
The result of the local mean_substraction.
|
||||
|
||||
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8_percentiles.mean_substraction, _crank16_percentiles.mean_substraction, image, selem, out=out,
|
||||
mask=mask, shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
|
||||
|
||||
|
||||
def percentile_morph_contr_enh(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
|
||||
"""Return greyscale local morph_contr_enh of an image.
|
||||
|
||||
morph_contr_enh is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, as the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
p0, p1 : float in [0.,...,1.]
|
||||
define the [p0,p1] percentile interval to be considered for computing the value.
|
||||
|
||||
Returns
|
||||
-------
|
||||
local morph_contr_enh : uint8 array or uint16 array depending on input image
|
||||
The result of the local morph_contr_enh.
|
||||
|
||||
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8_percentiles.morph_contr_enh, _crank16_percentiles.morph_contr_enh, image, selem, out=out,
|
||||
mask=mask, shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
|
||||
|
||||
|
||||
def percentile(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
|
||||
"""Return greyscale local percentile of an image.
|
||||
|
||||
percentile is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, as the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
p0, p1 : float in [0.,...,1.]
|
||||
define the [p0,p1] percentile interval to be considered for computing the value.
|
||||
|
||||
Returns
|
||||
-------
|
||||
local percentile : uint8 array or uint16 array depending on input image
|
||||
The result of the local percentile.
|
||||
|
||||
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8_percentiles.percentile, _crank16_percentiles.percentile, image, selem, out=out, mask=mask,
|
||||
shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
|
||||
|
||||
|
||||
def percentile_pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
|
||||
"""Return greyscale local pop of an image.
|
||||
|
||||
pop is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, as the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
p0, p1 : float in [0.,...,1.]
|
||||
define the [p0,p1] percentile interval to be considered for computing the value.
|
||||
|
||||
Returns
|
||||
-------
|
||||
local pop : uint8 array or uint16 array depending on input image
|
||||
The result of the local pop.
|
||||
|
||||
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8_percentiles.pop, _crank16_percentiles.pop, image, selem, out=out, mask=mask, shift_x=shift_x,
|
||||
shift_y=shift_y, p0=p0, p1=p1)
|
||||
|
||||
|
||||
def percentile_threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
|
||||
"""Return greyscale local threshold of an image.
|
||||
|
||||
threshold is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, as the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
p0, p1 : float in [0.,...,1.]
|
||||
define the [p0,p1] percentile interval to be considered for computing the value.
|
||||
|
||||
Returns
|
||||
-------
|
||||
local threshold : uint8 array or uint16 array depending on input image
|
||||
The result of the local threshold.
|
||||
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8_percentiles.threshold, _crank16_percentiles.threshold, image, selem, out=out, mask=mask,
|
||||
shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
|
||||
|
||||
|
||||
@@ -0,0 +1,703 @@
|
||||
"""rank.py - rankfilter for local (custom kernel) maximum, minimum, median, mean, auto-level, equalization, etc
|
||||
|
||||
The local histogram is computed using a sliding window similar to the method described in
|
||||
|
||||
Reference: Huang, T. ,Yang, G. ; Tang, G.. "A fast two-dimensional median filtering algorithm",
|
||||
IEEE Transactions on Acoustics, Speech and Signal Processing, Feb 1979. Volume: 27 , Issue: 1, Page(s): 13 - 18.
|
||||
|
||||
input image can be 8 bit or 16 bit with a value < 4096 (i.e. 12 bit),
|
||||
for 16 bit input images, the number of histogram bins is determined from the maximum value present in the image
|
||||
|
||||
result image is 8 or 16 bit with respect to the input image
|
||||
|
||||
"""
|
||||
|
||||
from skimage import img_as_ubyte
|
||||
import numpy as np
|
||||
from skimage.filter.rank import _crank8, _crank16
|
||||
|
||||
from skimage.filter.rank.generic import find_bitdepth
|
||||
|
||||
__all__ = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', 'mean', 'meansubstraction', 'median', 'minimum',
|
||||
'modal', 'morph_contr_enh', 'pop', 'threshold', 'tophat','noise_filter','entropy','otsu']
|
||||
|
||||
|
||||
def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y):
|
||||
selem = img_as_ubyte(selem)
|
||||
if mask is not None:
|
||||
mask = img_as_ubyte(mask)
|
||||
if image.dtype == np.uint8:
|
||||
if func8 is None:
|
||||
raise TypeError("not implemented for uint8 image")
|
||||
return func8(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, out=out)
|
||||
elif image.dtype == np.uint16:
|
||||
if func16 is None:
|
||||
raise TypeError("not implemented for uint16 image")
|
||||
bitdepth = find_bitdepth(image)
|
||||
if bitdepth > 11:
|
||||
raise ValueError("only uint16 <4096 image (12bit) supported!")
|
||||
return func16(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, bitdepth=bitdepth + 1, out=out)
|
||||
else:
|
||||
raise TypeError("only uint8 and uint16 image supported!")
|
||||
|
||||
|
||||
def autolevel(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Autolevel image using local histogram.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The result of the local autolevel.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from skimage import data
|
||||
>>> from skimage.morphology import disk
|
||||
>>> from skimage.filter.rank import autolevel
|
||||
>>> # Load test image
|
||||
>>> ima = data.camera()
|
||||
>>> # Stretch image contrast locally
|
||||
>>> auto = autolevel(ima, disk(20))
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8.autolevel, _crank16.autolevel, image, selem, out=out, mask=mask, shift_x=shift_x,
|
||||
shift_y=shift_y)
|
||||
|
||||
|
||||
def bottomhat(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Returns greyscale local bottomhat of an image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
local bottomhat : uint8 array or uint16 array depending on input image
|
||||
The result of the local bottomhat.
|
||||
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8.bottomhat, _crank16.bottomhat, image, selem, out=out, mask=mask, shift_x=shift_x,
|
||||
shift_y=shift_y)
|
||||
|
||||
|
||||
def equalize(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Equalize image using local histogram.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The result of the local equalize.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from skimage import data
|
||||
>>> from skimage.morphology import disk
|
||||
>>> from skimage.filter.rank import equalize
|
||||
>>> # Load test image
|
||||
>>> ima = data.camera()
|
||||
>>> # Local equalization
|
||||
>>> equ = equalize(ima, disk(20))
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8.equalize, _crank16.equalize, image, selem, out=out, mask=mask, shift_x=shift_x,
|
||||
shift_y=shift_y)
|
||||
|
||||
|
||||
def gradient(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Return greyscale local gradient of an image (i.e. local maximum - local minimum).
|
||||
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The local gradient.
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8.gradient, _crank16.gradient, image, selem, out=out, mask=mask, shift_x=shift_x,
|
||||
shift_y=shift_y)
|
||||
|
||||
|
||||
def maximum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Return greyscale local maximum of an image.
|
||||
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The local maximum.
|
||||
|
||||
See also
|
||||
--------
|
||||
skimage.morphology.dilation()
|
||||
|
||||
Note
|
||||
----
|
||||
* input image can be 8 bit or 16 bit with a value < 4096 (i.e. 12 bit)
|
||||
|
||||
* the lower algorithm complexity makes the rank.maximum() more efficient for larger images and structuring elements
|
||||
|
||||
"""
|
||||
|
||||
return _apply(_crank8.maximum, _crank16.maximum, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
|
||||
|
||||
|
||||
def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Return greyscale local mean of an image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The local mean.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from skimage import data
|
||||
>>> from skimage.morphology import disk
|
||||
>>> from skimage.filter.rank import mean
|
||||
>>> # Load test image
|
||||
>>> ima = data.camera()
|
||||
>>> # Local mean
|
||||
>>> avg = mean(ima, disk(20))
|
||||
"""
|
||||
|
||||
return _apply(_crank8.mean, _crank16.mean, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
|
||||
|
||||
|
||||
def meansubstraction(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Return image substracted from its local mean.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The result of the local meansubstraction.
|
||||
|
||||
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8.meansubstraction, _crank16.meansubstraction, image, selem, out=out, mask=mask,
|
||||
shift_x=shift_x, shift_y=shift_y)
|
||||
|
||||
|
||||
def median(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Return greyscale local median of an image.
|
||||
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The local median.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from skimage import data
|
||||
>>> from skimage.morphology import disk
|
||||
>>> from skimage.filter.rank import median
|
||||
>>> # Load test image
|
||||
>>> ima = data.camera()
|
||||
>>> # Local mean
|
||||
>>> avg = median(ima, disk(20))
|
||||
"""
|
||||
|
||||
return _apply(_crank8.median, _crank16.median, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
|
||||
|
||||
|
||||
def minimum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Return greyscale local minimum of an image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The local minimum.
|
||||
|
||||
See also
|
||||
--------
|
||||
skimage.morphology.erosion()
|
||||
|
||||
Note
|
||||
----
|
||||
* input image can be 8 bit or 16 bit with a value < 4096 (i.e. 12 bit)
|
||||
|
||||
* the lower algorithm complexity makes the rank.minimum() more efficient for larger images and structuring elements
|
||||
|
||||
"""
|
||||
|
||||
return _apply(_crank8.minimum, _crank16.minimum, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
|
||||
|
||||
|
||||
def modal(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Return greyscale local mode of an image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The local modal.
|
||||
|
||||
|
||||
"""
|
||||
|
||||
return _apply(_crank8.modal, _crank16.modal, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
|
||||
|
||||
|
||||
def morph_contr_enh(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Enhance an image replacing each pixel by the local maximum if pixel graylevel is closest to maximimum
|
||||
than local minimum OR local minimum otherwise.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The result of the local morph_contr_enh.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from skimage import data
|
||||
>>> from skimage.morphology import disk
|
||||
>>> from skimage.filter.rank import morph_contr_enh
|
||||
>>> # Load test image
|
||||
>>> ima = data.camera()
|
||||
>>> # Local mean
|
||||
>>> avg = morph_contr_enh(ima, disk(20))
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8.morph_contr_enh, _crank16.morph_contr_enh, image, selem, out=out, mask=mask, shift_x=shift_x,
|
||||
shift_y=shift_y)
|
||||
|
||||
|
||||
def pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Return the number (population) of pixels actually inside the neighborhood.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The number of pixels belonging to the neighborhood.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> # Local mean
|
||||
>>> from skimage.morphology import square
|
||||
>>> import skimage.filter.rank as rank
|
||||
>>> ima = 255*np.array([[0, 0, 0, 0, 0],
|
||||
... [0, 1, 1, 1, 0],
|
||||
... [0, 1, 1, 1, 0],
|
||||
... [0, 1, 1, 1, 0],
|
||||
... [0, 0, 0, 0, 0]], dtype=np.uint8)
|
||||
>>> rank.pop(ima, square(3))
|
||||
array([[4, 6, 6, 6, 4],
|
||||
[6, 9, 9, 9, 6],
|
||||
[6, 9, 9, 9, 6],
|
||||
[6, 9, 9, 9, 6],
|
||||
[4, 6, 6, 6, 4]], dtype=uint8)
|
||||
|
||||
|
||||
"""
|
||||
|
||||
return _apply(_crank8.pop, _crank16.pop, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
|
||||
|
||||
|
||||
def threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Return greyscale local threshold of an image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The result of the local threshold.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> # Local threshold
|
||||
>>> from skimage.morphology import square
|
||||
>>> from skimage.filter.rank import threshold
|
||||
>>> ima = 255*np.array([[0, 0, 0, 0, 0],
|
||||
... [0, 1, 1, 1, 0],
|
||||
... [0, 1, 1, 1, 0],
|
||||
... [0, 1, 1, 1, 0],
|
||||
... [0, 0, 0, 0, 0]], dtype=np.uint8)
|
||||
>>> threshold(ima, square(3))
|
||||
array([[0, 0, 0, 0, 0],
|
||||
[0, 1, 1, 1, 0],
|
||||
[0, 1, 0, 1, 0],
|
||||
[0, 1, 1, 1, 0],
|
||||
[0, 0, 0, 0, 0]], dtype=uint8)
|
||||
|
||||
|
||||
"""
|
||||
|
||||
return _apply(
|
||||
_crank8.threshold, _crank16.threshold, image, selem, out=out, mask=mask, shift_x=shift_x,
|
||||
shift_y=shift_y)
|
||||
|
||||
|
||||
def tophat(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Return greyscale local tophat of an image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The image tophat.
|
||||
|
||||
"""
|
||||
|
||||
return _apply(_crank8.tophat, _crank16.tophat, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
|
||||
|
||||
def noise_filter(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Returns the noise feature as described in [Hashimoto12]_
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's. Central element is removed during the filtering.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Reference
|
||||
----------
|
||||
|
||||
.. [Hashimoto12] N. Hashimoto et al. Referenceless image quality evaluation for whole slide imaging. J Pathol Inform 2012;3:9.
|
||||
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
The image noise .
|
||||
|
||||
"""
|
||||
# ensure that the central pixel in the structuring element is empty
|
||||
centre_r = int(selem.shape[0] / 2) + shift_y
|
||||
centre_c = int(selem.shape[1] / 2) + shift_x
|
||||
# make a local copy
|
||||
selem_cpy = selem.copy()
|
||||
selem_cpy[centre_r,centre_c] = 0
|
||||
|
||||
return _apply(_crank8.noise_filter, None, image, selem_cpy, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
|
||||
|
||||
def entropy(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Returns the entropy [wiki_entropy]_ computed locally. Entropy is computed using base 2 logarithm i.e.
|
||||
the filter returns the minimum number of bits needed to encode local greylevel distribution.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [wiki_entropy] http://en.wikipedia.org/wiki/Entropy_(information_theory)
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
entropy x10 (uint8 images) and entropy x1000 (uint16 images)
|
||||
|
||||
|
||||
Examples
|
||||
--------
|
||||
|
||||
>>> # Local entropy
|
||||
>>> from skimage import data
|
||||
>>> from skimage.filter.rank import entropy
|
||||
>>> from skimage.morphology import disk
|
||||
>>> # 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
|
||||
|
||||
"""
|
||||
|
||||
return _apply(_crank8.entropy, _crank16.entropy, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
|
||||
|
||||
def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
|
||||
"""Returns the image threshold using a the Otsu [otsu]_ locally .
|
||||
|
||||
References
|
||||
----------
|
||||
|
||||
.. [otsu] http://en.wikipedia.org/wiki/Otsu's_method
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram,
|
||||
an exception will be raised if image has a value > 4095
|
||||
selem : ndarray
|
||||
The neighborhood expressed as a 2-D array of 1's and 0's.
|
||||
out : ndarray
|
||||
If None, a new array will be allocated.
|
||||
mask : ndarray (uint8)
|
||||
Mask array that defines (>0) area of the image included in the local neighborhood.
|
||||
If None, the complete image is used (default).
|
||||
shift_x, shift_y : (int)
|
||||
Offset added to the structuring element center point.
|
||||
Shift is bounded to the structuring element sizes (center must be inside the given structuring element).
|
||||
|
||||
Returns
|
||||
-------
|
||||
out : uint8 array or uint16 array (same as input image)
|
||||
threshold image
|
||||
|
||||
|
||||
Examples
|
||||
--------
|
||||
|
||||
>>> # Local entropy
|
||||
>>> from skimage import data
|
||||
>>> from skimage.filter.rank import otsu
|
||||
>>> from skimage.morphology import disk
|
||||
>>> # defining a 8- and a 16-bit test images
|
||||
>>> a8 = data.camera()
|
||||
>>> loc_otsu = otsu(a8,disk(5))
|
||||
|
||||
"""
|
||||
|
||||
return _apply(_crank8.otsu, None, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
|
||||
@@ -0,0 +1,57 @@
|
||||
import sys
|
||||
print sys.path
|
||||
import skimage
|
||||
print skimage
|
||||
|
||||
import unittest
|
||||
|
||||
import numpy as np
|
||||
from skimage.filter import rank
|
||||
|
||||
from skimage import data
|
||||
from skimage.morphology import cmorph,disk
|
||||
from skimage.filter.rank import _crank8, _crank16
|
||||
from skimage.filter.rank import _crank16_percentiles
|
||||
|
||||
|
||||
class TestSequenceFunctions(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
pass
|
||||
|
||||
def test_trivial_selem(self):
|
||||
# check that min, max and mean returns identity if structuring element contains only central pixel
|
||||
|
||||
a = np.zeros((5,5),dtype='uint8')
|
||||
a[2,2] = 255
|
||||
a[2,3] = 128
|
||||
a[1,2] = 16
|
||||
elem = np.asarray([[0,0,0],[0,1,0],[0,0,0]],dtype='uint8')
|
||||
f = _crank8.mean(image=a,selem = elem,shift_x=0,shift_y=0)
|
||||
np.testing.assert_array_equal(a,f)
|
||||
f = _crank8.minimum(image=a,selem = elem,shift_x=0,shift_y=0)
|
||||
np.testing.assert_array_equal(a,f)
|
||||
f = _crank8.maximum(image=a,selem = elem,shift_x=0,shift_y=0)
|
||||
np.testing.assert_array_equal(a,f)
|
||||
|
||||
def test_smallest_selem(self):
|
||||
# check that min, max and mean returns identity if structuring element contains only central pixel
|
||||
|
||||
a = np.zeros((5,5),dtype='uint8')
|
||||
a[2,2] = 255
|
||||
a[2,3] = 128
|
||||
a[1,2] = 16
|
||||
elem = np.asarray([[1]],dtype='uint8')
|
||||
f = _crank8.mean(image=a,selem = elem,shift_x=0,shift_y=0)
|
||||
np.testing.assert_array_equal(a,f)
|
||||
f = _crank8.minimum(image=a,selem = elem,shift_x=0,shift_y=0)
|
||||
np.testing.assert_array_equal(a,f)
|
||||
f = _crank8.maximum(image=a,selem = elem,shift_x=0,shift_y=0)
|
||||
np.testing.assert_array_equal(a,f)
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions)
|
||||
unittest.TextTestRunner(verbosity=2).run(suite)
|
||||
@@ -0,0 +1,191 @@
|
||||
import sys
|
||||
print sys.path
|
||||
import skimage
|
||||
print skimage
|
||||
|
||||
import unittest
|
||||
|
||||
import numpy as np
|
||||
from skimage.filter import rank
|
||||
|
||||
from skimage import data
|
||||
from skimage.morphology import cmorph,disk
|
||||
from skimage.filter.rank import _crank8, _crank16
|
||||
from skimage.filter.rank import _crank16_percentiles
|
||||
|
||||
|
||||
class TestSequenceFunctions(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
pass
|
||||
|
||||
def test_random_sizes(self):
|
||||
# make sure the size is not a problem
|
||||
|
||||
niter = 10
|
||||
elem = np.asarray([[1,1,1],[1,1,1],[1,1,1]],dtype='uint8')
|
||||
for m,n in np.random.random_integers(1,100,size=(10,2)):
|
||||
a8 = np.ones((m,n),dtype='uint8')
|
||||
r = _crank8.mean(image=a8,selem = elem,shift_x=0,shift_y=0)
|
||||
self.assertTrue(a8.shape == r.shape)
|
||||
r = _crank8.mean(image=a8,selem = elem,shift_x=+1,shift_y=+1)
|
||||
self.assertTrue(a8.shape == r.shape)
|
||||
|
||||
for m,n in np.random.random_integers(1,100,size=(10,2)):
|
||||
a16 = np.ones((m,n),dtype='uint16')
|
||||
r = _crank16.mean(image=a16,selem = elem,shift_x=0,shift_y=0)
|
||||
self.assertTrue(a16.shape == r.shape)
|
||||
r = _crank16.mean(image=a16,selem = elem,shift_x=+1,shift_y=+1)
|
||||
self.assertTrue(a16.shape == r.shape)
|
||||
|
||||
for m,n in np.random.random_integers(1,100,size=(10,2)):
|
||||
a16 = np.ones((m,n),dtype='uint16')
|
||||
r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9)
|
||||
self.assertTrue(a16.shape == r.shape)
|
||||
r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=+1,shift_y=+1,p0=.1,p1=.9)
|
||||
self.assertTrue(a16.shape == r.shape)
|
||||
|
||||
def test_compare_with_cmorph_dilate(self):
|
||||
#compare the result of maximum filter with dilate
|
||||
|
||||
a = (np.random.random((500,500))*256).astype('uint8')
|
||||
|
||||
for r in range(1,20,1):
|
||||
elem = np.ones((r,r),dtype='uint8')
|
||||
# elem = (np.random.random((r,r))>.5).astype('uint8')
|
||||
rc = _crank8.maximum(image=a,selem = elem)
|
||||
cm = cmorph.dilate(image=a,selem = elem)
|
||||
self.assertTrue((rc==cm).all())
|
||||
|
||||
def test_compare_with_cmorph_erode(self):
|
||||
#compare the result of maximum filter with erode
|
||||
|
||||
a = (np.random.random((500,500))*256).astype('uint8')
|
||||
|
||||
for r in range(1,20,1):
|
||||
elem = np.ones((r,r),dtype='uint8')
|
||||
# elem = (np.random.random((r,r))>.5).astype('uint8')
|
||||
rc = _crank8.minimum(image=a,selem = elem)
|
||||
cm = cmorph.erode(image=a,selem = elem)
|
||||
self.assertTrue((rc==cm).all())
|
||||
|
||||
def test_bitdepth(self):
|
||||
# test the different bit depth for rank16
|
||||
|
||||
elem = np.ones((3,3),dtype='uint8')
|
||||
a16 = np.ones((100,100),dtype='uint16')*255
|
||||
r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=8)
|
||||
a16 = np.ones((100,100),dtype='uint16')*255*2
|
||||
r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=9)
|
||||
a16 = np.ones((100,100),dtype='uint16')*255*4
|
||||
r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=10)
|
||||
a16 = np.ones((100,100),dtype='uint16')*255*8
|
||||
r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=11)
|
||||
a16 = np.ones((100,100),dtype='uint16')*255*16
|
||||
r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=12)
|
||||
|
||||
def test_population(self):
|
||||
# check the number of valid pixels in the neighborhood
|
||||
|
||||
a = np.zeros((5,5),dtype='uint8')
|
||||
elem = np.ones((3,3),dtype='uint8')
|
||||
p = _crank8.pop(image=a,selem = elem)
|
||||
r = np.asarray([[4, 6, 6, 6, 4],
|
||||
[6, 9, 9, 9, 6],
|
||||
[6, 9, 9, 9, 6],
|
||||
[6, 9, 9, 9, 6],
|
||||
[4, 6, 6, 6, 4]])
|
||||
np.testing.assert_array_equal(r,p)
|
||||
|
||||
def test_structuring_element(self):
|
||||
# check the output for a custom structuring element
|
||||
|
||||
a = np.zeros((6,6),dtype='uint8')
|
||||
a[2,2] = 255
|
||||
elem = np.asarray([[1,1,0],[1,1,1],[0,0,1]],dtype='uint8')
|
||||
f = _crank8.maximum(image=a,selem = elem,shift_x=1,shift_y=1)
|
||||
r = np.asarray([[ 0, 0, 0, 0, 0, 0],
|
||||
[ 0, 0, 0, 0, 0, 0],
|
||||
[ 0, 0, 255, 0, 0, 0],
|
||||
[ 0, 0, 255, 255, 255, 0],
|
||||
[ 0, 0, 0, 255, 255, 0],
|
||||
[ 0, 0, 0, 0, 0, 0]])
|
||||
np.testing.assert_array_equal(r,f)
|
||||
|
||||
|
||||
@unittest.expectedFailure
|
||||
def test_fail_on_bitdepth(self):
|
||||
# should fail because data bitdepth is too high for the function
|
||||
|
||||
a16 = np.ones((100,100),dtype='uint16')*255
|
||||
elem = np.ones((3,3),dtype='uint8')
|
||||
f = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=4)
|
||||
|
||||
def test_output(self):
|
||||
#check rank function with external OUT output array
|
||||
|
||||
selem = disk(20)
|
||||
a = (np.random.random((500,500))*256).astype('uint8')
|
||||
out = np.zeros_like(a)
|
||||
f1 = rank.mean(a,selem,out=out)
|
||||
f2 = rank.mean(a,selem)
|
||||
np.testing.assert_array_equal(f1,f2)
|
||||
np.testing.assert_array_equal(out,f2)
|
||||
|
||||
@unittest.expectedFailure
|
||||
def test_inplace_output(self):
|
||||
#rank filters are not supposed to filter inplace
|
||||
|
||||
selem = disk(20)
|
||||
a = (np.random.random((500,500))*256).astype('uint8')
|
||||
out = a
|
||||
f = rank.mean(a,selem,out=out)
|
||||
np.testing.assert_array_equal(f,out)
|
||||
|
||||
|
||||
def test_compare_autolevels(self):
|
||||
# compare autolevel and percentile autolevel with p0=0.0 and p1=1.0
|
||||
# should returns the same arrays
|
||||
|
||||
image = data.camera()
|
||||
|
||||
selem = disk(20)
|
||||
loc_autolevel = rank.autolevel(image,selem=selem)
|
||||
loc_perc_autolevel = rank.percentile_autolevel(image,selem=selem,p0=.0,p1=1.)
|
||||
|
||||
assert (loc_autolevel==loc_perc_autolevel).all()
|
||||
|
||||
def test_compare_autolevels_16bit(self):
|
||||
# compare autolevel(16bit) and percentile autolevel(16bit) with p0=0.0 and p1=1.0
|
||||
# should returns the same arrays
|
||||
|
||||
image = data.camera().astype(np.uint16)*4
|
||||
|
||||
selem = disk(20)
|
||||
loc_autolevel = rank.autolevel(image,selem=selem)
|
||||
loc_perc_autolevel = rank.percentile_autolevel(image,selem=selem,p0=.0,p1=1.)
|
||||
|
||||
assert (loc_autolevel==loc_perc_autolevel).all()
|
||||
|
||||
def test_compare_8bit_vs_16bit(self):
|
||||
# filters applied on 8bit image ore 16bit image (having only real 8bit of dynamic)
|
||||
# should be identical
|
||||
|
||||
i8 = data.camera()
|
||||
i16 = i8.astype(np.uint16)
|
||||
assert (i8==i16).all()
|
||||
|
||||
methods = ['autolevel','bottomhat','equalize','gradient','maximum','mean'
|
||||
,'meansubstraction','median','minimum','modal','morph_contr_enh','pop','threshold', 'tophat']
|
||||
|
||||
for method in methods:
|
||||
func = eval('rank.%s'%method)
|
||||
f8 = func(i8,disk(3))
|
||||
f16 = func(i16,disk(3))
|
||||
assert (f8==f16).all()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions)
|
||||
unittest.TextTestRunner(verbosity=2).run(suite)
|
||||
+38
-2
@@ -14,11 +14,47 @@ def configuration(parent_package='', top_path=None):
|
||||
|
||||
cython(['_ctmf.pyx'], working_path=base_path)
|
||||
cython(['_denoise.pyx'], working_path=base_path)
|
||||
cython(['rank/_core8.pyx'], working_path=base_path)
|
||||
cython(['rank/_core16.pyx'], working_path=base_path)
|
||||
cython(['rank/_crank8.pyx'], working_path=base_path)
|
||||
cython(['rank/_crank8_percentiles.pyx'], working_path=base_path)
|
||||
cython(['rank/_crank16.pyx'], working_path=base_path)
|
||||
cython(['rank/_crank16_percentiles.pyx'], working_path=base_path)
|
||||
cython(['rank/_crank16_bilateral.pyx'], working_path=base_path)
|
||||
cython(['rank/rank.pyx'], working_path=base_path)
|
||||
cython(['rank/percentile_rank.pyx'], working_path=base_path)
|
||||
cython(['rank/bilateral_rank.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('_ctmf', sources=['_ctmf.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_denoise', sources=['_denoise.c'],
|
||||
include_dirs=[get_numpy_include_dirs(), '../_shared'])
|
||||
include_dirs=[get_numpy_include_dirs(), '../_shared'])
|
||||
config.add_extension('rank/_core8', sources=['rank/_core8.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('rank/_core16', sources=['rank/_core16.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('rank/_crank8', sources=['rank/_crank8.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension(
|
||||
'rank/_crank8_percentiles', sources=['rank/_crank8_percentiles.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('rank/_crank16', sources=['rank/_crank16.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension(
|
||||
'rank/_crank16_percentiles', sources=['rank/_crank16_percentiles.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension(
|
||||
'rank/_crank16_bilateral', sources=['rank/_crank16_bilateral.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension(
|
||||
'rank/rank', sources=['rank/rank.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension(
|
||||
'rank/percentile_rank', sources=['rank/percentile_rank.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension(
|
||||
'rank/bilateral_rank', sources=['rank/bilateral_rank.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
|
||||
Reference in New Issue
Block a user