# Source code for pysteg.features.texture

```##  \$Id\$
## -*- coding: utf-8 -*-
## (C) 2012: Hans Georg Schaathun <georg@schaathun.net>

"""
Functions to estimate the beta coefficient to quantify the
level of texture in an image.

:Module:    pysteg.features.texture
:Date:      \$Date\$
:Revision:  \$Revision\$
"""

print "[pysteg.features.texture] \$Id\$"

import numpy as np
import pywt
from math import gamma
__all__ = [ "ggdBeta" ]

[docs]def ggdBeta(im,wavelet="haar"):
"""Estimate the beta coefficient of the Generalised Gaussian
Distribution modelling the wavelet coefficients of an image im."""
(A,(H,V,D)) = pywt.dwt2(im,wavelet)
return [ ggdBetaComp(C) for C in (D,H,V) ]

def ggdBetaComp(C):
"""Estimate the beta coefficient of the Generalised Gaussian
Distribution modelling the given wavelet component C."""
return _beta(_m1(C),_m2(C))

def _m1(C):
"First moment of C."
return np.sum(np.abs(C)) / C.size
def _m2(C):
"Second moment of C."
return np.sum(C**2) / C.size
def _beta(m1,m2):
return _Finv(m1**2/m2)
def _alpha(m1,m2,beta):
return m1 * gamma(1.0/beta) / gamma(2.0/beta)

def _Finv(y,eps=10**(-5),xmin=0.025,xmax=4.0):
"""The inverse of _F(), calculated using binomial search."""
f = lambda x : _F(x) - y
(fmin,fmax) = (f(xmin),f(xmax))
print eps,xmin,xmax
while xmax-xmin > eps:
xmid = (xmin+xmax)/2.0
fmid = f(xmid)
if fmin*fmid < 0:
(xmax,fmax) = (xmid,fmid)
elif fmax*fmid < 0:
(xmin,fmin) = (xmid,fmid)
else:
raise Exception, "Neither interval contains the root."
return (xmin+xmax)/2.0
def _F(x):
return gamma(2.0/x)**2 / (gamma(3.0/x)*gamma(1.0/x))

class ggd(object):
def __init__(self,C):
(m1,m2) = (_m1(C),_m2(C))
beta = _beta(m1,m2)
self.beta = beta
alpha = _alpha(m1,m2,beta)
self.alpha = alpha
self.K = beta / ( 2*alpha*gamma(1.0/beta) )
def __call__(self,x):
return self.K * np.exp( -(np.abs(x)/self.alpha)**self.beta )

# Tests:
# 1. Compare alpha to standard deviation
# 2. Plot estimated PDF against actual histogram
```