SCM / SIM / MIDE

__init__.py

# Structural Cooccurrence Matrix - SCM
# OBS.: implementation of SCM for image analysis
# Geraldo Ramalho v0.8 ago/2016
#
# please cite: 
#    Ramalho et al. Rotation-invariant Feature Extraction using a Structural Co-occurrence Matrix. 
#    Measurement, v.94, p.406-415, dec/2016.
#    doi:10.1016/j.measurement.2016.08.012
#    lapisco.ifce.edu.br/?page_id=191
#
# In this version the SCM parameters k, P and W are matrices passed by the user.
# Default values for these parameters are provided.
# Some default partition functions Q are provided as well, but user can chose any other partition function.
# Weight matrix W is introduced since it was already used in MATLAB version of the code.
#
# known issues:
# - matrix and attributes values may vary sligthly from MATLAB version because of the differences in the partition function used
#
# status/to do:
# - 
#from scmclass import *
__all__=["scmclass"]

scm_example.py

########################################################################
### SCM EXAMPLE
### This code shows how to use scm_v??? module.
###
### Geraldo Ramalho - ago/2016
### Code corrections: 
###    1) a bug could make the matrix sum less than all pixels in the input image - Daniel Silva ([email protected]) - 4/jul/2019
### 
### please cite: 
###    Ramalho et al. Rotation-invariant Feature Extraction using a Structural Co-occurrence Matrix. 
###    Measurement, v.94, p.406-415, dec/2016.
###    doi:10.1016/j.measurement.2016.08.012
###    lapisco.ifce.edu.br/?page_id=191
########################################################################
from skimage import io, color # self-dependencies: Cython 2.3
import numpy as np
import time
#from skimage.filters import gaussian
from scm_v009_8 import *
#import scipy.ndimage
#from scipy.misc.pilutil import Image

### first input signal (f)
#print("Read image...")
#f = io.imread("2013-10-27 (1).jpg")

load_samples = np.loadtxt("filteredFeaturesPump1data.txt", delimiter = " ")#carrega o arquivo
#f = (np.transpose(load_samples[0:,1:2]))[0,:] #entropy
f = (np.transpose(load_samples[0:,0:1]))[0,:] #energy

#if f.ndim == 1: #if f is a vector, converts it into matrix
#    f=np.expand_dims(f,axis=0)

# print "Capture image..."
# from SimpleCV import Camera # dependencias: Pillow, pygame
# # fnitialize the camera
# cam = Camera()
# # Loop to continuously get images
# while True:
#     # Get fmage from camera
#     img = cam.getfmage()
#     # Make image black and white
#     img = img.binarize()
#     # Draw the text "Hello World" on image
#     img.drawText("Hello World!")
#     # Show the image
#     img.show()
# cam.getfmage().show()
# exit()
'''
# adjust first input signal
f = color.rgb2gray(f)
f = f.astype(float)
f = (f - f.min())/(f.max() - f.min()) # skimage already in [0,1]
'''
 # Pega somente a saída e transforma em um vetor

### second input signal (g)
#g = io.imread("teste.png")
g = f;
#g = (np.transpose(load_samples[0:,1:2]))[0,:] #entropy
# g = f + 3.6 * f.std() * np.random.random(f.shape) 
# g = np.clip(g, 0, 1)
# g = gaussian(f, sigma = 1, mode = 'reflect')
# g = mean(f, disk(1)) # usar disk apenas em morfologia (1 = 3x3)

# performing convolution - average filter
# g = scipy.ndimage.filters.convolve(f, np.ones((3,3))/9)

# adjust second input signal
#g = color.rgb2gray(g)
#g = g.astype(float)
#g = (g - g.min())/(g.max() - g.min()) # skimage already in [0,1]

