Script modification for testing.
[master-thesis.git] / Documentation / ma_autocount / plasmodium-autocount-0.1.py
index b2fa3d2..14ae21d 100644 (file)
@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
 #!/usr/bin/env python
 
 #    Copyright (C) 2010 Monash University
@@ -23,7 +24,6 @@ Identify red blood cells and red blood cells infected with plasmodium in a set o
 """
 
 from PIL import Image
-
 import sys, numpy, os, random
 from scipy import ndimage
 
@@ -59,19 +59,32 @@ N_processors = 1
 # What to show in diagnostic output images
 
 Show_edges = False
+Show_dark_stain = False
 Show_stain = False
 Show_infection = False
 
-Show_hand_marks = True
+Show_hand_marks = False
 Show_computer_marks = True
-Show_cell_shapes = False
+Show_cell_shapes = True
 
 
 # =================
 # Image scaling
 
-Resize_width = 680
-Resize_height = 512
+Resize = False
+
+# Tailles pour les images dans '08.10.2015'
+Resize_width = 1296 #680
+Resize_height = 972 #512
+
+# Tailles pour les images dans '22.05.2015' (50%)
+#Resize_width = 1296
+#Resize_height = 972
+
+
+# [GB] Hemato
+# Resize_width = 2864
+# Resize_height = 2885
 
 
 # ==========================================
@@ -81,9 +94,9 @@ Resize_height = 512
 #
 # gaussian_blur(image, a) - b * gaussian_blur(image, c)
 
-Enhance_a = 1.0
-Enhance_b = 0.75
-Enhance_c = 20.0
+Enhance_a = 2.0  # base: 1.0,  50x: 2.0, 100x: 4.0
+Enhance_b = 0.75 # base: 0.75
+Enhance_c = 40.0 # base: 20.0, 50x: 40.0, 100x: 80.0
 
 # Larger = enlarge background area, smaller = enlarge foreground area
 Background_weight = 0.915 #&median
@@ -108,8 +121,8 @@ Infection_pixels_required = 1
 # Cell size parameters
 
 # Size of cells, for Hough transform
-Hough_min_radius = 5  
-Hough_max_radius = 20 
+Hough_min_radius = 20 # 50x: 20, 100x: 40
+Hough_max_radius = 45 # 50x: 45, 100x: 90
 
 # Portion of cell used for classification as infected/weird
 # (may be abutting an infected cell, so be cautious)
@@ -117,7 +130,7 @@ Hough_max_radius = 20
 Cell_inside_radius = 0.9
 
 # Do not allow another cell center to be detected within this radius
-Cell_suppression_radius = 1.25
+Cell_suppression_radius = 0.8 # 1.05
 
 
 # =======================
@@ -127,17 +140,13 @@ Cell_suppression_radius = 1.25
 Max_dark_stain_proportion = 0.728
 
 # Minimum/maximum radius (as determined by Hough transform)
-Min_cell_radius = 9
-Max_cell_radius = 20
+Min_cell_radius = 25 # Base: 9,  50x: 25, 100x: 50
+Max_cell_radius = 40 # Base: 20, 50x: 40, 100x: 80
 
 # Proportion of radius that center of mass can be away from Hough determined center
 Max_offcenter = 0.5
 
 
-
-
-
-
 _offset_cache = { }
 def offsets(inner_radius, outer_radius):
     """
