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$
:Copyright: © 2012 Hans Georg Schaathun <georg@schaathun.net>
"""

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

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

 
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