back to iatexture

Function histstat

Synopse

The histstat function computes the statistics from the image histogram.

  • h,stats = histstat(f, mask=[])
    • Output
      • h: histogram
      • stats: array with histogram' statistics
    • Input
      • f: input image
      • mask: define the mask

Description

The first block of histstat function computes the histogram from a specific image. So, the second block

extracts some relevant statistics about the image from the histogram that was found before. These statistics are

called Mazda's attributes and they are:

  • Mean;
  • Variance;
  • Skewness;
  • Kurtosis;
  • Perc.01%;
  • Perc.10%;
  • Perc.50%;
  • Perc.90%;
  • Perc.99%;
  • Entropy;
  • Mode;
  • Median;

Function Code

 1 import ia636 as ia
 2 import numpy as np
 3 
 4 def histstat(f,mask=[]):
 5 
 6 
 7     ### compute histogram ###
 8 
 9     if mask !=[]:
10         h = ia.iahistogram(f[mask])
11     else:
12         h = ia.iahistogram(f)
13     hn = 1.0*h/h.sum()                                                                   # compute the normalized image histogram
14     stats = np.zeros(12)                                                                 # number of statistics
15 
16 
17     ### compute statistics ###
18 
19     n = len(h) # number of gray values
20     stats[0]  = np.sum((np.arange(n)*hn))                                                # mean
21     stats[1]  = np.sum(np.power((np.arange(n)-stats[0]),2)*hn)                           # variance
22     stats[2]  = np.sum(np.power((np.arange(n)-stats[0]),3)*hn)/(np.power(stats[1],1.5))  # skewness
23     stats[3]  = np.sum(np.power((np.arange(n)-stats[0]),4)*hn)/(np.power(stats[1],2))-3  # kurtosis
24     stats[4:9]  = np.round(ia.iah2percentile(h,np.array([1,10,50,90,99])))               # 1,10,50,90,99% percentile
25 
26 
27     ## extra attributes ##
28 
29     stats[9] =  -1*(hn*np.log10(hn+np.spacing(1))).sum()                                 # entropy
30     stats[10] = np.argmax(hn)                                                            # mode
31     stats[11] = np.where(np.cumsum(hn) >= 0.5)[0][0]                                     # median
32 
33     if  stats[1] == 0:                                                                   # if std equal to zero, both skewness and kurtosis are zero
34         stats[2] = 0
35         stats[3] = 0
36 
37     return h,stats

Examples

Numeric Example

 1 import ia636 as ia
 2 import numpy as np
 3 import scipy
 4 
 5 from iatexture.histstat import histstat
 6 
 7 f = np.array([1,1,1,0,1,2,2,2,1])
 8 h,stats = histstat(f)
 9 print 'histogram',h
10 print
11 print 'statistics =',stats
12 print scipy.stats.entropy(f)
histogram [1 5 3]

statistics = [ 1.22222222  0.39506173 -0.20992233 -0.62109375  0.          1.          1.
  2.          2.          0.40688542  1.          1.        ]
2.01981499249

Numeric Example using a mask

 1 import numpy as np
 2 import ia636 as ia
 3 
 4 from iatexture.histstat import histstat
 5 
 6 f = np.array([1,1,1,0,1,2,2,2,1])
 7 mask = f>0
 8 h,stats = histstat(f,mask)
 9 print 'histogram',h
10 print
11 print 'statistics =',stats)
------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 11
    print 'statistics =',stats)
                              ^
SyntaxError: invalid syntax

------------------------------------------------------------

Numeric Example 3D

 1 import numpy as np
 2 import ia636 as ia
 3 
 4 from iatexture.histstat import histstat
 5 
 6 f = np.array([1,1,1,0,1,2,2,2,1,0,0,0,1,1,1,1,2,0,1,1,2,1,0,0]).reshape(2,3,4)
 7 h,stats = histstat(f)
 8 print 'histogram',h
 9 print
10 print 'statistics =', stats
histogram [ 7 12  5]

statistics = [ 0.91666667  0.49305556  0.11700664 -0.97242611  0.          0.          1.
  2.          2.          0.44851494  1.          1.        ]

Numeric Example 3D

 1 import numpy as np
 2 import ia636 as ia
 3 
 4 from iatexture.histstat import histstat
 5 
 6 f = array([[[1,2,2,0,0,1],[0,0,1,2,2,1],[1,1,0,0,0,2],[1,1,1,2,2,2],[1,1,2,2,0,0]],
 7     [[1,2,2,0,0,1],[0,0,1,2,2,1],[1,1,0,0,0,2],[1,1,1,2,2,2],[1,1,2,2,0,0]],
 8     [[1,2,2,0,0,1],[0,0,1,2,2,1],[1,1,0,0,0,2],[1,1,1,2,2,2],[1,1,2,2,0,0]]], dtype=uint8)
 9 
10 h,stats = histstat(f)
11 print 'histogram',h
12 print
13 print 'statistics =', stats
histogram [27 33 30]

statistics = [ 1.03333333  0.63222222 -0.05953097 -1.41606308  0.          0.          1.
  2.          2.          0.47567118  1.          1.        ]