# crop - used only to speedy tests and to make debug easier
# f = f[0:5,0:5] 
# g = g[0:5,0:5] 
'''
### set user callback functions
def cb_Q_quant(img, num_levels, *args):
    print(args)
    img = (img - img.min())/(img.max() - img.min())
    return np.around(img * (num_levels - 1))

def cb_k_average(img, *args):
    print(args)
    return scipy.ndimage.filters.convolve(img, np.ones((3,3))/9)

def cb_P_absdif(img1, img2, num_levels):
    # Code correction #1:
    #   return abs(img1 - img2) < num_levels
    return abs(img1.astype(float) - img2.astype(float)) < num_levels

def cb_W_uniform(num_levels):
    tmp = np.ones((num_levels, num_levels))
    return (tmp + 1)/num_levels
'''
### Compute SCM
#print("Computing SCM for image analysis...")
ini = time.time()
SCM_PARAM = {'N': 8,'k':'LowPass','Verbose': False}#, 'W_callback': cb_W_uniform}
#SCM_PARAM = {'N': 8,'k':'HighPass','Verbose': False}#, 'W_callback': cb_W_uniform}

M, p = scmclass.SCMimg().scm(f, g, SCM_PARAM)
'''print("SCM attributes:")
print("  group I  :\n  cor = {:0.4f} idm = {:0.4f}".format(p[0], p[1]))
print("  group II :\n  ent = {:0.4f}".format(p[2]))
print("  group III:\n  csd = {:0.4f} csr = {:0.4f} mdr = {:0.4f} dkl = {:0.4f} cad = {:0.4f}".format(p[3], p[4], p[5], p[6], p[7]))
print("  group NEW:\n  net = {:0.4f} mut = {:0.4f}".format(p[8], p[9]))
#print M
print("DONE!!!")
print("TOTAL TIME = {:0.4f}".format(time.time() - ini))
print("")
'''
print(M)
print("Grupo I:\nCOR = {:0.4f}; IDM = {:0.4f}".format(p[0], p[1]))
print("Grupo II:\nENT = {:0.4f}".format(p[2]))
print("Grupo III:\nCSD = {:0.4f}; CSR = {:0.4f}; MDR = {:0.4f}; DLK = {:0.4f}; CAD = {:0.4f}".format(p[3], p[4], p[5], p[6], p[7]))
'''
#from scipy.misc import toimage
#im = toimage(outImage)
#im.save("outImage6.png")
from matplotlib import pyplot as plt
plt.imshow(outImage, interpolation='nearest')
'''

io.imshow(M)
io.show()

scmclass.py

########################################################################
### SCM framework for image analysis
### input	- signals: f and g
### output - matrix: m_(i,j)=#{(i,j)|P(i,j), i=Q(f)_p, j=Q(k(g))_p+d}
###        - attributes:	cor, idm, ent, csd, csr, mdr, dkl, cad
###
### please cite: 
###    Ramalho et al. Rotation-invariant Feature Extraction using a Structural Co-occurrence Matrix. 
###    Measurement, v.94, p.406-415, dec/2016.
###    doi:10.1016/j.measurement.2016.08.012
###    lapisco.ifce.edu.br/?page_id=191
###
### Code correction: 
###    1) a bug could make the matrix sum less than all pixels in the input image - Daniel Silva ([email protected]) - 4/jul/2019
###
### default SCM PARAMETERS
### PAR_SCM_N = 8								# N	= num of levels (partitions)
### PAR_SCM_Q = "Quantization" 					# Q = partition function: Quantization,...
### PAR_SCM_k = scm_default_k('none')				# k = convolutional filter: none, Low Pass, High Pass, Soft Low Pass...
### PAR_SCM_d = (0,0)							# d = row and col displacement of g
### PAR_SCM_P = np.array([])						# P	= must be computed for Q(f) and Q(k(g))
### PAR_SCM_W = scm_default_W(PAR_SCM_N)			# W	= matrix of weigths
########################################################################

### SCM DEPENDENCIES
import numpy as np
import scipy.ndimage
import os
import csv
from skimage import io, color
from skimage.restoration import denoise_tv_chambolle#, denoise_bilateral #Total Variation
from skimage.morphology import disk #used for Average Quantization
from skimage.filters.rank import mean #used for Average Quantization

