diff --git a/code/lung_feature_patch_allPatients.py b/code/lung_feature_patch_allPatients.py
index f138775b265b2fe0ae44672d73d305039ba4d3a9..8b8e9d42e634d5b23c110bbcc33ef2abbcfc9158 100755
--- a/code/lung_feature_patch_allPatients.py
+++ b/code/lung_feature_patch_allPatients.py
@@ -16,20 +16,20 @@ Extended functionality by Monica J. Emerson.
 - Study of several diseases. 
 - Added functionality for supporting the analysis of different data versions.
 - Computation of health scores per image, sample and patient.
-- Visualisation of cluster centres as a grid.
+- Visualisation of cluster centers as a grid.
 - Boxplots to compare probabilities across samples and to clinical values.
-- Normalisation of intensities across images and channels.
+- Normalisation of int_sum_list across images and channels.
 - Possibility to ignore the background (air phase).
 - Visualisation of assignment images to inspect results and support the development of the approach to ignore background.
 - Implementation and investigation of feature variations (colour, bnw, bnw+colour).
 - Study of the parameters (nr clusters and scale - relative patch/image size) 
-- Extended visualisation of cluster centres to support the comparison across diseases and parameters.
-    a) Visualisation of cluster centres split into channels.
+- Extended visualisation of cluster centers to support the comparison across diseases and parameters.
+    a) Visualisation of cluster centers split into channels.
     b) Compute a population (p) and condition probability (c) value for each cluster.
     c) Identify the presence of "weak" clusters. If they exist, rerun kmeans.
     d) Select and visualise characteristic clusters for the conditions based on p and c.
     e) Plot all clusters in the population/condition probability space.
