back to iatexture

Function gradstat

Synopse

The gradstat function computes statistics based on the image gradient

  • g = gradimg(f)
    • Output
      • g: ndarray, 2D, gradient image
    • Input
      • f: ndarray, 2D, input image
  • g = gradstat(f, mask=[])
    • Output
      • g: gradient statistics
    • Input
      • f: ndarray, 2D, input image
      • mask: define the mask

Description

The * gradimg* function computes the image gradient, while the gradstat extracts some relevant statistics from the image gradient.

Function Code

001. from scipy.stats import *
002. from ia636 import ianormalize
003. from ia636 import iacontour
004. 
005. def gradimg(f):
006.     from numpy import *
007.     import ia636 as ia
008. 
009.     if len(f.shape) == 2:                                    # 2D case
010.         h1 = array([[0,1,0],                                 # Defining the horizontal mask
011.                    [0,0,0],
012.                    [0,-1,0]])
013.         h2 = array([[0,0,0],                                 # Defining the vertical mask
014.                    [1,0,-1],
015.                    [0,0,0]])
016. 
017.         aux1 = ia.iaconv(f,h1)[1:-1,1:-1].astype(int)        # Make the convolution between horizontal mask and image
018.         aux2 = ia.iaconv(f,h2)[1:-1,1:-1].astype(int)        # Make the convolution between vertical mask and image
019.         g = sqrt(aux1**2 + aux2**2)                       # Use the equation to compute the gradient of an image
020. 
021.         return g
022. 
023.     else:                                                    # 3D case
024.         h1 = array([[[0,0,0],                                # Defining the horizontal mask
025.                      [0,0,0],
026.                      [0,0,0]],
027.                     [[0,1,0],
028.                      [0,0,0],
029.                      [0,-1,0]],
030.                     [[0,0,0],
031.                      [0,0,0],
032.                      [0,0,0]]])
033.         h2 = array([[[0,0,0],                                # Defining the vertical mask
034.                      [0,0,0],
035.                      [0,0,0]],
036.                     [[0,0,0],
037.                      [1,0,-1],
038.                      [0,0,0]],
039.                     [[0,0,0],
040.                      [0,0,0],
041.                      [0,0,0]]])
042.         h3 = array([[[0,0,0],                                # Defining the depth mask
043.                      [0,1,0],
044.                      [0,0,0]],
045.                     [[0,0,0],
046.                      [0,0,0],
047.                      [0,0,0]],
048.                     [[0,0,0],
049.                      [0,-1,0],
050.                      [0,0,0]]])
051. 
052.         aux1 = ia.iaconv(f,h1)[1:-1,1:-1,1:-1].astype(int)    # Make the convolution between horizontal mask and image
053.         aux2 = ia.iaconv(f,h2)[1:-1,1:-1,1:-1].astype(int)    # Make the convolution between vertical mask and image
054.         aux3 = ia.iaconv(f,h3)[1:-1,1:-1,1:-1].astype(int)    # Make the convolution between depth mask and image
055.         grad = sqrt(aux1**2 + aux2**2 + aux3**2)              # Use the equation to compute the gradient of an image
056.         return grad
057. 
058. 
059. def gradstat(f,mask=[]):
060.     import ia636 as ia
061.     import ia870
062.     from numpy import *
063.     import morph
064.     from scipy.stats import skew,kurtosis
065. 
066.     grad = gradimg(f)
067.     b = ia870.iasecross(1)
068. 
069.     if len(f.shape)==3:                                                   #3D case
070. 
071.         if mask!=[]:
072.             for i in arange(mask.shape[0]):
073.                 mask_aux = zeros((mask.shape[1]+2,mask.shape[2]+2))
074.                 mask_aux[1:-1,1:-1] = mask[i]                             # zero bounder
075.                 mask[i] = ia870.iaero(mask_aux,b)[1:-1,1:-1]              # extract the outside border
076.                 #mask[i] = mask[i] - ia.iacontour(mask[i])
077.             grad = grad[mask>0]
078. 
079.         else:
080.             m = grad>-1 # all pixels
081.             for i in arange(m.shape[0]):
082.                 m_aux = zeros((m.shape[1]+2,m.shape[2]+2))
083.                 m_aux[1:-1,1:-1] = m[i]                                   # zero bounder
084.                 m[i] = ia870.iaero(m_aux,b)[1:-1,1:-1]                    # extract the outside border
085.                 #m[i] = m[i] - ia.iacontour(m[i])                         # extract the outside border
086.             grad = grad[m>0]
087. 
088. 
089.     else:
090.         if mask!=[]:
091.             mask_aux = zeros((mask.shape[0]+2,mask.shape[1]+2))
092.             mask_aux[1:-1,1:-1] = mask                                    # zero bounder
093.             mask = ia870.iaero(mask_aux,b)[1:-1,1:-1]                     # extract the outside border
094.             #mask = mask - ia.iacontour(mask)
095.             grad = grad[mask>0]
096. 
097.         else:
098.             m = grad>-1                                                   # all pixels
099.             m_aux = zeros((m.shape[0]+2,m.shape[1]+2))
100.             m_aux[1:-1,1:-1] = m                                          # zero bounder
101.             m = ia870.iaero(m_aux,b)[1:-1,1:-1]                           # extract the outside border
102.             #m = m - ia.iacontour(m)                                      # extract the outside border
103.             grad = grad[m>0]
104. 
105. 
106.     M = size(grad)
107.     percNonZeros = 1.0*count_nonzero(grad)/M
108.     return grad,[mean(grad), var(grad), skew(ravel(grad)), kurtosis(ravel(grad)), percNonZeros]

