Grey-Level Cooccurrence Matrix ( GLCM )

GLCM_Ignore_Background.py

# -*- coding: utf-8 -*-
"""
Created on Wed Jul 04 11:07:04 2018

@author: Baikal

E-mail: lxianbnu1@yeah.net
If you have any question or some advice about the program,please don't hesitate to contact me.

Introduction:
    Due to GLCM calculation procedure is going on every pixels of the whole image and Python processing efficiency is
    not very well, so the program is a little time-consuming.
"""

# Reference:
#   https://www.harrisgeospatial.com/docs/backgroundtexturemetrics.html (Texture Metrics Background---Harris Forum)
#    

#Purpose:To calculate GLCM ignore background.

import os 
import gdal
import numpy as np
import time
from collections import Counter

def Grey_Scale_Convert( S1_DIR, Grey_Scale ):
    
    HH_file = 'Aim_Beta0_HH_Correction_RB.tif'
    S1_HH = os.path.join( S1_DIR, HH_file )
    S1_HH_Image = gdal.Open( S1_HH )
    S1_HH_Data = S1_HH_Image.ReadAsArray(  )
    ## Remove data value of 255 pixels.
    S1_HH_Data_Remove0 = S1_HH_Data[ S1_HH_Data != 255 ]
    max = S1_HH_Data_Remove0.max( )
    min = S1_HH_Data_Remove0.min( )
    Xsize = S1_HH_Image.RasterXSize # Columns
    Ysize = S1_HH_Image.RasterYSize # Rows
    Geotransform = S1_HH_Image.GetGeoTransform(  )
    Projection = S1_HH_Image.GetProjection(  )
    
    ##Create Output file.
    Driver = gdal.GetDriverByName( "GTiff" )
    S1_Scale_Data = np.zeros( ( Ysize , Xsize ), dtype = np.int8 )
    out_file = 'grey_level_convention.tif'
    Out_S1_HH = os.path.join( S1_DIR, out_file )
    if os.path.exists( Out_S1_HH ):
        os.remove( Out_S1_HH )
    Output_S1_HH_Tiff = Driver.Create( Out_S1_HH ,Xsize ,Ysize ,1 ,gdal.GDT_Int16 )
    Output_S1_HH_Tiff.SetGeoTransform( Geotransform )
    Output_S1_HH_Tiff.SetProjection( Projection )

    for row in range( Ysize ):
        for col in range( Xsize ):
            if S1_HH_Data[ row, col ] != 255:
                S1_Scale_Data[ row, col ] = ( S1_HH_Data[ row, col ] * 1.0 - min ) / ( max - min ) * int( Grey_Scale )

    Output_S1_HH_Tiff.GetRasterBand(1).WriteArray( S1_Scale_Data )
    Output_S1_HH_Tiff.FlushCache() #Save to disk
    S1_HH_Image = None #Close it.
    Output_S1_HH_Tiff = None #Close the file.
    
    return Out_S1_HH
            
