Initial commit

Signed-off-by: jeonghwa yoo <[email protected]>

Author: jeongHwarr

MFCM.py

@@ -0,0 +1,261 @@
+#from skimage import data
+#from scipy import signal
+import sys
+import os
+from os import listdir
+from os.path import isfile, join
+import cv2
+import numpy as np 
+import matplotlib.pyplot as plt
+import time
+
+
+def make_spatial_distance_window(windowSize, center_coordi):    
+    center_y = center_coordi[1]
+    center_x = center_coordi[0]
+    
+    distance_window = np.zeros(windowSize)
+    for y in range(windowSize[1]):
+        for x in range(windowSize[0]):           
+            distance_window[y][x] = np.sqrt((y-center_y)**2+(x-center_x)**2) # Euclidean distance
+    return distance_window
+
+def find_reliable_set(window, windowSize):
+    median = np.median(window)    
+    median_window = np.full((windowSize[1],windowSize[0]), median) # shape:(windowSize[1],windowSize[0]), value=median
+        
+    difference_nei_med = window-median_window # neighbour - med
+        
+    deviation = np.sqrt((np.mean(pow(difference_nei_med,2))))
+            
+    reliable_mat = (difference_nei_med<=deviation).astype(int) # 1: reliable, 0: unreliable
+    return reliable_mat
+         
+ 
+def sliding_window(image, neighbour_effect, stepSize, windowSize):
+    '''slide a window across the image'''     
+     
+    center_coordi = (int(windowSize[1]/2),int(windowSize[0]/2)) # center coordination of the window
+    distance_window = make_spatial_distance_window(windowSize, center_coordi) # Euclidean distance matrix between center coordination and neighbor cordination
+    
+    filtered_image = []
+    
+    start = time.time()
+    for y in range(0, image.shape[0]-int(windowSize[0]/2), stepSize):
+        for x in range(0, image.shape[1]-int(windowSize[1]/2), stepSize):            
+            cur_window = image[y:y + windowSize[1], x:x + windowSize[0]] # neighbour window 
+            
+            if cur_window.shape[0]<windowSize[0] or cur_window.shape[1]<windowSize[1]:
+                continue
+                
+            if np.count_nonzero(cur_window)==0: # Pass if all values in window are zero
+                filtered_image.append(0)
+                continue
+            
+            #--------------Find reliable set matrix--------------        
+            reliable_mat = find_reliable_set(cur_window, windowSize)
+      
+            #--------------Weighting coefficients about window--------------           
+            center_window = np.full((windowSize[1],windowSize[0]), cur_window[center_coordi[1],center_coordi[0]])
+            difference_nei_cent = cur_window - center_window
+            
+            gray_deviation = np.sqrt(np.sum(pow(difference_nei_cent*reliable_mat,2))/np.count_nonzero(reliable_mat))             
+            if gray_deviation==0: gray_deviation = sys.float_info.epsilon
+            
+            coeff_gray = np.exp(-(pow(difference_nei_cent,2)*reliable_mat)/(neighbour_effect*gray_deviation))
+            coeff_distance = np.exp(-distance_window*reliable_mat)
+            
+            coeff_window = coeff_gray*coeff_distance
+                   
+            #--------------New intensity value of the center point in the window--------------
+            new_gray = np.sum(coeff_window*cur_window)/np.sum(coeff_window)               
+            filtered_image.append(round(new_gray))
+            
+    filtered_image = np.array(filtered_image)
+    print("Time :", time.time() - start)
+    return filtered_image
+    
+
+class MFCM():
+    def __init__(self, image, n_clusters, m, neighbour_effect, epsilon, max_iter, kernel_size):
+        '''Modified Fuzzy C-means clustering
+
+        <image>: 2D array, grey scale image.
+        <n_clusters>: int, number of clusters/segments to create.
+        <m>: float > 1, fuzziness parameter. A large <m> results in smaller
+             membership values and fuzzier clusters. Commonly set to 2.
+        <kernel_size>: int >= 1, size of neighborhood.
+        <neighbour_effect>: float, parameter which controls the influence extent of neighbouring pixels.
+        <max_iter>: int, max number of iterations.
+        '''
+
+        #-------------------Check inputs-------------------
+        if np.ndim(image) != 2:
+            raise Exception("<image> needs to be 2D (gray scale image).")
+        if n_clusters <= 0 or n_clusters != int(n_clusters):
+            raise Exception("<n_clusters> needs to be positive integer.")
+        if m < 1:
+            raise Exception("<m> needs to be >= 1.")
+        if kernel_size <=0 or kernel_size != int(kernel_size):
+            raise Exception("<kernel_size> needs to be positive integer.")
+        if epsilon <= 0:
+            raise Exception("<epsilon> needs to be > 0")
+
+        self.image = image
+        self.n_clusters = n_clusters
+        self.m = m
+        self.neighbour_effect = neighbour_effect
+        self.epsilon = epsilon
+        self.max_iter = max_iter
+        self.kernel_size = kernel_size
+
+        self.shape = image.shape # image shape
+        self.X = image.flatten().astype('float') # flatted image shape: (number of pixels,1) 
+        self.numPixels = image.size
+        
+    def initial_U(self):
+        U=np.zeros((self.num_gray, self.n_clusters))
+        idx = np.arange(self.num_gray)
+        for ii in range(self.n_clusters):
+            idxii = idx%self.n_clusters==ii
+            U[idxii,ii] = 1        
+        return U
+    
+    def update_U(self):
+        '''Compute weights'''
+        idx = np.arange(self.num_gray)
+        c_mesh,idx_mesh = np.meshgrid(self.C,idx)
+        power = -2./(self.m-1)
+        numerator = abs(idx_mesh-c_mesh)**power
+        denominator = np.sum(abs(idx_mesh-c_mesh)**power,axis=1)
+        return numerator/denominator[:,None]
+
+    def update_C(self):
+        '''Compute centroid of clusters'''
+        idx = np.arange(self.num_gray)
+        idx_reshape = idx.reshape(len(idx),1)
+        numerator = np.sum(self.histogram*idx_reshape*pow(self.U,self.m),axis=0)
+        denominator = np.sum(self.histogram*pow(self.U,self.m),axis=0)
+        return numerator/denominator
+   
+    def get_filtered_image(self):
+         # Create padding image
+        print("Getting filtered image..(This process can be time consuming.)") 
+        pad_size_y = int(self.kernel_size/2)
+        pad_size_x = int(self.kernel_size/2)   
+        pad_img = cv2.copyMakeBorder(self.image, pad_size_y, pad_size_y, pad_size_x, pad_size_x, cv2.BORDER_CONSTANT, value=0 ) # zero padding
+             
+        filtered_image = sliding_window(pad_img, self.neighbour_effect, stepSize=1, windowSize=(self.kernel_size,self.kernel_size))
+        dtype = self.image.dtype
+        self.filtered_image = filtered_image.reshape(self.shape).astype(dtype)
+    
+    def calculate_histogram(self):        
+        hist = cv2.calcHist([self.filtered_image],[0],None,[256],[0,256])
+        self.num_gray = len(hist)
+        self.histogram = hist
+                       
+    def form_clusters(self):
+        self.get_filtered_image() 
+        self.calculate_histogram()
+        
+        '''Iterative training'''        
+        d = 100
+        self.U = self.initial_U()
+        if self.max_iter != -1:
+            i = 0
+            while True:             
+                self.C = self.update_C()
+                old_u = np.copy(self.U)
+                self.U = self.update_U()
+                d = np.sum(abs(self.U - old_u))
+                print("Iteration %d : cost = %f" %(i, d))
+
+                if d < self.epsilon or i > self.max_iter:
+                    break
+                i+=1
+        else:
+            i = 0
+            while d > self.epsilon:
+                self.C = self.update_C()
+                old_u = np.copy(self.U)
+                self.U = self.update_U()
+                d = np.sum(abs(self.U - old_u))
+                print("Iteration %d : cost = %f" %(i, d))
+
+                if d < self.epsilon or i > self.max_iter:
+                    break
+                i+=1
+        self.segmentImage()
+
+
+    def deFuzzify(self):
+        return np.argmax(self.U, axis = 1)
+
+    def segmentImage(self):
+        '''Segment image based on max weights'''
+
+        result = self.deFuzzify()
+        
+        self.result = np.array(self.filtered_image, copy=True)
+        for i in range(len(result)):
+            self.result[self.result==i]=result[i]
+            
+        self.result = self.result.reshape(self.shape).astype('int')
+
+        return self.result
+    
+    
+def main(DIRECTORY, args):
+    IMG_PATH = DIRECTORY['IMG_PATH']
+    OUTPUT_FILT_IMG_PATH = DIRECTORY['OUTPUT_FILT_IMG_PATH']
+    OUTPUT_PLOT_PATH = DIRECTORY['OUTPUT_PLOT_PATH']
+    
+    IS_PLOT = args.plot_show 
+    IS_SAVE = args.plot_save
+    
+    files = [f for f in listdir(IMG_PATH) if isfile(join(IMG_PATH, f))] # read all files in IMG_PATH
+    
+    for file in files:
+        target_img_path = os.path.join(IMG_PATH,file)
+        try:
+            #--------------Lord image file--------------  
+            img= cv2.imread(target_img_path, cv2.IMREAD_GRAYSCALE) # cf. 8bit image-> 0~255
+
+            #--------------Clustering--------------  
+            cluster = MFCM(img, n_clusters=args.num_cluster, m=args.fuzziness, neighbour_effect=args.neighbour_effect, epsilon=args.epsilon, max_iter=args.max_iteration, kernel_size=args.win_size)
+            cluster.form_clusters()
+            result=cluster.result
+                     
+            #-------------------Plot and save result------------------------
+            if IS_PLOT:      
+                
+                fig=plt.figure(figsize=(12,8),dpi=100)
+            
+                ax1=fig.add_subplot(1,2,1)
+                ax1.imshow(img,cmap='gray')
+                ax1.set_title('image')
+            
+                ax2=fig.add_subplot(1,2,2)
+                ax2.imshow(result)
+                ax2.set_title('segmentation')
+                
+                plt.show(block=False)
+                plt.close()
+                
+            if IS_SAVE:
+                            
+                seg_result_path = os.path.join(OUTPUT_PLOT_PATH,"%s.png"%(os.path.splitext(file)[0]))
+                plt.imshow(result)
+                plt.savefig(seg_result_path, dpi=300)
+                plt.close()
+                
+                filtered_img_path = os.path.join(OUTPUT_FILT_IMG_PATH,"%s.png"%(os.path.splitext(file)[0]))
+                plt.imshow(cluster.filtered_image,cmap='gray')
+                plt.savefig(filtered_img_path, dpi=300)
+                plt.close()
+            
+        except IOError:
+            print("Error")
+
+if __name__ == '__main__':
+    main()