### SCM PROCEDURES
class SCMimg(object):

    # scm_default_W: default weigth matrix
    def scm_default_W(self, N):
        tmp = np.zeros((N,N))
        for i in range(N):
            for j in range(N):
                tmp[i,j] = abs(i - j);
        return (tmp + 1)/N

    # scm_default_P: default property matrix
    def scm_default_P(self, img1, img2, NL):
        # Code correction #1:
        #   return (abs(img1 - img2)) < NL
        return abs(img1.astype(float) - img2.astype(float)) < NL

    # scm_default_k: filtering of the second input signal
    def scm_default_k(self, img, *args):
        if args[0] == 'High Pass': # High Pass = 3x3 laplacian filter with alpha = 0.2
            filter = np.matrix([[0.1667, 0.6667, 0.1667], [0.6667, -3.3333, 0.6667], [0.1667, 0.6667, 0.1667]]) #original line
            #filter = np.asarray([[0.1667, 0.6667, 0.1667], [0.6667, -3.3333, 0.6667], [0.1667, 0.6667, 0.1667]]) #gscm_map
        elif args[0] == 'Low Pass': # Low Pass = 3x3 average filter
            filter = np.ones((3,3))/9
        elif args[0] == 'Soft Low Pass': # Soft Low Pass = 3x3 gaussian filter with sigma = 0.5 #gscm_map
            filter = np.matrix([[0.0113, 0.0838, 0.0113], [0.0838, 0.6193, 0.0838], [0.0113, 0.0838, 0.0113]]) #original line
            #filter = np.asarray([[0.0113, 0.0838, 0.0113], [0.0838, 0.6193, 0.0838], [0.0113, 0.0838, 0.0113]])
        else:
            filter = np.array([])
        if filter.size > 0:
            """
            #for gscm_map
            if img.shape[1] == 3:
                filter = filter.reshape((1,3,3)) #reshape to calculate the convolution over the windows
            """
            tmp = scipy.ndimage.filters.convolve(img, filter)
            
            """
            if (tmp.max() - tmp.min() == 0):
                tmp = np.ones(img) #for gscm_map
            #else:
            #    tmp = (tmp - tmp.min())/(tmp.max() - tmp.min())
            """
            tmp = (tmp - tmp.min())/(tmp.max() - tmp.min())
            return tmp
        else:
            return img

    # scm_Q - Q: partitioning function
    def scm_default_Q(self, img, *args):
        numl = args[0]
        if args[1] == "Quantization":
            img = (img - img.min())/(img.max() - img.min())
            return np.around(img * (numl - 1))
        elif args[1] == "Total Variation":
            img = denoise_tv_chambolle(img, weight = 0.1, multichannel = True)
            img = (img - img.min())/(img.max() - img.min())
            return np.round(img * (numl - 1))
        elif args[1] == "Average Quantization":
            img = mean(img, disk(3))
            img = (img - img.min())/(img.max() - img.min())
            return np.round(img * (numl - 1))
        else:
            return img

    # scm_d - d: displacement of the second input signal
    def scm_d(self, img, d):
        if sum(d) > 0:
            tmp = np.roll(img, d[0], axis = 0)
            return np.roll(tmp, d[1], axis = 1)
        else:
            return img

    # scm_m - m: assymetric matrix of frequencies of matching pairs
    def scm_M(self, img1, img2, prop):
        numl = np.floor(img1.max()) + 1 # original line
        #numl = np.floor(img2.max()) + 1
        if prop.size == 0: # default property
            prop = SCMimg.scm_default_P(img1 - img2)
        #(cols,rows) = np.meshgrid(list(range(img1.shape[1])), list(range(img1.shape[0])))
        (cols,rows) = np.meshgrid(range(img1.shape[1]), range(img1.shape[0]))
        rows = rows[np.where(prop)]
        cols = cols[np.where(prop)]
        
        i = np.asarray(img1[rows,cols])
        j = np.asarray(img2[rows,cols])
        h = np.concatenate(([i],[j]), axis = 0)
        h = h.astype(int)
        hs = h[0,:] + h[1,:] * 255
        
        ua, uind = np.unique(hs, return_inverse = True)
        uacount = np.bincount(uind)
        m = np.zeros((int(numl), int(numl)), dtype = np.uint32) #warning (non-integer)
        i = 0
        for n in ua:
            n1 = n//255
            n0 = n - n1 * 255
            m[n0,n1] = uacount[i]
            i += 1
        return m

    # scm_props - computes attributes from SCM
    def scm_props(self, m, w):
        m = m * w
        m = m.astype(float)/m.sum()
        mv = np.squeeze(m) # vectorize m
        mvnz = mv[np.where(mv > 0)] # m with non-zero values

        # group I (cor, idm) and group II (ent)
 
        # implementation of SCM as in MATLAB version
        # row and col subscripts => pixel values in the SCM
        (c,r) = np.meshgrid(list(range(m.shape[0])), list(range(m.shape[0])))
        c = np.squeeze(c) + 1 # vectorize c
        r = np.squeeze(r) + 1 # vectorize r

        # marginal probabilities
        P = np.sum(m, axis = 0) # sum of rows
        P /= P.sum() # marginal probabilities of the rows
        Q = np.sum(m, axis = 1) # sum of cols
        Q /= Q.sum() # marginal probabilities of the cols
        Pnz = P[np.where(P > 0)]
        Qnz = Q[np.where(Q > 0)]

        # mr = mean of rows, mc = mean of cols
        # sr = std of rows, sc = std of cols
        mr = r * np.squeeze(m); mr = mr.sum()
        sr = np.square(r - mr) * mv; sr = np.sqrt(sr.sum())
        mc = c * np.squeeze(m); mc = mc.sum()
        sc = np.square(c - mc) * mv; sc = np.sqrt(sc.sum())

        # COR
       
        term1 = (r - mr) * (c - mc) * mv
        term2 = term1.sum()
        #decision to avoid NaN
        if term2 == 0 or (sr * sc) == 0:
            cor = 0.0
        else:
            cor = term2/(sr * sc)

        # IDM
        term1 = (1 + abs(r - c))
        term2 = mv/term1
        idm = term2.sum()

        # ENT
        term1 = mvnz * np.log(mvnz)
        ent = -term1.sum()

        # CSD
        O = m.diagonal()
        E = np.sum(m, axis = 1)
        term2 = O - E
        term1 = np.square(term2[np.where(E > 0)])/E[np.where(E > 0)]
        csd = term1.sum()

        # CSR
        mv2 = mv + 1e-6
        term1 = mv2[0: m.shape[0]//2, 0: m.shape[1]//2]; term1 = np.squeeze(term1)
        term2 = mv2[m.shape[0]//2:, m.shape[1]//2:]; term2 = np.squeeze(term2)
        term1 /= term1.sum()
        term2 /= term2.sum()
        term2 = (term1 + term2)/2
        term1 = np.square(term1 - term2)/term2
        csr = term1.sum()

        # MDR
        idxp = np.squeeze(np.argwhere(P > 0))
        idxq = np.squeeze(np.argwhere(Q > 0))
        coefp = 0
        coefq = 0
        if len(idxp.shape) > 0:
            for i in idxp:
                for j in idxp:
                    coefp += P[i] * abs(i - j)
        if len(idxq.shape) > 0:
            for i in idxq:
                for j in idxq:
                    coefq += Q[i] * abs(i - j)
        mdr = 0
        if (coefp > 0) & (coefq > 0):
            mdr = coefp/coefq
            if mdr > 1:
                mdr = coefq/coefp

        # DKL
        idx = (P > 0) & (Q > 0)
        term1 = P[np.where(idx)]/Q[np.where(idx)]
        term2 = P[np.where(idx)]
        term1 = np.log(term1) * term2;
        dkl = term1.sum()

        # CAD
        term1 = abs(P - Q)
        cad = 1 - term1.sum()

        # NET
        term1 = mvnz * np.log(mvnz)/mv.shape[0]
        net = -term1.sum()

        # MUT
        term = mvnz * np.log2(mvnz) # joint entropy
        term1 = Pnz * np.log2(Pnz)
        term2 = Qnz * np.log2(Qnz)
        mut = (-term1.sum()) + (-term2.sum()) - (-term.sum());

        return (cor, idm, ent, csd, csr, mdr, dkl, cad, net, mut)
	
    
    # scm - computes matrix and attributes given f, g and parameters
    def scm(self, f, g, par):
        
       #if f and/or g are vectors, converts them into matrix
        if f.ndim == 1:
            f=np.expand_dims(f,axis=0)
        if g.ndim == 1:
            g=np.expand_dims(g,axis=0)
        
        #
        if 'Verbose' in par:
            verbose = par['Verbose']
        else:
            verbose = False
        #
        if 'N' in par:
            N = par['N']
        else:
            N = 8
        #
        if 'd' in par:
            d = par['d']
        else:
            d = (0,0)
        #
        if 'k' in par:
            k_name = par['k']
        else:
            k_name = 'none'
        #
        if 'k_callback' in par:
            kcallback = par['k_callback']
        else:
            kcallback = self.scm_default_k
        #
        if 'Q' in par:
            Q_name = par['Q']
        else:
            Q_name = 'Quantization'
        #
        if 'Q_callback' in par:
            Qcallback = par['Q_callback']
        else:
            Qcallback = self.scm_default_Q
        #	
        if 'P_callback' in par:
            Pcallback = par['P_callback']
        else:
            Pcallback = self.scm_default_P
        #	
        if 'W_callback' in par:
            Wcallback = par['W_callback']
        else:
            Wcallback = self.scm_default_W
        
            

        # W: weigth matrix	
        W = Wcallback(N)

        # Q: partitioning Q(f)
        Qf = Qcallback(f, N, Q_name)
        

        # Qkgd: partitioning Q(k(g))_p+d
        Qkgd = self.scm_d(Qcallback(kcallback(g, k_name), N, Q_name), d)

        # P: similarity property between pixels P(i,j), i = Q(f)_p, j = Q(k(g))_p+d
        P = Pcallback(Qf, Qkgd, N)

        # M: computing matrix elements m(i,j)=#{(i,j)|P(i,j), i = Q(f)_p, j = Q(k(g))_p+d}
        M = self.scm_M(Qf, Qkgd, P)

        # computing attributes
        p = self.scm_props(M, W)

        #
        if verbose:
            print ("Q(f)_p =\n", Qf)
            print ("Q(k(g))_p+d =\n", Qkgd)
            print ("P(i,j), i = Q(f)_p, j = Q(k(g))_p+d =\n", P)
            print ("W =\n", W)
            print ("M = #{(i,j)|P(i,j), i=Q(f)_p, j=Q(k(g))_p+d} =\n", M)
            print ("p =\n", p)
        return M, p
        ########################################################################
        #SCM Map

class SCMio(object):
    
    ### SCM PROCEDURES
    def scm_feature_extraction_from_image_file(self, path, param):
        X = []
        fn = []
        for dirname, dirnames, filenames in os.walk(path):  
            for filename in filenames:
                if filename[0] != '.':	
                    fullfilename = os.path.join(dirname, filename)
                    f = io.imread(fullfilename)
                    f = color.rgb2gray(f)
                    f = (f - f.min())/(f.max() - f.min())
                    f = f.astype(float)
                    M, p = self.scm(f, f, param)
                    X.append(p)
                    fn.append(filename)
        return X, fn
    
    #
    def scm_feature_extraction_from_image_list(self, list, param):
        X = []
        for f in list:
            f = color.rgb2gray(f)
            f = (f - f.min())/(f.max() - f.min())
            f = f.astype(float)
            M, p = self.scm(f, f, param)
            X.append(p)
        return X

    #
    def scm_save_data_to_csv_file(self, filename, columndelimiter, data):
        with open(filename, 'wb') as csvfile:
            spamwriter = csv.writer(csvfile, delimiter = columndelimiter, quotechar = '|', quoting = csv.QUOTE_MINIMAL)
            for x in data:
                spamwriter.writerow(x)
        return 0

    #
    def scm_load_data_from_csv_file(self, filename, columndelimiter):
        d = []
        with open(filename, 'rb') as csvfile:
            spamreader = csv.reader(csvfile, delimiter = columndelimiter, quotechar = '|')
            for line in spamreader:
                d.append(line)
        return np.array(d)