Examples

Numerical example:

01. from gradstat import gradstat,gradimg
02. from iatexture import histstat
03. 
04. f = array( [[2,2,0,1,1,1,0,0,0],
05.             [1,2,1,1,1,1,1,1,1],
06.             [1,1,1,1,1,1,1,1,1],
07.             [0,0,2,2,2,1,1,0,0]], dtype=uint8)
08. 
09. g = gradimg(f)
10. set_printoptions(precision=2)
11. print 'Etapa 1: gradiente usando as mascaras acima:\n',g
12. grad,stats = gradstat(f)
13. print 'Etapa 2: gradiente - pixels do contorno externo:\n', grad
14. print
15. print
16. print
Etapa 1: gradiente usando as mascaras acima:
[[ 2.24  2.83  1.41  1.41  1.    1.41  1.41  1.    1.  ]
 [ 2.24  1.    1.41  0.    0.    0.    1.    1.    1.41]
 [ 1.41  2.    1.    1.    1.    0.    0.    1.    1.41]
 [ 1.    2.24  2.24  1.    1.41  1.41  1.41  1.41  1.  ]]
Etapa 2: gradiente - pixels do contorno externo:
[ 1.    1.41  0.    0.    0.    1.    1.    2.    1.    1.    1.    0.    0.
  1.  ]

Numerical example:

01. from gradstat import gradstat,gradimg
02. from iatexture import histstat
03. 
04. f = array( [[0,0,0,1,1,1,0,0,0],
05.             [1,1,1,1,1,1,1,1,1],
06.             [1,1,1,1,1,1,1,1,1],
07.             [0,0,0,1,1,1,1,0,0]], dtype=uint8)
08. 
09. g = gradimg(f)
10. print 'gradient[f]=',g
11. print 'g: ', g[f>0]
12. print
13. _,gstat = gradstat(f,f>0)
14. print 'mean absolute gradient',gstat[0]
15. print 'variance absolute gradient',gstat[1]
16. print 'skewness absolute gradient',gstat[2]
17. print 'kurtosis absolute gradient',gstat[3]
18. print 'percentage of pixels with nonzero gradient',gstat[4]
19. 
20. print
21. print 'g: ', g
22. _,gstat = gradstat(f)
23. print 'mean absolute gradient',gstat[0]
24. print 'variance absolute gradient',gstat[1]
25. print 'skewness absolute gradient',gstat[2]
26. print 'kurtosis absolute gradient',gstat[3]
27. print 'percentage of pixels with nonzero gradient',gstat[4]
gradient[f]= [[ 1.    1.    1.41  1.41  1.    1.41  1.41  1.    1.  ]
 [ 1.41  1.    1.    0.    0.    0.    1.    1.    1.41]
 [ 1.41  1.    1.    0.    0.    0.    0.    1.    1.41]
 [ 1.    1.    1.41  1.41  1.    1.    1.41  1.41  1.  ]]
g:  [ 1.41  1.    1.41  1.41  1.    1.    0.    0.    0.    1.    1.    1.41
  1.41  1.    1.    0.    0.    0.    0.    1.    1.41  1.41  1.    1.
  1.41]

mean absolute gradient 0.0
variance absolute gradient 0.0
skewness absolute gradient 0.0
kurtosis absolute gradient -3.0
percentage of pixels with nonzero gradient 0.0

g:  [[ 1.    1.    1.41  1.41  1.    1.41  1.41  1.    1.  ]
 [ 1.41  1.    1.    0.    0.    0.    1.    1.    1.41]
 [ 1.41  1.    1.    0.    0.    0.    0.    1.    1.41]
 [ 1.    1.    1.41  1.41  1.    1.    1.41  1.41  1.  ]]