-    f) Order cluster centres according to condition probability 
+    f) Order cluster centers according to condition probability 
 """
 
 import numpy as np
@@ -43,6 +43,7 @@ import os
 from datetime import datetime
 import sys
 from matplotlib.offsetbox import OffsetImage, AnnotationBbox
+from math import ceil
 
 startTime = datetime.now()
 
@@ -51,14 +52,15 @@ plt.close('all')
 sc_fac = 0.25 #25 #25 #0.5 # scaling factor of the image
 patch_size = 17
 nr_clusters = 100 # number of clusters
-
+colour_mode = 'colour' #'bnw' #colour
+    
 #%% Directories
 
 version = 'corrected_bis'
-preprocessing = ''  #'/preprocessed_ignback/' #''(none)
-colour_mode = 'colour' #'bnw' #colour
-disease = ['sarcoidosis']#diseased (mix, 2 of each condition)' #''emphysema' 'sarcoidosis'
-
+preprocessing = '/preprocessed_ignback/' #''(none)
+disease = 'diseased' #diseased (mix, 2 of each condition)' #''emphysema' 'sarcoidosis'
+conditions = [disease]
+conditions.insert(0,'control')
 
 # input directories - images start with the name 'frame'
 dir_in = '../maxProjImages_'+version + preprocessing 
@@ -66,16 +68,15 @@ dir_in = '../maxProjImages_'+version + preprocessing
 # dir_sick = dir_in + disease + '/' #191216_100a/'
 base_name = 'frame'
 
-
 #output directories
 dir_results = '../results_monj/patches/data_' + version+'/' #rerun just to make sure not to muck up the results
 os.makedirs(dir_results, exist_ok = True)  
 
 dir_probs = dir_results + disease + '_' + colour_mode + '_%dclusters_%ddownscale_%dpatchsize/'%(nr_clusters,1/sc_fac,patch_size)
-ma.make_output_dirs(dir_probs, disease)
+ma.make_output_dirs(dir_probs, conditions)
 
-dir_probs_withBack = dir_probs + 'withBackground'
-ma.make_output_dirs(dir_probs_withBack, disease)
+dir_probs_noBack = dir_probs + 'noBackground/'
+ma.make_output_dirs(dir_probs_noBack, conditions)
 
 # if not os.path.exists(dir_probs):
 #         os.mkdir(dir_probs)
@@ -88,41 +89,42 @@ ma.make_output_dirs(dir_probs_withBack, disease)
 #Read maximum projection images (runtime ~= 20 s )
 print('Reading maximum projection images')
 
-max_img_list_control = ma.read_max_imgs(dir_in+ 'control', base_name)
-max_img_list_sick = ma.read_max_imgs(dir_in+ disease, base_name)
+max_img_list_control = ma.read_max_imgs(dir_in + 'control', base_name)
+max_img_list_sick = ma.read_max_imgs(dir_in + disease, base_name)
 
- 
-# max_im_list_control = ma.read_max_ims(dir_in_max + 'control/' ,base_name,sc_fac,colour_mode)
-
-# max_im_list_sick = []    
-# for sample in dir_list_sick:
-#     in_dir = dir_sick + sample + '/'
-#     dir_list = [dI for dI in sorted(os.listdir(in_dir)) if dI[0:len(base_name)]==base_name]
-#     frames_list = []
-#     for ind, frame in enumerate(dir_list):
-#         frame_path = in_dir + frame
-#         if colour_mode == 'bnw':
-#             img = skimage.color.rgb2gray(skimage.io.imread(frame_path).astype(float))
-#             max_im_list_sick += [skimage.transform.rescale(img, sc_fac)] 
-#         else:
-#             img = skimage.io.imread(frame_path).astype(float)
-#             max_im_list_sick += [skimage.transform.rescale(img, sc_fac, multichannel=True)] 
-
-#TO DO: Rescale max proj. images, overwrite original variables
-#TO DO: Compute bnw version, but keep the colour one for displaying it at the end
+#TO DO: Consider removing colour mode and rescaling from the read function.
+
+#%% Prepare maximum projection images for feature extraction
+
+#Flatten and merge conditions
+max_img_list_control_flat = ma.flatten_list(max_img_list_control) 
+max_img_list_sick_flat = ma.flatten_list(max_img_list_sick)
+
+#Rescale
+max_img_list_control_processed = []
+for max_img in ma.flatten_list(max_img_list_control) :
+    # if colour_mode == 'bnw':
+    #     max_img_list_flat_processed += [skimage.transform.rescale(skimage.color.rgb2gray(max_img), sc_fac, preserve_range = True).astype('uint8')]
+    # else:
+    max_img_list_control_processed += [skimage.transform.rescale(max_img, sc_fac, preserve_range = True, multichannel=True).astype('uint8')]
+
+max_img_list_sick_processed = []
+for max_img in ma.flatten_list(max_img_list_sick) :
+    max_img_list_sick_processed += [skimage.transform.rescale(max_img, sc_fac, preserve_range = True, multichannel=True).astype('uint8')]
 
 #%% Compute patches
 patch_feat_list_control = []
-for max_im in max_im_list_control:
-    patch_feat_list_control += [ma.ndim2col_pad(max_im, (patch_size, patch_size),norm=False).transpose()]
+for max_img in max_img_list_control_processed:
+    patch_feat_list_control += [ma.ndim2col_pad(max_img, (patch_size, patch_size),norm=False).transpose()]
 
 patch_feat_list_sick = []
-for max_im in max_im_list_sick:
-    patch_feat_list_sick += [ma.ndim2col_pad(max_im, (patch_size, patch_size),norm=False).transpose()]
+for max_img in max_img_list_sick_processed:
+    patch_feat_list_sick += [ma.ndim2col_pad(max_img, (patch_size, patch_size),norm=False).transpose()]
+
+patch_feat_total = patch_feat_list_control + patch_feat_list_sick
 
-patch_feat_total = []
-patch_feat_total += patch_feat_list_control
-patch_feat_total += patch_feat_list_sick
+#TO DO: Consider turning patches to bnw here instead
+#so we can read the protein content for each patch as the sum of the patch int_sum_list
 
 #%% features for clustering
 nr_keep = 10000 # number of features randomly picked for clustering 
@@ -137,101 +139,155 @@ for patch_feat in patch_feat_total:
 
 #%% k-means clustering
 batch_size = 1000
-th_nr_pathesINcluster = 5
+th_nr_pathesINcluster = 10
 
-if os.path.exists(dir_probs + 'array_cluster_centres'+colour_mode+'.npy'):
-     cluster_centres = np.load(dir_probs + 'array_cluster_centres'+colour_mode+'.npy')
-     #kmeans = sklearn.cluster.MiniBatchKMeans(n_clusters=nr_clusters, init = cluster_centres, batch_size = batch_size)
+#If cluster centers were already computed, use those
+if os.path.exists(dir_probs + 'array_cluster_centers_'+colour_mode+'.npy'):
+     cluster_centers = np.load(dir_probs + 'array_cluster_centers_'+colour_mode+'.npy')
+     #kmeans = sklearn.cluster.MiniBatchKMeans(n_clusters=nr_clusters, init = cluster_centers, batch_size = batch_size)
      #kmeans.fit(patch_feat_to_cluster)
      kmeans = sklearn.cluster.MiniBatchKMeans(n_clusters=nr_clusters, batch_size = batch_size)
-     kmeans.cluster_centers_=cluster_centres
-     reusing_clusters = True
-else:
+     kmeans.cluster_centers_ = cluster_centers
+
+#Compute cluster centers and save if there is no weak clusters
+else: 
     kmeans = sklearn.cluster.MiniBatchKMeans(n_clusters=nr_clusters, batch_size = batch_size)
     kmeans.fit(patch_feat_to_cluster)
-    all_cluster_centres = kmeans.cluster_centers_
     
-    #Cluster statistics
+    #Nr. patches that have contributed to each cluster
     features_in_cluster = []
     for cluster in range(0,nr_clusters):
         features_in_cluster += [[ind for ind,i in enumerate(kmeans.labels_) if i==cluster]]  
     nr_feat_in_cluster =  [len(i) for i in features_in_cluster]  
     
+    #Weak clusters are those composed by very few patches
     nr_weakClusters = len([1 for i in nr_feat_in_cluster if i<th_nr_pathesINcluster])
     
+    #If there are weak clusters, the clustering should be recomputed
     if nr_weakClusters!=0:
         sys.exit(str(nr_weakClusters)+ " clusters composed of less than "+str(th_nr_pathesINcluster)+" images") 
     else:
-        np.save(dir_probs + 'array_cluster_centres'+colour_mode+'.npy', all_cluster_centres)    # .npy extension is added if not given
+        np.save(dir_probs + 'array_cluster_centers_'+colour_mode+'.npy', kmeans.cluster_centers_)    # .npy extension is added if not given
     
-     
-#%% Read background pixels
-dir_background = dir_in + 'background/'
-dir_background_list = [dI for dI in os.listdir(dir_background) if os.path.isdir(os.path.join(dir_background,dI))]
-
-fig, axs = plt.subplots(2,len(dir_background_list), sharex=True, sharey=True)
-patch_feat_back = []
-for ind, directory in enumerate(dir_background_list):
-    #load images and corresponding background labels
-    file_names = [f for f in os.listdir(dir_background +directory) if f.endswith('.png')]
-    im_file = [f for f in file_names if not f.startswith('back')].pop()
-    label_file = [f for f in file_names if f.startswith('back')].pop() 
-    im_back = skimage.io.imread(dir_background + directory + '/' + im_file).astype('uint8')
-    label_back = skimage.color.rgb2gray(skimage.io.imread(dir_background + directory + '/' + label_file).astype('float'))
-    label_back += -np.min(label_back)
-    label_back = label_back.astype('bool')
-    #plot imagesand corresonding labels
-    axs[0][ind].imshow(im_back)
-    axs[1][ind].imshow(label_back,'gray')
-    plt.show()
-    #compute features
-    if colour_mode == 'bnw':
-        im_back = skimage.transform.rescale(skimage.color.rgb2gray(im_back.astype(float)), sc_fac, multichannel=False)
-    else :
-        im_back = skimage.transform.rescale(im_back.astype(float), sc_fac, multichannel=True)
-    im_feat = ma.ndim2col_pad(im_back, (patch_size, patch_size)).transpose()
-    patch_feat_back += [im_feat[(skimage.transform.rescale(label_back, sc_fac, multichannel=False)==True).ravel(),:]]
-
-#%% Plot all cluster centres
-
-plot_grid_cluster_centres(kmeans.cluster_centers_)
-
-if nr_clusters==100:
-    fig, axs = plt.subplots(10,10, figsize=(5,5), sharex=True, sharey=True)
-if nr_clusters==200:
-    w, h = plt.figaspect(2.)
-    fig, axs = plt.subplots(20,10, figsize=(w,h), sharex=True, sharey=True)  
-if nr_clusters==1000:
-    fig, axs = plt.subplots(100,100, figsize=(5,5), sharex=True, sharey=True)  
-
-intensities = []
-for ax, cluster_nr in zip(axs.ravel(), np.arange(0,nr_clusters)):
-    if colour_mode == 'bnw':
-        cluster_centre = np.reshape(kmeans.cluster_centers_[cluster_nr,:],(patch_size,patch_size))
-        intensities += [sum((cluster_centre).ravel())] 
-        ax.imshow(cluster_centre.astype('uint8'),cmap='gray')
+#%% Plot all cluster centers (grid view)
+
+#grid dimensions
+size_x = round(nr_clusters**(1/2))
+size_y = ceil(nr_clusters/size_x)
+
+#figure format
+w, h = plt.figaspect(size_x/size_y)
+fig, axs = plt.subplots(size_x,size_y, figsize=(w,h), sharex=True, sharey=True)
+
+print('Grid size: ', size_x, size_y, 'Figure size: ', w, h)
+
+#plot cluster centers
+int_sum_list = []
+ax_list = axs.ravel()
+for ind in np.arange(0,nr_clusters):
+    if colour_mode == 'bnw': #in bnw + colour give the clusters a uniform colour
+        cluster_centre = np.reshape(kmeans.cluster_centers_[ind,:],(patch_size,patch_size))
+        int_sum_list += [sum((cluster_centre).ravel())] 
+        ax_list[ind].imshow(cluster_centre.astype('uint8'),cmap='gray')
     else:
-        cluster_centre = np.transpose((np.reshape(kmeans.cluster_centers_[cluster_nr,:],(3,patch_size,patch_size))),(1,2,0))
-        intensities += [sum((np.max(cluster_centre,2)).ravel())] 
-        ax.imshow(cluster_centre.astype('uint8'))
+        cluster_centre = np.transpose((np.reshape(kmeans.cluster_centers_[ind,:],(3,patch_size,patch_size))),(1,2,0))
+        int_sum_list += [sum((np.max(cluster_centre,2)).ravel())] 
+        ax_list[ind].imshow(cluster_centre.astype('uint8'))
     
 plt.setp(axs, xticks=[], yticks=[])
-plt.savefig(dir_probs + 'clusterCentres_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+plt.savefig(dir_probs + 'clustercenters_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
  
-
 #%% Histograms for background, healthy and sick
-#hist_control, assignment_list_control, hist_back, assignment_back = ma.compute_assignment_hist(patch_feat_list_control, kmeans, background_feat=im_feat_back)
-
-hist_background, assignment_list_background = ma.compute_assignment_hist(patch_feat_back, kmeans)
 hist_control, assignment_list_control= ma.compute_assignment_hist(patch_feat_list_control, kmeans)
 hist_sick, assignment_list_sick = ma.compute_assignment_hist(patch_feat_list_sick, kmeans)
 
-#%% Cluster centres in the 2d space determined by the relationshop between histogram
-occurrence_ratio = hist_control/hist_sick
-occurrence_ratio[occurrence_ratio<1] = -1/occurrence_ratio[occurrence_ratio<1] 
+#%% show bar plot of healthy and sick
+
+fig, ax = plt.subplots(1,1)
+ax.bar(np.array(range(0,nr_clusters))-0.25, hist_control, width = 0.5, label='Control', color = 'r')
+ax.bar(np.array(range(0,nr_clusters))+0.25, hist_sick, width = 0.5, label='Sick', color = 'b')
+ax.legend()
+plt.savefig(dir_probs + 'assignmentHistograms_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+
+#%% Compute populated, occurrence and importance cluster measures from the histogram
+
+#OCCURRENCE RATIO: How more often does it occur in one of the conditions
+occurrence_ratio = np.empty((nr_clusters))
+ind_control = hist_control>hist_sick
+ind_sick = ~ind_control
+control_cond = (hist_sick!=0) & ind_control
+occurrence_ratio[control_cond] = hist_control[control_cond]/hist_sick[control_cond]
+sick_cond = (hist_control!=0) & ind_sick
+occurrence_ratio[sick_cond] = -hist_sick[sick_cond]/hist_control[sick_cond]
+
+#POPULATED How trustworthy is the cluster? 
 populated = hist_control+hist_sick
 
-plt.hist2d(occurrence_ratio,populated)
+weights = populated/populated.max()
+
+#Cluster IMPORTANCE: occurrence weighed by POPULATED
+importance = ((2*np.exp(populated/populated.max())-1)*occurrence_ratio).reshape((size_x,size_y))
+importance[importance>2*importance.std()] = 2*importance.std()
+importance[importance<-2*importance.std()] = -2*importance.std()
+
+
+#%% Plot measures as images
+
+#OCCURRENCE RATIO
+fig, ax = plt.subplots(figsize=(5,5))
+plt.imshow(occurrence_ratio.reshape((size_x,size_y)),cmap = 'PiYG')
+plt.colorbar()
+plt.setp(axs, xticks=[], yticks=[])
+plt.title('Condition dominance \n Dominant condition (+ control, - disease)')
+
+#POPULATED
+fig, ax = plt.subplots(figsize=(5,5))
+plt.imshow(populated.reshape((size_x,size_y)),cmap = 'PiYG')
+plt.colorbar()
+plt.setp(axs, xticks=[], yticks=[])
+plt.title('Cluster population ratio \n Distribution of patches across clusters')
+
+#Cluster IMPORTANCE: occurrence times POPULATED
+fig, ax = plt.subplots(figsize=(5,5))
+plt.imshow((populated*occurrence_ratio).reshape((size_x,size_y)),cmap = 'PiYG')
+plt.colorbar() 
+plt.setp(ax, xticks=[], yticks=[])
+plt.title('Cluster importance \n Occurrence ratio*population')
+
+#Saturated cluster IMPORTANCE: occurrence times POPULATED Saturated
+fig, ax = plt.subplots(figsize=(5,5))
+
+plt.imshow(importance, cmap = 'PiYG', vmin=-2*importance.std(), vmax=2*importance.std())
+colorbar = plt.colorbar() 
+limit_cmap = max(occurrence_ratio)
+ticks = np.linspace(-limit_cmap,limit_cmap,len(colorbar.get_ticks()))
+tick_labels = [str(i) for i in ticks.tolist()]
+colorbar.set_ticklabels(np.linspace(-limit_cmap,limit_cmap,len(colorbar.get_ticks())))
+plt.setp(ax, xticks=[], yticks=[])
+plt.title('Cluster importance for the condition \n Healthy (green), ' + str(disease) + ' (pink)')
+
+plt.savefig(dir_probs + 'clusterCentreImportance_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+ 
+#%% Plot grid of ordered clusters
+
+#Order the clusters according to importance
+clusters_control_importance_sorted = [np.arange(0,100)[i] for i in np.argsort(importance.flatten()) if (importance.flatten())[i]>0]
+clusters_control_importance_sorted.reverse()
+
+clusters_sick_importance_sorted = [np.arange(0,100)[i] for i in np.argsort(importance.flatten()) if (importance.flatten())[i]<0]
+
+clusters_importance_sorted = [np.arange(0,100)[i] for i in np.argsort(importance.flatten())]
+clusters_importance_sorted.reverse()
+
+occurrence_ratio_sorted = [occurrence_ratio[i] for i in np.argsort(importance.flatten()) ]
+occurrence_ratio_sorted.reverse()
+
+ma.plot_grid_cluster_centers(kmeans.cluster_centers_, clusters_importance_sorted, patch_size, colour_mode = 'colour', occurrence = occurrence_ratio_sorted)
+
+plt.suptitle('Cluster centers \n from more characteristic of control (0, ' + str(len(clusters_control_importance_sorted))  + ') to more characteristic of sick')
+plt.savefig(dir_probs + 'clustercenters_sorted_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+ 
+#%% Plot cluster centers in the 2d populated/occurrence space
 
 fig_control, ax_control = plt.subplots(figsize=(15,15))
 ax_control.scatter(populated[occurrence_ratio>1], occurrence_ratio[occurrence_ratio>1]) 
@@ -254,74 +310,82 @@ for x0, y0, cluster_nr in zip(populated, occurrence_ratio, np.arange(0,nr_cluste
         else:
             ab = AnnotationBbox(OffsetImage(cluster_centre.astype('uint8')), (x0, y0), frameon=False)
         ax_control.add_artist(ab)
-        plt.savefig(dir_probs + 'controlClusterCentres_2Dspace_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+        plt.savefig(dir_probs + 'controlClustercenters_2Dspace_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
     else:
         if colour_mode == 'bnw':
             ab = AnnotationBbox(OffsetImage(cluster_centre.astype('uint8'),cmap='gray'), (x0, -y0), frameon=False)
         else:
            ab = AnnotationBbox(OffsetImage(cluster_centre.astype('uint8')), (x0, -y0), frameon=False)
         ax_sick.add_artist(ab)
-        plt.savefig(dir_probs + 'sickClusterCentres_2Dspace_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+        plt.savefig(dir_probs + 'sickClustercenters_2Dspace_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
 
     
-    
-#%% show bar plot of healthy and sick
+#%% Find background clusters
+th_back = 10000
+#background clusters
+clusters_background_intBased = [i for i in range(len(int_sum_list)) if int_sum_list[i] < 8000]
+#clusters_background_annotBased = [ind for ind, value in enumerate(hist_background)  if value>0]
 
-fig, ax = plt.subplots(1,1)
-ax.bar(np.array(range(0,nr_clusters)), hist_background, width = 1)
-plt.show()
-plt.savefig(dir_probs + 'backgroundHistogram_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
-  
-fig, ax = plt.subplots(1,1)
-ax.bar(np.array(range(0,nr_clusters))-0.25, hist_control, width = 0.5, label='Control', color = 'r')
-ax.bar(np.array(range(0,nr_clusters))+0.25, hist_sick, width = 0.5, label='Sick', color = 'b')
-ax.legend()
-plt.savefig(dir_probs + 'assignmentHistograms_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
-  
+clusters_background = clusters_background_intBased
+#clusters_background = list(set(clusters_background_intBased) | set(clusters_background_annotBased)) 
+print('Background clusters'+str(clusters_background))
 
-#%% Find background and characteristic clusters
+#%% Plot grid of ordered, non-background clusters that are populated over a thresh
+th_populated = 0.4/nr_clusters
+th_occurr = 1.25
 
-#background clusters
-clusters_background_intBased = [i for i in range(len(intensities)) if intensities[i] < 7000]
-clusters_background_annotBased = [ind for ind, value in enumerate(hist_background)  if value>0]
+#control
+clusters_control_importance_sorted_noBackground = [i for i in clusters_control_importance_sorted if (i not in clusters_background)&(populated[i]>th_populated)&(occurrence_ratio[i]>th_occurr)]
 
-clusters_background = list(set(clusters_background_intBased) | set(clusters_background_annotBased)) 
-print('Background clusters'+str(clusters_background))
+ma.plot_grid_cluster_centers(kmeans.cluster_centers_, clusters_control_importance_sorted_noBackground, patch_size, colour_mode = 'colour')
+
+plt.suptitle('Characteristic clusters for control, excluding background, \n from most important to least')
+plt.savefig(dir_probs_noBack + 'clustercenters_control_sorted_noBack%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+ 
+#disease
+clusters_sick_importance_sorted_noBackground = [i for i in clusters_sick_importance_sorted if (i not in clusters_background)&(populated[i]>th_populated)&(abs(occurrence_ratio[i])>th_occurr)]
 
-#characteristic clusters
-th_proportion = 2#2 #2.4
-th_populated = 0.01#0.005#0.015
+ma.plot_grid_cluster_centers(kmeans.cluster_centers_, clusters_sick_importance_sorted_noBackground, patch_size, colour_mode = 'colour')
 
-clusters_sick = (hist_sick>th_populated)&(hist_sick>th_proportion*hist_control)
-clusters_sick = [ind for ind,value in enumerate(clusters_sick) if value == True]
-clusters_control = (hist_control>th_populated)&(hist_control>th_proportion*hist_sick)
-clusters_control = [ind for ind,value in enumerate(clusters_control) if value == True]
+plt.suptitle('Charactertistic clusters for ' + disease + ', excluding background, \n from most important to least')
+plt.savefig(dir_probs_noBack + 'clustercenters_sick_sorted_noBack%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+ 
+# #%% Find characteristic clusters
+# th_proportion = 2#2 #2.4
+# th_populated_condition = 0.01#0.005#0.015
 
-#eliminate backgorund clusters if contained here
-clusters_control = [i for i in clusters_control if i not in clusters_background]
-clusters_sick = [i for i in clusters_sick if i not in clusters_background]
+# clusters_sick = (hist_sick>th_populated_condition)&(hist_sick>th_proportion*hist_control)
+# clusters_sick = [ind for ind,value in enumerate(clusters_sick) if value == True]
+# clusters_control = (hist_control>th_populated_condition)&(hist_control>th_proportion*hist_sick)
+# clusters_control = [ind for ind,value in enumerate(clusters_control) if value == True]
 
-print('Clusters characteristic of the ' + disease + ' tissue',clusters_sick)
-print('Clusters characteristic of the control tissue',clusters_control)
+# #eliminate background clusters if contained here
+# clusters_control = [i for i in clusters_control if i not in clusters_background]
+# clusters_sick = [i for i in clusters_sick if i not in clusters_background]
 
-#%% Plot centres of the characteristic clusters 
-cluster_centres_control = kmeans.cluster_centers_[clusters_control]
+# print('Clusters characteristic of the ' + disease + ' tissue',clusters_sick)
+# print('Clusters characteristic of the control tissue',clusters_control)
+
+
+#%% Plot centers of the characteristic clusters 
+nr_clusters_display = 10
+clusters_control = clusters_control_importance_sorted_noBackground[1:10] 
 
 #control clusters and contrast enhanced
-cluster_centres_control = kmeans.cluster_centers_[clusters_control]
+cluster_centers_control = kmeans.cluster_centers_[clusters_control]
 fig, axs = plt.subplots(1,len(clusters_control), figsize=(len(clusters_control)*3,3), sharex=True, sharey=True)
-fig.suptitle('Cluster centres for control')
+fig.suptitle('Cluster centers for control')
 if colour_mode!='bnw':
     fig_split, axs_split = plt.subplots(3,len(clusters_control), figsize=(len(clusters_control)*3,9), sharex=True, sharey=True)
-    fig_split.suptitle('Control centres split channels (contrast enhanced)')
+    fig_split.suptitle('Control centers split channels (contrast enhanced)')
 for l in np.arange(0,len(clusters_control)):
     if colour_mode == 'bnw':
-        cluster_centre = np.reshape(cluster_centres_control[l,:],(patch_size,patch_size))
+        cluster_centre = np.reshape(cluster_centers_control[l,:],(patch_size,patch_size))
         axs[l].imshow(cluster_centre.astype(np.uint8),cmap='gray')
         axs[l].axis('off')
         axs[l].set_title(clusters_control[l])
     else:    
-        cluster_centre = np.transpose((np.reshape(cluster_centres_control[l,:],(3,patch_size,patch_size))),(1,2,0))
+        cluster_centre = np.transpose((np.reshape(cluster_centers_control[l,:],(3,patch_size,patch_size))),(1,2,0))
         axs_split[0][l].imshow(cluster_centre[...,0].astype(np.uint8),cmap='gray')
         axs_split[0][l].axis('off')
         axs_split[0][l].set_title(clusters_control[l])
@@ -330,28 +394,30 @@ for l in np.arange(0,len(clusters_control)):
         axs_split[2][l].imshow(cluster_centre[...,2].astype(np.uint8),cmap='gray')
         axs_split[2][l].axis('off')
         plt.figure(fig_split.number),
-        plt.savefig(dir_probs + 'clusterCentres_control_splitContrastEnhanced_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+        plt.savefig(dir_probs + 'clustercenters_control_splitContrastEnhanced_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
         axs[l].imshow(cluster_centre.astype(np.uint8))
         axs[l].axis('off')
         axs[l].set_title(clusters_control[l])
     plt.figure(fig.number),
-    plt.savefig(dir_probs + 'clusterCentres_control_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+    plt.savefig(dir_probs + 'clustercenters_control_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
     
 #sick clusters and contrast enhanced
-cluster_centres_sick = kmeans.cluster_centers_[clusters_sick]
+clusters_sick = clusters_sick_importance_sorted_noBackground[1:10] 
+
+cluster_centers_sick = kmeans.cluster_centers_[clusters_sick]
 fig, axs = plt.subplots(1,len(clusters_sick), figsize=(len(clusters_sick)*3,3), sharex=True, sharey=True)
-fig.suptitle('Cluster centres for sick')
+fig.suptitle('Cluster centers for sick')
 if colour_mode!='bnw':
     fig_split, axs_split = plt.subplots(3,len(clusters_sick), figsize=(len(clusters_sick)*3,9), sharex=True, sharey=True)
-    fig_split.suptitle('sick centres split channels (contrast enhanced)')
+    fig_split.suptitle('sick centers split channels (contrast enhanced)')
 for l in np.arange(0,len(clusters_sick)):
     if colour_mode == 'bnw':
-        cluster_centre = np.reshape(cluster_centres_sick[l,:],(patch_size,patch_size))
+        cluster_centre = np.reshape(cluster_centers_sick[l,:],(patch_size,patch_size))
         axs[l].imshow(cluster_centre.astype(np.uint8),cmap='gray')
         axs[l].axis('off')
         axs[l].set_title(clusters_sick[l])
     else:
-        cluster_centre = np.transpose((np.reshape(cluster_centres_sick[l,:],(3,patch_size,patch_size))),(1,2,0))
+        cluster_centre = np.transpose((np.reshape(cluster_centers_sick[l,:],(3,patch_size,patch_size))),(1,2,0))
         axs_split[0][l].imshow(cluster_centre[...,0].astype(np.uint8),cmap='gray')
         axs_split[0][l].axis('off')
         axs_split[0][l].set_title(clusters_sick[l])
@@ -360,24 +426,24 @@ for l in np.arange(0,len(clusters_sick)):
         axs_split[2][l].imshow(cluster_centre[...,2].astype(np.uint8),cmap='gray')
         axs_split[2][l].axis('off')
         plt.figure(fig_split.number),
-        plt.savefig(dir_probs + 'clusterCentres_'+disease+'_splitContrastEnhanced_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+        plt.savefig(dir_probs + 'clustercenters_'+disease+'_splitContrastEnhanced_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
         axs[l].imshow(cluster_centre.astype(np.uint8))
         axs[l].axis('off')
         axs[l].set_title(clusters_sick[l])
     plt.figure(fig.number),
-    plt.savefig(dir_probs + 'clusterCentres_'+disease+'_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+    plt.savefig(dir_probs + 'clustercenters_'+disease+'_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
     
 if colour_mode!='bnw':
     #plot sick and control clusters with the same intensity range
     max_value = []
     min_value = []
     for l in np.arange(0,len(clusters_control)):
-        cluster_centre = np.transpose((np.reshape(cluster_centres_control[l,:],(3,patch_size,patch_size))),(1,2,0))
+        cluster_centre = np.transpose((np.reshape(cluster_centers_control[l,:],(3,patch_size,patch_size))),(1,2,0))
         max_value += [cluster_centre.max()]
         min_value += [cluster_centre.min()]
         
     for l in np.arange(0,len(clusters_sick)):
-        cluster_centre = np.transpose((np.reshape(cluster_centres_sick[l,:],(3,patch_size,patch_size))),(1,2,0))
+        cluster_centre = np.transpose((np.reshape(cluster_centers_sick[l,:],(3,patch_size,patch_size))),(1,2,0))
         max_value += [cluster_centre.max()]
         min_value += [cluster_centre.min()]
        
@@ -385,9 +451,9 @@ if colour_mode!='bnw':
     range_min = min(min_value)
     
     fig_split, axs_split = plt.subplots(3,len(clusters_control), figsize=(len(clusters_control)*3,9),sharex=True, sharey=True)
-    fig_split.suptitle('Control centres split channels (fixed intensity range)')
+    fig_split.suptitle('Control centers split channels (fixed intensity range)')
     for l in np.arange(0,len(clusters_control)):
-        cluster_centre = np.transpose((np.reshape(cluster_centres_control[l,:],(3,patch_size,patch_size))),(1,2,0))
+        cluster_centre = np.transpose((np.reshape(cluster_centers_control[l,:],(3,patch_size,patch_size))),(1,2,0))
         im = np.zeros(cluster_centre.shape)
         im[...,0] = cluster_centre[...,0]
         axs_split[0][l].imshow(im.astype(np.uint8))
@@ -401,13 +467,13 @@ if colour_mode!='bnw':
         im[...,2] = cluster_centre[...,2]
         axs_split[2][l].imshow(im.astype(np.uint8))
         axs_split[2][l].axis('off')
-        plt.savefig(dir_probs + 'clusterCentres_control_split_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+        plt.savefig(dir_probs + 'clustercenters_control_split_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
         
         
     fig_split, axs_split = plt.subplots(3,len(clusters_sick), figsize=(len(clusters_sick)*3,9),sharex=True, sharey=True)
-    fig_split.suptitle('sick centres split channels (fixed intensity range)')
+    fig_split.suptitle('sick centers split channels (fixed intensity range)')
     for l in np.arange(0,len(clusters_sick)):
-        cluster_centre = np.transpose((np.reshape(cluster_centres_sick[l,:],(3,patch_size,patch_size))),(1,2,0))
+        cluster_centre = np.transpose((np.reshape(cluster_centers_sick[l,:],(3,patch_size,patch_size))),(1,2,0))
         im = np.zeros(cluster_centre.shape)
         im[...,0] = cluster_centre[...,0]
         axs_split[0][l].imshow(im.astype(np.uint8),vmin = range_min, vmax = range_max)
@@ -421,7 +487,7 @@ if colour_mode!='bnw':
         im[...,2] = cluster_centre[...,2]
         axs_split[2][l].imshow(im.astype(np.uint8),vmin = range_min, vmax = range_max)
         axs_split[2][l].axis('off')
-        plt.savefig(dir_probs + 'clusterCentres_'+disease+'_split_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)
+        plt.savefig(dir_probs + 'clustercenters_'+disease+'_split_%dclusters_%ddownscale_%dpatchsize.png'%(nr_clusters,1/sc_fac,patch_size), dpi=300)