X-Git-Url: http://git.euphorik.ch/?p=master-thesis.git;a=blobdiff_plain;f=Documentation%2Fma_autocount%2Fplasmodium-autocount-0.1.py;h=14ae21d8d7c0dfb6e53ce023ce190b6b8eec95b4;hp=b2fa3d2ac8c707f6066a1d78ca0ff30bdae6e2b1;hb=32ef4bbb0381f822c4df124595d401ad2f502c1b;hpb=c6fe9d606eff98fc73b75f23e02a9fd436458ed0 diff --git a/Documentation/ma_autocount/plasmodium-autocount-0.1.py b/Documentation/ma_autocount/plasmodium-autocount-0.1.py index b2fa3d2..14ae21d 100644 --- a/Documentation/ma_autocount/plasmodium-autocount-0.1.py +++ b/Documentation/ma_autocount/plasmodium-autocount-0.1.py @@ -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):