mean absolute gradient 0.5
variance absolute gradient 0.25
skewness absolute gradient 0.0
kurtosis absolute gradient -2.0
percentage of pixels with nonzero gradient 0.5

Numerical example 2:

01. from gradstat import gradstat,gradimg
02. 
03. f = array( [[0,0,0,1,1],
04.             [1,1,1,1,1],
05.             [1,1,1,1,1],
06.             [0,0,0,1,1],
07.             [1,1,1,1,1]], dtype=uint8)
08. 
09. adshow(f,'f')
10. g = gradimg(f)
11. 
12. print histstat(f)
13. print
14. print 'gradient[f]=',g
15. 
16. _,gstat = gradstat(f)
17. print 'mean absolute gradient',gstat[0]
18. print 'variance absolute gradient',gstat[1]
19. print 'skewness absolute gradient',gstat[2]
20. print 'kurtosis absolute gradient',gstat[3]
21. print 'percentage of pixels with nonzero gradient',gstat[4]
(array([ 6, 19]), array([ 0.76,  0.18, -1.22, -0.52,  0.  ,  0.  ,  1.  ,  1.  ,  1.  ,
        0.24,  1.  ,  1.  ]))

gradient[f]= [[ 1.    1.    1.41  1.41  1.41]
 [ 1.41  1.    1.    0.    1.  ]
 [ 1.41  1.    1.    0.    1.  ]
 [ 0.    0.    1.    1.    1.  ]
 [ 1.    0.    0.    1.    1.41]]
mean absolute gradient 0.666666666667
variance absolute gradient 0.222222222222
skewness absolute gradient -0.707106781187
kurtosis absolute gradient -1.5
percentage of pixels with nonzero gradient 0.666666666667

f

3D examples

Example 1

01. from gradstat import gradstat,gradimg
02. from iatexture import histstat
03. import ia636 as ia
04. 
05. f = array( [[[0,0,0,1,1],
06.             [1,1,1,1,1],
07.             [1,1,1,1,1],
08.             [0,0,0,1,1]],
09.             [[0,0,0,1,1],
10.             [1,1,1,1,1],
11.             [1,1,1,1,1],
12.             [0,0,0,1,1]]], dtype=uint8)
13. print 'input image: \n',f
14. g = gradimg(f)
15. print 'gradiente[f]=',g
16. _,gstat = gradstat(f)
17. m = g>-1 # all pixels
18. for i in arange(m.shape[0]):
19.     m[i] = m[i] - ia.iacontour(m[i]) # extract the outside border
20.     grad = g[m]
21. print grad
22. print 'mean absolute gradient',gstat[0]
23. print 'variance absolute gradient',gstat[1]
24. print 'skewness absolute gradient',gstat[2]
25. print 'kurtosis absolute gradient',gstat[3]
26. print 'percentage of pixels with nonzero gradient',gstat[4]
input image: 
[[[0 0 0 1 1]
  [1 1 1 1 1]
  [1 1 1 1 1]
  [0 0 0 1 1]]

 [[0 0 0 1 1]
  [1 1 1 1 1]
  [1 1 1 1 1]
  [0 0 0 1 1]]]
gradiente[f]= [[[ 1.    1.    1.41  1.73  1.73]
  [ 1.73  1.41  1.41  1.    1.41]
  [ 1.73  1.41  1.41  1.    1.41]
  [ 1.    1.    1.41  1.73  1.73]]

 [[ 1.    1.    1.41  1.73  1.73]
  [ 1.73  1.41  1.41  1.    1.41]
  [ 1.73  1.41  1.41  1.    1.41]
  [ 1.    1.    1.41  1.73  1.73]]]
[ 1.41  1.41  1.    1.41  1.41  1.    1.41  1.41  1.    1.41  1.41  1.  ]
mean absolute gradient 1.27614237492
variance absolute gradient 0.038127305612
skewness absolute gradient -0.707106781187
kurtosis absolute gradient -1.5
percentage of pixels with nonzero gradient 1.0
1. from iatexture import *
2. import ia636
3. 
4. f = adread('cameraman.tif')
5. g = gradimg(f)
6. adshow(ia636.ianormalize(f),'original image')
7. adshow(ia636.ianormalize(g),'gradient image')

original image

gradient image

Example 2