def GLCM( Out_S1_HH, Grey_Scale, Window_Size, Step, S1_DIR ):
    

    Start_time = time.time()
    S1_HH_glc_Image = gdal.Open( Out_S1_HH )
    S1_HH_glc_Data = S1_HH_glc_Image.ReadAsArray(  )    
    Xsize = S1_HH_glc_Image.RasterXSize # Columns
    Ysize = S1_HH_glc_Image.RasterYSize # Rows    
    Geotransform = S1_HH_glc_Image.GetGeoTransform(  )
    Projection = S1_HH_glc_Image.GetProjection(  ) 
    # To create Textural Features TIFF files.
    Driver = gdal.GetDriverByName( "GTiff" )
    Output_Classification_Dir =  S1_DIR + r'/' + Grey_Scale + '_' + Window_Size + '_' + Step + '/'
    if not os.path.exists( Output_Classification_Dir ):
        os.mkdir( Output_Classification_Dir )
    Output_Classification_Dir_Name = Output_Classification_Dir + r'Aim_Beta0_HH_Correction_Statistics_' + Grey_Scale + '_' + Window_Size + '_' + Step + '.tif'
    
    Grey_Scale = int( Grey_Scale )
    Window_Size = int( Window_Size )
    Step = int( Step )
    
    if os.path.exists( Output_Classification_Dir_Name ):
        os.remove( Output_Classification_Dir_Name )   
    Output_Tiff = Driver.Create( Output_Classification_Dir_Name ,Xsize ,Ysize ,10 ,gdal.GDT_Float32 )
    Output_Tiff.SetGeoTransform( Geotransform )
    Output_Tiff.SetProjection( Projection )
    
    # Define Reference Matrix.
    # Define Target Matrix.
    ## Initializing.  
    Mean = np.zeros( [ Ysize, Xsize ] )
    Variance = np.zeros( [ Ysize, Xsize ] )
    Energy = np.zeros( [ Ysize, Xsize ] )
    Homogeneity = np.zeros( [ Ysize, Xsize ] )
    Contrast = np.zeros( [ Ysize, Xsize ] )
    Entropy = np.zeros( [ Ysize, Xsize ] )
    Skewness = np.zeros( [ Ysize, Xsize ] )
    Kurtosis = np.zeros( [ Ysize, Xsize ] )
    Correlation = np.zeros( [ Ysize, Xsize ] )
    Cluster_Prominence = np.zeros( [ Ysize, Xsize ] )    
    
    ## Few parameters
    Half = ( Window_Size - 1 ) / 2
    Prob = np.zeros( [ Grey_Scale, Grey_Scale ] ) ## Grey Level Co-occurrence Matrix
    
    Dif_row_col = np.zeros( [ Grey_Scale, Grey_Scale ] )
    Mul_row_col = np.zeros( [ Grey_Scale, Grey_Scale ] )
    Add_row_col = np.zeros( [ Grey_Scale, Grey_Scale ] )
    
    # Compute Intermediate Variable.
    for row in range( Grey_Scale ):
        for col in range( Grey_Scale ):
            Dif_row_col[ row, col ] = row - col
            Mul_row_col[ row, col ] = ( row + 1 ) * ( col + 1 )
            Add_row_col[ row, col ] = ( row + 1 ) + ( col + 1 )
    

    
    # Calculate GLCM.
    for row in range( Ysize ):
        for col in range( Xsize ):
            # Control Overflow. 0 degree
            if S1_HH_glc_Data[ row, col ] != 0  \
            and ( row - Step - Half ) > 0 and ( col - Step - Half ) > 0 \
            and ( row + Step + Half ) < ( Ysize ) and ( col + Step + Half ) < ( Xsize ): 
               
                M_0 = np.zeros( [ Grey_Scale, Grey_Scale ] )
#                M_45 = np.zeros( [ Grey_Scale, Grey_Scale ] )
#                M_90 = np.zeros( [ Grey_Scale, Grey_Scale ] )
#                M_135 = np.zeros( [ Grey_Scale, Grey_Scale ] )        
                
                # Define Source Window.
                Ref_Value_Matrix = S1_HH_glc_Data[ row - Half : row + Half + 1, col - Half : col + Half + 1 ]
                Ref_Str_Matrix = Ref_Value_Matrix.astype( str )
                # Define Target Window.
                Tar_0_Value_Matrix = S1_HH_glc_Data[ row + Step - Half : row + Step + Half + 1, col + Step - Half  : col + Step + Half + 1 ]
                Tar_0_Str_Matrix = Tar_0_Value_Matrix.astype( str )
                # Connect Source and Target Window.
                Ref_Tar_Matrix = np.core.defchararray.add( np.core.defchararray.add( Ref_Str_Matrix ,'_' ), Tar_0_Str_Matrix )
                # To one-dimensional.
                Ref_Tar_int_Array  = Ref_Tar_Matrix.reshape( Window_Size * Window_Size )
                Ref_Tar_list_Matrix = Ref_Tar_int_Array.tolist()
                result = Counter( Ref_Tar_list_Matrix )
                
                ##
                for pixel_value, num in result.items():
                    ref = int( pixel_value[ : pixel_value.find( '_' ) ] )
                    tar = int( pixel_value[ pixel_value.find( '_' ) + 1 : ] )
                    M_0[ ref - 1, tar - 1 ] = num
                    
                if M_0.sum() == 0:
                    Prob_0 = np.zeros( [ Grey_Scale, Grey_Scale ] )
                else:
                    Prob_0 = M_0 / ( M_0.sum() * 1.0 )
