%matplotlib inline
from IPython.display import HTML,Image,SVG,YouTubeVideo
Texture features¶
Here is an example of a texture found in ultrasound medical imaging, the image exhibits a clear striation pattern due to muscle fiber orientation.
Image('../data/muscle.jpg')
The following video sequience illustrates how fiber orientation can vary during the motricity exercice.
from IPython.display import YouTubeVideo
YouTubeVideo('PUcz11MLxUk', start=0, autoplay=1, theme="light", color="blue",)
Here are some classical textures extracted from an atlas (Broadatz). Texture can be natural, such as grains, sand, soil, stone, biological tissue or resulting from a human artifact such as fabric, wall, pavements ...
Texture in an image can be characterized by a repetitive pattern or by a statistical geometrical distribution.
Image('https://s3-eu-west-1.amazonaws.com/ppreviews-plos-725668748/569742/preview.jpg')
Texture is relative to a region content, it can by due to a specific gray-level dynamic, a gray-level organisation one define several oreder for texture description:
- first order : the pixel alone is used ti compute the texture descriptor
- second order : the pixel and one of its neighboor are used
- third order : more than 1 neighboor are used
The coocurence matrix¶
see also:
- the coocurence matrix in the digital image section
Coorence matrix can be defined such as:
$$ C^{\Delta x, \Delta y}_{i,j}=\sum_{p=1}^n\sum_{q=1}^m\begin{cases} 1, & \mbox{if }I(p,q)=i\mbox{ and }I(p+\Delta x,q+\Delta y)=j \\ 0, &\mbox{otherwise}\end{cases} $$for a given $(\Delta x, \Delta y)$ defining the neighborhood relationship.
From the coocurence matrix $C$ one can extract several descriptors that are sensitive to the texture:
Normalized coocurence matrix:
$$ P_{i,j} = \frac{C_{i,j}}{\sum^{N-1}_{i,j=0} C_{i,j}}$$contrast measures
$$Contrast = \sum^{N-1}_{i,j=0} P_{i,j}(i-j)^2$$$$Dissimilarity = \sum^{N-1}_{i,j=0} P_{i,j}|i-j|$$$$Homogeneity = \sum^{N-1}_{i,j=0} \frac{P_{i,j}}{1+(i+j)^2}$$organization measures
$$ASM = \sum^{N-1}_{i,j=0} P_{i,j}^2$$$$ Energy = \sqrt{ASM}$$$$entropy = - \sum^{N-1}_{i,j=0} P_{i,j} \, \log(P_{i,j})$$statistical measures
$$ \mu_i = \sum^{N-1}_{i,j=0} i\,P_{i,j} $$$$ \mu_j = \sum^{N-1}_{i,j=0} j\,P_{i,j} $$$$ \sigma_i = \sum^{N-1}_{i,j=0} (i-\mu_i)^2\,P_{i,j} $$$$ \sigma_j = \sum^{N-1}_{i,j=0} (j-\mu_j)^2\,P_{i,j} $$$$Correlation = \sum^{N-1}_{i,j=0} P_{i,j} \, \frac{(i-\mu_i)(j-\mu_j)}{\sqrt{\sigma_i^2 \sigma_j^2}}$$The Gabor filters¶
Gabor filters are a family of 2D filter kernels based on the product of a cosine(a sine for imaginary part) function and a gaussian function such that:
for the real part:
$$g(x,y;\lambda,\theta,\psi,\sigma,\gamma) = \exp\left(-\frac{x'^2+\gamma^2y'^2}{2\sigma^2}\right)\cos\left(2\pi\frac{x'}{\lambda}+\psi\right)$$and for the imaginary part:
$$g(x,y;\lambda,\theta,\psi,\sigma,\gamma) = \exp\left(-\frac{x'^2+\gamma^2y'^2}{2\sigma^2}\right)\sin\left(2\pi\frac{x'}{\lambda}+\psi\right)$$where
$$x' = x \cos\theta + y \sin\theta$$and
$$y' = -x \sin\theta + y \cos\theta$$The Gabor parameters are:
- $\lambda$ the wave length
- $\theta$ the orientation
- $\sigma$ the gaussian envelope
- $\gamma$ the form factor
- $\psi$ the phase
import numpy as np
from scipy.ndimage.filters import convolve
import matplotlib.pyplot as plt
def gabor(sigma,theta,lam,psi,gamma,X,Y):
Xt = X * np.cos(theta) + Y * np.sin(theta)
Yt = -X * np.sin(theta) + Y * np.cos(theta)
n = Xt**2. + gamma**2. * Yt**2.
d = 2.0 * sigma**2
E1 = np.exp(- n/d)
a = 2.0*np.pi* Xt/lam + psi
E2 = np.exp(1j * a)
return E1 * E2
def getGaborKernel(ksize,sigma,theta,lambd,gamma,psi):
X,Y = np.meshgrid(np.linspace(-1.,1.,ksize[0]),np.linspace(-1.,1.,ksize[1]))
G = gabor(sigma=sigma,theta=theta,lam=lambd,psi=psi,gamma=gamma,X=X,Y=Y)
return G
ksize = [31,31]
theta = 0.
lambd = 10.
gamma = .15
psi = np.pi
plt.figure();
fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(6, 6))
for i,sigma in enumerate(np.linspace(.1,.5,4)):
kernel = getGaborKernel(ksize,sigma,theta,lambd,gamma,psi)
axes[i].imshow(np.abs(kernel),cmap=plt.cm.gray);
axes[i].set_xlabel('sigma = %.2f'%sigma)
plt.figure();
sigma = .1
fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(6, 6))
for i,theta in enumerate(np.linspace(.1,.5,4)):
kernel = getGaborKernel(ksize,sigma,theta,lambd,gamma,psi)
axes[i].imshow(np.abs(kernel),cmap=plt.cm.gray)
axes[i].set_xlabel('theta = %.2f'%theta)
plt.figure();
theta = 0
fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(6, 6))
for i,gamma in enumerate(np.linspace(.1,1,4)):
kernel = getGaborKernel(ksize,sigma,theta,lambd,gamma,psi)
axes[i].imshow(np.abs(kernel),cmap=plt.cm.gray)
axes[i].set_xlabel('gamma = %.2f'%gamma)
plt.show()
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
Hurst coefficient¶
Hurst coefficient are related to the local fractal dimension, the algorithm is the following:
- compute local max and min (typ. on a 7 pixel octogonal neighbourhood)
- plot the log (difference) vs log (distance)
- least square fit the curve
| pixel | # | distance |
|---|---|---|
| a | 1 | 0 |
| b | 4 | 1 |
| c | 4 | 1.41 |
| d | 4 | 2 |
| e | 8 | 2.236 |
| f | 4 | 2.828 |
| g | 4 | 3 |
| h | 8 | 3.162 |
| h | g | h | ||||
| f | e | d | e | f | ||
| h | e | c | b | c | e | h |
| g | d | b | a | b | d | g |
| h | e | c | b | c | e | h |
| f | e | d | e | f | ||
| h | g | h |
Fractal analysis¶
Similarly we can analyse the value of the image filtered by a gaussian filter in a log/log plot of the pixel value with respect to the sigma of the filter. In the following example we compute the pixel gray level in five points located on the central line, the gray level is evaluated after applying a gaussian filter of increasing sigma.
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as npy
from skimage.io import imread
from skimage.data import camera
from scipy import ndimage,interpolate
from scipy.ndimage.filters import convolve,gaussian_filter
from os.path import join
from mpl_toolkits.mplot3d import Axes3D
im = camera()[-1::-2,::2].astype(float)
sigmas = npy.asarray([.5,1.,1.5,2.,2.5])
C = im.copy().astype(float)
C = C[:,:,npy.newaxis]
for s in sigmas:
c = gaussian_filter(im,s)
C = npy.dstack((C,c[:,:,npy.newaxis]))
sample_n = 120
samples_m = [230,186,146,72,20]
log_scale = npy.log(sigmas)
plt.figure(figsize=[5,5])
for m in samples_m:
value = npy.ravel(C[m,sample_n,1:])
log_value = npy.log(value)
plt.plot(log_scale,log_value,label='y=%d'%m)
plt.legend(loc=3)
plt.figure(figsize=[5,5])
plt.subplot(2,2,1)
plt.imshow(im,interpolation='nearest',cmap=cm.gray,origin='lower')
plt.title('original')
plt.vlines(sample_n,0,256,color=[.8,0,0])
for m in samples_m:
plt.hlines(m,sample_n-5,sample_n+5,color=[0,0.8,0])
plt.subplot(2,2,2)
plt.imshow(C[:,sample_n,:],interpolation='nearest',cmap=cm.gray,origin='lower',extent=[0,256,0,256])
plt.show()
Some location have a constant value with respect to the gaussian radius, suach as the pixel located in the sky, it means that this part of the image is very homogeneous and insensitive to gaussian filtering, on the opposite, a pixel located in the tripod ($y=72$) exhibits a variable profile, because it is situated into a part of the image having more high frequencies, therefore more sensitive to the gaussian filter.
see also:
- Texture [DIP] p506
- Fractal analysis: Hurst coeficient [IPH] p262,208
- Spectral approach [DIP] p511
- Gabor filters [HCVA] vol2 p71
- Edge-flow [Ma00] vol2 p71
- ... [DIPCASA] p185