1. 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]],
2.     [[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]],
3.     [[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)
4. _,gstat = gradstat(f)
5. print 'grad stats',gstat
grad stats [1.6462436205539284, 0.81765971956327144, -0.3824762471527743, -0.7156357814566463, 0.8611111111111112]

Equations

Let be the gradient matrix computed with 3x3 neighbor.

Mean absolute gradient

Variance of absolute gradient

Skewness of absolute gradient

Kurtosis of absolute gradient

References

See Also

Contributions

  • Mariana, 24oct2013: initial function.
  • Mariana, 19set2014: 3d implementation.

Memory test

01. import multiprocessing as mp
02. import numpy as np
03. import psutil
04. from ia636 import iacontour,iameshgrid,iaroi,iagshow
05. import glob, urllib, os
06. from glcmdesc import glcmdesc
07. import dicom
08. 
09. foldername = '/awmedia/www/media/p/LesionMRI/CAINdata/Original_Data/CAIN10320013/BASELINE - CAIN1/FLAIR_LongTR_601/*'
10. folder = glob.glob(foldername)
11. folder.sort()
12. 
13. img = zeros((len(folder),560,560),dtype=uint8)
14. 
15. for i in arange(len(folder)):
16.      aux = dicom.read_file(find_attachment_file(folder[i]))
17.      img[i] = aux.pixel_array.astype(float32)
18. 
19. # retalhar a imagem em retangulos
20. sizex = 50
21. sizey = 50
22. sizez = 20
23. 
24. z = arange(0,img.shape[0],sizez)
25. x = arange(0,img.shape[1],sizex)
26. y = arange(0,img.shape[2],sizey)
27. print 'initial memory use:  %9.3f MB' % (psutil.phymem_usage().used/1e06,)
28. for k in arange(z.shape[0]/2):
29.     for i in arange(x.shape[0]/2):
30.         for j in arange(y.shape[0]/2):
31.             roi = np.zeros(img.shape, dtype = bool)
32.             roi[z[k]:z[k+1],x[i]:x[i+1],y[j]:y[j+1]] =1
33.             print 'memory used before:  %9.3f MB' % (psutil.phymem_usage().used/1e06,)
34.             desc = gradstat(img,mask = roi)
35.             print 'memory used after:  %9.3f MB' % (psutil.phymem_usage().used/1e06,)
36. 
37. 
38. 
39. print 'final memory use:  %9.3f MB' % (psutil.phymem_usage().used/1e06,)
initial memory use:   1114.984 MB
memory used before:   1114.984 MB
memory used after:   1076.732 MB
memory used before:   1070.764 MB
memory used after:   1076.986 MB
memory used before:   1091.715 MB
memory used after:   1092.731 MB
memory used before:   1074.954 MB
memory used after:   1018.585 MB
memory used before:   1033.568 MB
memory used after:   1052.877 MB
memory used before:   1037.386 MB
memory used after:   1018.696 MB
memory used before:   1033.679 MB
memory used after:   1034.441 MB
memory used before:   1018.696 MB
memory used after:   1018.823 MB
memory used before:   1033.552 MB
memory used after:   1034.441 MB
memory used before:   1018.696 MB
memory used after:   1018.696 MB
memory used before:   1033.679 MB
memory used after:   1034.568 MB
memory used before:   1018.823 MB
memory used after:   1018.696 MB
memory used before:   1033.679 MB
memory used after:   1034.441 MB
memory used before:   1018.696 MB
memory used after:   1018.696 MB
memory used before:   1033.679 MB
memory used after:   1034.695 MB
memory used before:   1018.950 MB
memory used after:   1018.823 MB
memory used before:   1033.806 MB
memory used after:   1034.568 MB
memory used before:   1018.823 MB
memory used after:   1018.696 MB
memory used before:   1033.679 MB
memory used after:   1034.441 MB
memory used before:   1018.696 MB
memory used after:   1018.696 MB
memory used before:   1033.679 MB
memory used after:   1034.568 MB
memory used before:   1018.823 MB
memory used after:   1018.696 MB
memory used before:   1033.679 MB
memory used after:   1059.099 MB
memory used before:   1043.100 MB
memory used after:   1078.772 MB
memory used before:   1093.755 MB
memory used after:   1034.330 MB
memory used before:   1018.585 MB
memory used after:   1084.772 MB
memory used before:   1099.502 MB
memory used after:   1100.517 MB
memory used before:   1084.772 MB
memory used after:   1088.836 MB
memory used before:   1103.819 MB
memory used after:   1106.485 MB
memory used before:   1090.740 MB
memory used after:   1094.169 MB
memory used before:   1109.152 MB
memory used after:   1110.168 MB
memory used before:   1094.423 MB
memory used after:   1094.930 MB
memory used before:   1109.914 MB
memory used after:   1116.135 MB
memory used before:   1100.390 MB
memory used after:   1101.279 MB
memory used before:   1116.008 MB
memory used after:   1117.278 MB
memory used before:   1101.533 MB
memory used after:   1109.914 MB
final memory use:   1109.914 MB