#                Prob_45 = M_0 / ( M_45.sum() * 1.0 )
#                Prob_90 = M_0 / ( M_90.sum() * 1.0 )
#                Prob_135 = M_0 / ( M_135.sum() * 1.0 )
#                Prob = ( Prob_0 + Prob_45 + Prob_90 + Prob_135 ) / 4
                Prob = Prob_0
                
                # Calculate textural features.
                Energy[ row, col ] = np.square( Prob ).sum(  )
                Arbitrarily_Number = 2.2e-16
                Entropy[ row, col ] = - ( Prob * np.log( Prob + Arbitrarily_Number ) ).sum()

   
                Range = np.arange( Grey_Scale )
                
                Mean_X = ( ( Range + 1 ) * Prob.sum( axis = 1 ) ).sum( )# axis = 1 means rows 
                Mean_Y = ( ( Range + 1 ) * Prob.sum( axis = 0 ) ).sum( )# axis = 0 means columns
                Stan_X = np.sqrt( ( np.square( Range + 1 - Mean_X ) * Prob.sum( axis = 1 ) ).sum( ) )
                Stan_Y = np.sqrt( ( np.square( Range + 1 - Mean_Y ) * Prob.sum( axis = 0 ) ).sum( ) )                
                
                Mean[ row, col ] = ( ( Range + 1 ) * Prob.sum( axis = 1 ) ).sum()
                Variance[ row, col ] = ( np.square( Range + 1 - Mean[ row, col ] ) * Prob.sum( axis = 1 ) ).sum()
                if ( Variance[ row, col ] ).sum() == 0:
                    Skewness[ row, col ] = 0
                    Kurtosis[ row, col ] = 0
                else:
                    Skewness[ row, col ] = ( ( Prob - Mean[ row, col ] ) ** 3 / ( Variance[ row, col ] * np.sqrt( Variance[ row, col ] ) ) ).sum()
                    Kurtosis[ row, col ] = ( ( Prob - Mean[ row, col ] ) ** 4 / ( Variance[ row, col ] ** 2 ) ).sum()                                        
                Homogeneity[ row, col ] = ( Prob / ( 1 + np.square( Dif_row_col ) ) ).sum()                
                Contrast[ row, col ] = ( np.square( Dif_row_col ) * Prob ).sum()
                Correlation[ row, col ] = ( ( Mul_row_col * Prob ).sum() - Mean_X * Mean_Y ) / ( Stan_X * Stan_Y )
                Cluster_Prominence[ row, col ] = ( ( Add_row_col - Mean_X - Mean_Y ) ** 4 * Prob ).sum()
                
        print 'Have finished ' + str( row * 1.0 / Ysize * 100 ) + ' %'
    
    # Write Classification Files.
    Output_Tiff.GetRasterBand( 1 ).WriteArray( Mean )
    Output_Tiff.GetRasterBand( 2 ).WriteArray( Variance )
    Output_Tiff.GetRasterBand( 3 ).WriteArray( Homogeneity )
    Output_Tiff.GetRasterBand( 4 ).WriteArray( Contrast )
    Output_Tiff.GetRasterBand( 5 ).WriteArray( Entropy )
    Output_Tiff.GetRasterBand( 6 ).WriteArray( Energy )
    Output_Tiff.GetRasterBand( 7 ).WriteArray( Correlation )
    Output_Tiff.GetRasterBand( 8 ).WriteArray( Skewness )
    Output_Tiff.GetRasterBand( 9 ).WriteArray( Kurtosis )
    Output_Tiff.GetRasterBand( 10 ).WriteArray( Cluster_Prominence )
    Output_Tiff.FlushCache() #Save to disk
    S1_HH_glc_Image = None #Close it.
    Output_Tiff = None #Close the file.

    End_time = time.time()
    Using_time = End_time - Start_time    
    print 'Time:' + str( "%.*f" % (3, Using_time / 60 ) ) + 'min' + str( "%.*f" % (3, Using_time % 60 ) ) + 's'
    print 'Classification TIFF has written down\r\n'



### Main
IN_ROOT_DIR = r'G:\201606Sen\scihub.copernicus.eu\201619/'
Direction = 0

GLCM_Paramaters_Folder =                                     \
                       [                                     \
#                          '16_07_03','16_07_05',             \
#                          '16_15_03','16_15_05','16_15_07',  \
#                          '32_07_03','32_07_05',             \
                          '32_79_03'
#                          '32_15_03','32_15_05','32_15_07',  \
#                          '64_07_03','64_07_05',             \
#                          '64_15_03','64_15_05','64_15_07'   \
                       ]
                       
#GLCM_Paramaters_Folder =                                     \
#                       [                                     \
#                          '64_07_05',
#                       ]                       

files = os.listdir( IN_ROOT_DIR )
for file in files:
    if os.path.isdir( os.path.join( IN_ROOT_DIR, file ) ):
        
        for GLCM_ID in range( len( GLCM_Paramaters_Folder ) ): 
            
            GLCM_Pars = GLCM_Paramaters_Folder[ GLCM_ID ]
            Grey_Scale = GLCM_Pars[ :2 ]
            Window_Size = GLCM_Pars[ 3:5 ]
            Step = GLCM_Pars[ 6:8 ]
            
            FOLDER = file
            S1_DIR = os.path.join( IN_ROOT_DIR, FOLDER, 'Composition' )
            Out_S1_HH = Grey_Scale_Convert( S1_DIR, Grey_Scale )
            print 'Grey Scale Convention has down'
            GLCM( Out_S1_HH, Grey_Scale, Window_Size, Step, S1_DIR  )