Numeric Example 3D using a mask

 1 import numpy as np
 2 import ia636 as ia
 3 
 4 from iatexture.histstat import histstat
 5 
 6 f = np.array([1,1,1,0,1,2,2,2,1,0,0,0,1,1,1,1,2,0,1,1,2,1,0,0]).reshape(2,3,4)
 7 h,stats = histstat(f,mask = f>0)
 8 print 'histogram',h
 9 print
10 print 'statistics =', stats
histogram [ 0 12  5]

statistics = [ 1.29411765  0.20761246  0.90369611 -1.18333333  1.          1.          1.
  2.          2.          0.26309451  1.          1.        ]

Image Example

 1 import numpy as np
 2 import ia636 as ia
 3 
 4 from iatexture.histstat import histstat
 5 
 6 f = adread('cameraman.tif')
 7 adshow(f,'input image')
 8 h,stats = histstat(f)
 9 
10 print 'histogram',h
11 print
12 print 'statistics =',stats
histogram [   1  280 1229  952 1216  829  824  722  758  817  980  838 1172 1208 1353
 1123  921  648  502  319  265  228  157  118  103   95   83   68   76   61
   81   73   65   79   66   58   72   81   55   47   66   65   72   74   63
   76   70   73   82   78   73   82   86   71   64   80   72   49   54   49
   46   54   51   39   39   49   45   37   37   53   35   47   43   42   56
   42   61   54   54   52   40   45   49   48   37   49   52   45   60   67
   68   48   49   49   52   51   47   43   38   24   34   33   35   31   31
   49   31   31   35   32   24   41   36   39   40   46   38   31   41   38
   46   49   55   31   45   48   59   47   65   71   68   67   57   83   82
   95   96   97   96   97   95  104  126  121  109  139  144  133  148  153
  147  160  150  184  157  175  178  183  171  203  208  192  231  234  228
  225  268  268  279  285  278  298  274  317  319  346  339  308  348  362
  363  373  380  367  383  407  360  330  439  471  539  545  584  644  647
  603  702  682  702  680  775  721  689  628  702  754  739  681  812  807
  856  813  875  824  811  742  744  597  618  465  570  562  572  452  494
  479  442  343  353  297  274  223  235  220  235  204  164  149  153  124
  100   90   92   69   79   68   60   62   57   60   48  776]

statistics = [  1.37065933e+02   7.60425522e+03  -5.54976386e-01  -1.41520647e+00
   2.00000000e+00   8.00000000e+00   1.80000000e+02   2.22000000e+02
   2.51000000e+02   2.16964814e+00   1.40000000e+01   1.80000000e+02]

input image

Equation

The follow equations represents the statistics equations, where N denotes the number of intensity levels in the imagem ( usually between 0 and 255) and is the probability density function, i.e. the normalized histogram (h) of the image, and F is a random variable (any image) that represents the pixels intensities .

Mean

Variance

Skewness

Kurtosis

Percentile

The k-th percentile is the value that corresponds to the cumulative frequency of , where N is the amostral size

Entropy

Mode

The value that appears most often.

Median

It is the numerical value separating the higher half of a data sample or a probability distribution.

Contributions

  • Mariana Bento, August 2013: initial code
  • Mariana Leite, September 2014: 3D examples

Memory test

 1 import multiprocessing as mp
 2 import numpy as np
 3 import psutil
 4 from ia636 import iacontour,iameshgrid,iaroi,iagshow
 5 import glob, urllib, os
 6 
 7 foldername = '/awmedia/www/media/p/LesionMRI/CAINdata/Original_Data/CAIN10320013/BASELINE - CAIN1/FLAIR_LongTR_601/*'
 8 folder = glob.glob(foldername)
 9 folder.sort()
10 
11 img = zeros((len(folder),560,560),dtype=uint8)
12 
13 for i in arange(len(folder)):
14      aux = dicom.read_file(find_attachment_file(folder[i]))
15      img[i] = aux.pixel_array.astype(float32)
16 
17 # retalhar a imagem em retangulos
18 sizex = 30
19 sizey = 30
20 sizez = 30
21 
22 z = arange(0,img.shape[0],sizez)
23 x = arange(0,img.shape[1],sizex)
24 y = arange(0,img.shape[2],sizey)
25 print 'initial memory use:  %9.3f MB' % (psutil.phymem_usage().used/1e06,)
26 for k in arange(z.shape[0]-1):
27     for i in arange(x.shape[0]-1):
28         for j in arange(y.shape[0]-1):
29             roi = np.zeros(img.shape, dtype = bool)
30             roi[z[k]:z[k+1],x[i]:x[i+1],y[j]:y[j+1]] =1
31             print 'memory used before:  %9.3f MB' % (psutil.phymem_usage().used/1e06,)
32             att_rec = histstat(img,mask = roi)
33             print 'memory used after:  %9.3f MB' % (psutil.phymem_usage().used/1e06,)
34 
35 
36 
37 print 'final memory use:  %9.3f MB' % (psutil.phymem_usage().used/1e06,)
------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 14, in <module>
NameError: name 'dicom' is not defined

------------------------------------------------------------