README.md

@@ -1,2 +1,79 @@
-# MFCM_segmentation
-Tissue Segmentation Using Modified Fuzzy C-Means Algorithm on Mammography (Image segmentation)
+
+# MFCM_segmentation (Python 3.5)
+### Tissue Segmentation Using Modified Fuzzy C-Means Algorithm on Mammography (Image segmentation)
+
+This code uses modified fuzzy c-means algorithm (MFCM) to do tissue segmentation on mammography. This code is a implementation of The paper:
+
+[1] Z. Chen and R. Zwiggelaar, 2010. 'A Modified Fuzzy C Means Algorithm for Breast Tissue Density Segmentation in Mammograms’, IEEE/Information Technology and Applications in Biomedicine (ITAB) (https://ieeexplore.ieee.org/document/5687751)
+
+## How to run?
+
+### Run on mini data as default option 
+You just run <u>main.py<\u> in editor or enter <u>python main.py<\u> in command prompt.
+(In <i>img<\i> directory, there are 2 images for the test. These images are part of mini-MIAS database. (http://peipa.essex.ac.uk/info/mias.html) ) 
+
+### User mode
+If you want to replace the mini data with your own data, put your images to img directory or edit path for your direrectory in <u>main.py<\u>.   
+If you change parameters of your experiment, you can change parameters by changing the default value of the argument in <u>main.py<\u> or you can you the command line in command prompt.
+- You can see all the adjustable parameters and usage. </pre> python main.py --help</pre>
+(For example, if you want to cluster with 5 clusters, enter this command in your command prompt.  </u>python main.py -c 5<\u>)
+
+## Results
+
+<pre>
+Getting filtered image..(This process can be time consuming.)
+Time : 53.16970419883728
+Iteration 0 : cost = 383.503160
+Iteration 1 : cost = 130.096540
+Iteration 2 : cost = 52.057745
+Iteration 3 : cost = 36.463114
+Iteration 4 : cost = 33.625992
+Iteration 5 : cost = 40.136341
+Iteration 6 : cost = 42.462534
+Iteration 7 : cost = 29.869172
+Iteration 8 : cost = 16.829851
+Iteration 9 : cost = 10.073875
+Iteration 10 : cost = 6.756716
+Iteration 11 : cost = 4.957472
+Iteration 12 : cost = 3.912657
+Iteration 13 : cost = 3.158448
+Iteration 14 : cost = 2.573333
+Iteration 15 : cost = 2.098461
+Iteration 16 : cost = 1.710719
+Iteration 17 : cost = 1.394245
+Iteration 18 : cost = 1.136231
+Iteration 19 : cost = 0.925994
+Iteration 20 : cost = 0.754842
+Iteration 21 : cost = 0.615388
+Iteration 22 : cost = 0.501762
+Iteration 23 : cost = 0.409169
+Iteration 24 : cost = 0.333715
+Iteration 25 : cost = 0.272221
+Iteration 26 : cost = 0.222079
+Iteration 27 : cost = 0.181186
+Iteration 28 : cost = 0.147834
+Iteration 29 : cost = 0.120628
+Iteration 30 : cost = 0.098434
+Iteration 31 : cost = 0.080326
+Iteration 32 : cost = 0.065552
+Iteration 33 : cost = 0.053496
+Iteration 34 : cost = 0.043659
+</pre>
+
+|Image| Result |
+|:---:|:---: |
+|mdb321.jpg|<img src = 'result.png'>|
+|mdb322.jpg|<img src = 'result2.png'>|  
+
+
+## Troubleshooting:
+1. If it is not an 8-bit image, the code needs to be modified.
+2. When there is a problem with the environment, you can try this command line in your command prompt. <pre> pip install -r requirements.txt </pre>
+
+## References
+* [1] [Z. Chen and R. Zwiggelaar "A Modified Fuzzy C Means Algorithm for Breast Tissue Density Segmentation in Mammograms." IEEE/Information Technology and Applications in Biomedicine (ITAB) 2010.](https://ieeexplore.ieee.org/document/5687751)
+* [2] [J. Song and Z. Zhang "A Modified Robust FCM Model with Spatial Constraints for Brain MR Image Segmentation." Information 2019.](https://www.researchgate.net/publication/331278874_A_Modified_Robust_FCM_Model_with_Spatial_Constraints_for_Brain_MR_Image_Segmentation)
+* https://github.com/ab93/SIFCM
+
+## TODO:
+- [ ] It takes a lot of time to get the filtered image. This problem can be solved by using numpy.lib.stride_tricks.as_strided.  

main.py

@@ -0,0 +1,49 @@
+import os
+import argparse
+
+import MFCM
+
+WORKSPACE = os.path.split(os.getcwd())[-1] # current workspace 
+IMG_PATH = os.path.join(os.getcwd(),"img") # path for input image directory
+OUTPUT_PATH = os.path.join(os.getcwd(),"output") # path for output directory
+OUTPUT_PLOT_PATH = os.path.join(OUTPUT_PATH,'segmentation') # path for output (plot directory)
+OUTPUT_FILT_IMG_PATH = os.path.join(OUTPUT_PATH,'filtered_img') # path for output (filtered image directory)
+    
+def get_args():
+    parser = argparse.ArgumentParser(description="Tissue Segmentation Using Modified Fuzzy C-Means Algorithm on Mammography.")
+   
+    #-----------------Fundamental parameter-----------------
+    parser.add_argument('-c', '--num_cluster', default='4', type=int,
+                         help="Number of cluster")
+    parser.add_argument('-m', '--fuzziness', default='2', type=int,
+                         help="fuzziness degree")
+    parser.add_argument('-i', '--max_iteration', default='100', type=int,
+                         help="max number of iterations.")
+    parser.add_argument('-e', '--epsilon', default='0.05', type=float,
+                         help="threshold to check convergence.")
+    #-----------------User parameter-----------------
+    parser.add_argument('--plot_show', default=1, choices=[0,1],
+                         help="Show plot about result")
+    parser.add_argument('--plot_save', default=1, choices=[0,1],
+                         help="Save plot about result")
+    #-----------------MFCM parameter-----------------
+    parser.add_argument('-w', '--win_size', default='5', type=int,
+                         help="Window size of MFCM algorithm")
+    parser.add_argument('-n', '--neighbour_effect', default='3', type=float,
+                         help="Effect factor of the graylevel which controls the influence extent of neighbouring pixels.")
+    
+    args = parser.parse_args() 
+    return args
+
+if __name__ == '__main__':
+  
+    args = get_args()
+    
+    DIRECTORY = {}   
+    DIRECTORY['WORKSPACE'] = WORKSPACE
+    DIRECTORY['IMG_PATH'] = IMG_PATH
+    DIRECTORY['OUTPUT_PATH'] = OUTPUT_PATH
+    DIRECTORY['OUTPUT_PLOT_PATH'] = OUTPUT_PLOT_PATH
+    DIRECTORY['OUTPUT_FILT_IMG_PATH'] = OUTPUT_FILT_IMG_PATH
+    
+    MFCM.main(DIRECTORY, args)