@@ -225,7 +234,7 @@ def gaussian_blur(array, amount):
     """
     output = numpy.empty(array.shape, array.dtype)
     for i in xrange(array.shape[2]):
-        ndimage.gaussian_filter(array[:,:,i],amount,output=output[:,:,i])
+        ndimage.gaussian_filter(array[:,:,i],amount, output=output[:,:,i])
     return output
 
 
@@ -244,14 +253,14 @@ def process(filename_in, filename_out, filename_coords, name):
     print 'Processing', filename_in, '->', filename_out
 
     # Resize (mostly done for speed increase)
-    i = i.resize((Resize_width,Resize_height), Image.ANTIALIAS)
+    i = i.resize((Resize_width,Resize_height), Image.ANTIALIAS) if Resize else i
 
     # Convert to numpy array
-    a = numpy.asarray(i).astype('float64')
+    a = numpy.asarray(i).astype('float64') # [GB] Shape: (512, 680, 3)
     height = a.shape[0]
     width = a.shape[1]
 
-    a_gray = numpy.sum(a, 2)
+    a_gray = numpy.sum(a, 2) # [GB] Somme de chaque composante RGB (non-normalisé)
     
     # Make a copy of the array to doodle various diagnostic markings on
     a_output = a.copy()    
@@ -262,33 +271,40 @@ def process(filename_in, filename_out, filename_coords, name):
 
 
     # Denoise/unsharp mask
-    a_enhanced = gaussian_blur(a, Enhance_a) - Enhance_b*gaussian_blur(a,Enhance_c)
-    
+    # [GB] (Concerne les trois channels RGB).
+    a_enhanced = gaussian_blur(a, Enhance_a) - Enhance_b * gaussian_blur(a, Enhance_c)
+        
     
     # Split into foreground and background using k-medians
     status('fg/bg')
     
+    # [GB] Chaque canal de couleur est linéarisé dans un tableau.
     a_raveled = numpy.reshape(a_enhanced,(width*height,a_enhanced.shape[2]))
-    average = numpy.average(a_raveled, axis=0)
-    # Initial guess
-    color_bg = average*1.5
-    color_fg = average*0.5
     
+    # [GB] Moyenne de chaque canal.
+    average = numpy.average(a_raveled, axis=0) 
+    
+    # Initial guess of the medians
+    color_bg = average * 1.5 # [GB] Plus clair
+    color_fg = average * 0.5 # [GB] Plus sombre
+    
+    # Fifty times: (0 to 4).
     for i in xrange(5):
-        d_bg = a_raveled-color_bg[None,:]
+        d_bg = a_raveled - color_bg[None,:]
         d_bg *= d_bg
-        d_bg = numpy.sum(d_bg,1)
-        d_fg = a_raveled-color_fg[None,:]
+        d_bg = numpy.sum(d_bg, 1)  # [GB] Somme des carrés des couleurs
+        d_fg = a_raveled - color_fg[None,:]
         d_fg *= d_fg
-        d_fg = numpy.sum(d_fg,1)
-        fg = d_fg*Background_weight < d_bg
-        color_fg = numpy.median(a_raveled[fg,:],axis=0)
+        d_fg = numpy.sum(d_fg, 1)
+        fg = d_fg * Background_weight < d_bg
+        color_fg = numpy.median(a_raveled[fg,:],axis=0) # [GB] axis = 0 -> par colonne (ordonnées). Le résultat est alors 3 valeurs, une pour chaque couleur
         color_bg = numpy.median(a_raveled[~fg,:],axis=0)
-        
-    fg = numpy.reshape(fg, (height,width))
+
+    fg = numpy.reshape(fg, (height,width)) # [GB] Image binaire du premier plan.
     
-    edge = ndimage.maximum_filter(fg,size=3) != ndimage.minimum_filter(fg,size=3)
+    edge = ndimage.maximum_filter(fg,size=3) != ndimage.minimum_filter(fg,size=3) # [GB] image binaire des contours du premier plan.
 
+    # [GB] délinéarisation des tableaux résultant du k-median (niveau de gris).
     d_fg = numpy.sqrt(numpy.reshape(d_fg, (height,width)))
     d_bg = numpy.sqrt(numpy.reshape(d_bg, (height,width)))
     
@@ -298,14 +314,18 @@ def process(filename_in, filename_out, filename_coords, name):
     # (should be invariant under changes in brightness) 
     mag_bg = numpy.average(color_bg)
     
+    # [GB] Les "dark stain" :
+    # * sont d'une intensité plus élevée que la moyenne du bg multiplié par un facteur (Dark_stain_level)
+    # * les valeurs rouges sont plus petites que la mediane du premier-plan pour les rouges
+    # * les valeurs vertes sont plus petites que la médiane du premier-plan pour les verts
+    # * appartiennent au foreground
     dark_stain = (
         (d_fg > mag_bg*Dark_stain_level) &
-        (a_enhanced[:,:,0] < color_fg[0]) & 
-        (a_enhanced[:,:,1] < color_fg[1]) &
+        (a_enhanced[:,:,0] < color_fg[0]) & # [GB] Red channel
+        (a_enhanced[:,:,1] < color_fg[1]) & # [GB] Green channel
         fg
     )
-    
-    
+
     fg_float = fg.astype('float64')
     green = a[:,:,1] * fg_float
     green_local_bg = ndimage.gaussian_filter(green, Stain_bg_blur)
@@ -331,6 +351,10 @@ def process(filename_in, filename_out, filename_coords, name):
         a_output[stain,2] += 128
     if Show_infection:
         a_output[infection,0] += 128
+    if Show_dark_stain:
+        a_output[infection,0] += 128
+        
+    #a_output[stain, 1] += 128
         
         
     # Hough transform
@@ -388,7 +412,7 @@ def process(filename_in, filename_out, filename_coords, name):
         # (more blurring allows greater deviation from circularity)
         # (if you alter this you will also need to alter thresh, below)
         
-        hough = ndimage.gaussian_filter(hough, 1.0)
+        hough = ndimage.gaussian_filter(hough, 1.5) # 1.0
         
         return hough
 
@@ -405,13 +429,14 @@ def process(filename_in, filename_out, filename_coords, name):
     #a_output[:,:,:] = hough[:,:,None]*20
 
     # Find and classify peaks in Hough transform image
-    thresh = 1.5 #Hmm
+    thresh = 0.8 #Hmm 1.5
     candidates = [ ]
     for y in xrange(height):
         for x in xrange(width):
-            if hough[y,x] > thresh:
-                candidates.append((hough[y,x],y,x))
-    candidates.sort(reverse=True)
+            if hough[y, x] > thresh:
+                candidates.append((hough[y,x], y, x))
+    candidates.sort(reverse = True)
+        
     
     status('classify')
     
@@ -420,21 +445,20 @@ def process(filename_in, filename_out, filename_coords, name):
     
     suppress = numpy.zeros((height,width), 'bool')
     claim_strength = numpy.zeros((height,width), 'float64') + 1e30
-    cells = [ ] # [(y,x)]
+    cells = [] # [(y,x)]
     n = 0
     n_infected = 0
-    n_late = 0
     for _,y,x in candidates:
         radius = hough_radius[y,x]        
-
+        
         # Make sure candidate is not near a previously detected cell,
         # and does not abut the edge of the image        
         if not suppress[y,x] and \
            y >= radius and x >= radius and y < height-radius and x < width-radius:
-            coords = clipped_coordinates(0,radius*Cell_inside_radius, y,x,height,width)
+            coords = clipped_coordinates(0, radius * Cell_inside_radius, y, x, height, width)
             
             #Only include foreground
-            this_fg = fg[coords[:,0],coords[:,1]]
+            this_fg = fg[coords[:,0], coords[:,1]] # [GB] Booleans
             coords = coords[this_fg]
             
             this_strength = ((coords[:,0]-y)**2 + (coords[:,1]-x)**2)
@@ -443,13 +467,12 @@ def process(filename_in, filename_out, filename_coords, name):
                 this_strength
             )
             
-            cells.append((y,x, coords,this_strength))
+            cells.append((y, x, coords, this_strength))
 
-           
+       
         # Suppress detection of cells too close to here
-        # (note: we do this irrespective of whether we treated this location as a hit)
-        
-        coords = clipped_coordinates(0,radius*Cell_suppression_radius, y,x,height,width)
+        # (note: we do this irrespective of whether we treated this location as a hit)        
+        coords = clipped_coordinates(0,radius * Cell_suppression_radius, y,x,height,width)
         suppress[coords[:,0],coords[:,1]] = True
         
         #for oy, ox in offsets(0, radius*Cell_suppression_radius): #suppression_offsets:
@@ -525,17 +548,18 @@ def process(filename_in, filename_out, filename_coords, name):
             color = (0,0,255)
             file_coords.write('%d,%d,Uninfected\n' % (x,y))
         
+        if Show_cell_shapes:            
+            a_output[coords[:,0],coords[:,1]] += [[ random.random()*64+32 for i in (0,1,2) ]]
+            
         # Mark detected cell on output image
         if Show_computer_marks:
-            set(a_output,y,x,color)
+            set(a_output, y, x, color)
             for i in xrange(1,3):
                 set(a_output,y-i,x,color)  
                 set(a_output,y+i,x,color)  
                 set(a_output,y,x-i,color)  
                 set(a_output,y,x+i,color)
 
-        if Show_cell_shapes:            
-            a_output[coords[:,0],coords[:,1]] += [[ random.random()*64+32 for i in (0,1,2) ]]
             
         
     
@@ -546,6 +570,10 @@ def process(filename_in, filename_out, filename_coords, name):
     i2 = Image.fromarray(a_output.astype('uint8'))    
     i2.save(filename_out)
     
+    thresh_output = numpy.clip(fg, 0, 255)
+    i3 = Image.fromarray(thresh_output.astype('uint8'))    
+    i3.save(filename_out + " - thresholded.tiff")
+    
 
 
 def read_ac(filename):