1 # -*- coding: utf-8 -*-
4 # Copyright (C) 2010 Monash University
5 # This program is free software; you can redistribute it and/or modify
6 # it under the terms of the GNU General Public License as published by
7 # the Free Software Foundation; either version 2 of the License, or
8 # (at your option) any later version.
10 # This program is distributed in the hope that it will be useful,
11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 # GNU General Public License for more details.
15 # You should have received a copy of the GNU General Public License
16 # along with this program; if not, write to the Free Software
17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 Identify red blood cells and red blood cells infected with plasmodium in a set of images.
27 import sys
, numpy
, os
, random
28 from scipy
import ndimage
34 python %s [parameter adjustments] <image_directory>
36 Images should be in TIFF, PNG, or BMP format.
38 Results will be placed in a new directory named <image_directory>_output.
40 If the images were annotated using Cell Counting Aid, the hand annotations
41 will be shown as circles in the output.
43 Parameter adjustments take the form
47 See the top of the source code for parameters and their default values.
48 These have been tuned for our microscope setup, optimal parameters for
49 your setup may differ.
53 # How many subprocesses to run at once
59 # What to show in diagnostic output images
62 Show_dark_stain
= False
64 Show_infection
= False
66 Show_hand_marks
= False
67 Show_computer_marks
= True
68 Show_cell_shapes
= True
76 # Tailles pour les images dans '08.10.2015'
77 Resize_width
= 1296 #680
78 Resize_height
= 972 #512
80 # Tailles pour les images dans '22.05.2015' (50%)
87 # Resize_height = 2885
90 # ==========================================
91 # Image enhancement and fg/bg discrimination
93 # Blur slightly to reduce noise, and unsharp-mask
95 # gaussian_blur(image, a) - b * gaussian_blur(image, c)
97 Enhance_a
= 2.0 # base: 1.0, 50x: 2.0, 100x: 4.0
98 Enhance_b
= 0.75 # base: 0.75
99 Enhance_c
= 40.0 # base: 20.0, 50x: 40.0, 100x: 80.0
101 # Larger = enlarge background area, smaller = enlarge foreground area
102 Background_weight
= 0.915 #&median
105 # =============================
106 # Stain and infection detection
108 Dark_stain_level
= 0.839
112 Stain_spread_required
= 3.0
114 Infection_level
= 0.762
115 Infection_bg_blur
= 2.0
116 Infection_pixels_required
= 1
120 # ====================
121 # Cell size parameters
123 # Size of cells, for Hough transform
124 Hough_min_radius
= 20 # 50x: 20, 100x: 40
125 Hough_max_radius
= 45 # 50x: 45, 100x: 90
127 # Portion of cell used for classification as infected/weird
128 # (may be abutting an infected cell, so be cautious)
129 # (as a proportion of the true radius)
130 Cell_inside_radius
= 0.9
132 # Do not allow another cell center to be detected within this radius
133 Cell_suppression_radius
= 0.8 # 1.05
136 # =======================
137 # Peculiar cell filtering
139 # White blood cells appear to be composed entirely of dark staining, filter them out
140 Max_dark_stain_proportion
= 0.728
142 # Minimum/maximum radius (as determined by Hough transform)
143 Min_cell_radius
= 25 # Base: 9, 50x: 25, 100x: 50
144 Max_cell_radius
= 40 # Base: 20, 50x: 40, 100x: 80
146 # Proportion of radius that center of mass can be away from Hough determined center
151 def offsets(inner_radius
, outer_radius
):
153 Generate a list of offsets for pixels within a circular region
156 key
= (inner_radius
,outer_radius
)
157 if key
not in _offset_cache
:
159 int_radius
= int(outer_radius
)+1
160 for y
in xrange(-int_radius
,int_radius
+1):
161 for x
in xrange(-int_radius
,int_radius
+1):
162 if inner_radius
*inner_radius
<= y
*y
+x
*x
<= outer_radius
*outer_radius
:
164 _offset_cache
[key
] = numpy
.array(result
)
165 return _offset_cache
[key
]
168 def clipped_coordinates(inner_radius
,outer_radius
, y
,x
, height
,width
):
170 Coordinates contained within a ring or circle, clipped to rectangular bounds.
173 coords
= offsets(inner_radius
,outer_radius
) + [[y
,x
]]
176 (coords
[:,0] < height
) &
178 (coords
[:,1] < width
)
183 def set(array
, y
, x
, value
):
185 Set an element of an array, if it is in bounds, otherwise silently do nothing
188 if y
>= 0 and x
>= 0 and y
< array
.shape
[0] and x
< array
.shape
[1]:
192 def mark_hand(array
, filename_in
):
194 If there is a file of hand marked cells from Cell Counting Aid, mark them.
197 hand_name
= os
.path
.splitext(filename_in
)[0] + '.txt'
198 if not os
.path
.exists(hand_name
): return
200 f
= open(hand_name
,'rU')
206 parts
= line
.rstrip('\n').split(',')
207 x
,y
,infected
= parts
[:3]
209 x
= int( int(x
)/1005.0*array
.shape
[1] +0.5)
210 y
= int( int(y
)/592.0*array
.shape
[0] +0.5)
211 if infected
== 'True':
215 infection_type
= parts
[3]
216 if infection_type
== '1':
217 color
= (196, 196, 0)
218 elif infection_type
== '3':
221 raise Exception('Unknown infection type')
224 for i
in (-2,-1,0,1,2):
225 set(array
,y
-3,x
+i
,color
)
226 set(array
,y
+3,x
+i
,color
)
227 set(array
,y
+i
,x
-3,color
)
228 set(array
,y
+i
,x
+3,color
)
231 def gaussian_blur(array
, amount
):
233 Gaussian blur each channel of a color image.
235 output
= numpy
.empty(array
.shape
, array
.dtype
)
236 for i
in xrange(array
.shape
[2]):
237 ndimage
.gaussian_filter(array
[:,:,i
],amount
, output
=output
[:,:,i
])
243 #sys.stdout.write(word+' ')
247 def process(filename_in
, filename_out
, filename_coords
, name
):
248 # print 'File: ', filename_in
251 i
= Image
.open(filename_in
)
253 print 'Processing', filename_in
, '->', filename_out
255 # Resize (mostly done for speed increase)
256 i
= i
.resize((Resize_width
,Resize_height
), Image
.ANTIALIAS
) if Resize
else i
258 # Convert to numpy array
259 a
= numpy
.asarray(i
).astype('float64') # [GB] Shape: (512, 680, 3)
263 a_gray
= numpy
.sum(a
, 2) # [GB] Somme de chaque composante RGB (non-normalisé)
265 # Make a copy of the array to doodle various diagnostic markings on
268 # If there are hand-marked cells, show them in the output image
270 mark_hand(a_output
, filename_in
)
273 # Denoise/unsharp mask
274 # [GB] (Concerne les trois channels RGB).
275 a_enhanced
= gaussian_blur(a
, Enhance_a
) - Enhance_b
* gaussian_blur(a
, Enhance_c
)
278 # Split into foreground and background using k-medians
281 # [GB] Chaque canal de couleur est linéarisé dans un tableau.
282 a_raveled
= numpy
.reshape(a_enhanced
,(width
*height
,a_enhanced
.shape
[2]))
284 # [GB] Moyenne de chaque canal.
285 average
= numpy
.average(a_raveled
, axis
=0)
287 # Initial guess of the medians
288 color_bg
= average
* 1.5 # [GB] Plus clair
289 color_fg
= average
* 0.5 # [GB] Plus sombre
291 # Fifty times: (0 to 4).
293 d_bg
= a_raveled
- color_bg
[None,:]
295 d_bg
= numpy
.sum(d_bg
, 1) # [GB] Somme des carrés des couleurs
296 d_fg
= a_raveled
- color_fg
[None,:]
298 d_fg
= numpy
.sum(d_fg
, 1)
299 fg
= d_fg
* Background_weight
< d_bg
300 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
301 color_bg
= numpy
.median(a_raveled
[~fg
,:],axis
=0)
303 fg
= numpy
.reshape(fg
, (height
,width
)) # [GB] Image binaire du premier plan.
305 edge
= ndimage
.maximum_filter(fg
,size
=3) != ndimage
.minimum_filter(fg
,size
=3) # [GB] image binaire des contours du premier plan.
307 # [GB] délinéarisation des tableaux résultant du k-median (niveau de gris).
308 d_fg
= numpy
.sqrt(numpy
.reshape(d_fg
, (height
,width
)))
309 d_bg
= numpy
.sqrt(numpy
.reshape(d_bg
, (height
,width
)))
311 # Staining is taken as foreground pixels that differ markedly from the mean foreground color
312 # (also, they must be darker in the red and green channels)
313 # The background intensity is simply a useful basis for scaling
314 # (should be invariant under changes in brightness)
315 mag_bg
= numpy
.average(color_bg
)
317 # [GB] Les "dark stain" :
318 # * sont d'une intensité plus élevée que la moyenne du bg multiplié par un facteur (Dark_stain_level)
319 # * les valeurs rouges sont plus petites que la mediane du premier-plan pour les rouges
320 # * les valeurs vertes sont plus petites que la médiane du premier-plan pour les verts
321 # * appartiennent au foreground
323 (d_fg
> mag_bg
*Dark_stain_level
) &
324 (a_enhanced
[:,:,0] < color_fg
[0]) & # [GB] Red channel
325 (a_enhanced
[:,:,1] < color_fg
[1]) & # [GB] Green channel
329 fg_float
= fg
.astype('float64')
330 green
= a
[:,:,1] * fg_float
331 green_local_bg
= ndimage
.gaussian_filter(green
, Stain_bg_blur
)
332 green_local_bg_divisor
= ndimage
.gaussian_filter(fg_float
, Stain_bg_blur
)
334 (green
* green_local_bg_divisor
< green_local_bg
* Stain_level
) &
338 # Detect infection as lumpiness in the green channel (green is most affected by staining)
339 # Must also be within a stained region
340 green_local_bg
= ndimage
.gaussian_filter(green
, Infection_bg_blur
)
341 green_local_bg_divisor
= ndimage
.gaussian_filter(fg_float
, Infection_bg_blur
)
343 (green
* green_local_bg_divisor
< green_local_bg
* Infection_level
) &
347 # Show edges and infection on output image
349 a_output
[edge
,:] += 64
351 a_output
[stain
,2] += 128
353 a_output
[infection
,0] += 128
355 a_output
[infection
,0] += 128
357 #a_output[stain, 1] += 128
363 x_edge
= ndimage
.convolve(a_gray
,
367 y_edge
= ndimage
.convolve(a_gray
,
371 mag
= numpy
.sqrt(x_edge
*x_edge
+y_edge
*y_edge
)
375 # Vectorized for enhanced performance and incomprehensibility:
378 for y
in xrange(height
):
379 for x
in xrange(width
):
385 indexes
= numpy
.arange(len(ys
))
387 def do_hough(radius
):
388 hough
= numpy
.zeros((height
,width
), 'float64')
390 oys
= (ys
+ y_dir
[ys
,xs
]*radius
+ 0.5).astype('int')
391 oxs
= (xs
+ x_dir
[ys
,xs
]*radius
+ 0.5).astype('int')
398 for i
in indexes
[valid
]:
399 hough
[oys
[i
],oxs
[i
]] += 1.
401 #for y in xrange(height):
402 # for x in xrange(width):
404 # oy = y + y_dir[y,x]*radius
405 # ox = x + x_dir[y,x]*radius
408 # if ioy>0 and iox>0 and ioy<height and iox<width:
409 # hough[ioy,iox] += 1.
411 # Blur the transformed image a little
412 # (more blurring allows greater deviation from circularity)
413 # (if you alter this you will also need to alter thresh, below)
415 hough
= ndimage
.gaussian_filter(hough
, 1.5) # 1.0
419 hough
= numpy
.zeros((height
,width
),'float64')
420 hough_radius
= numpy
.zeros((height
,width
),'int')
421 for i
in xrange(Hough_min_radius
,Hough_max_radius
+1):
422 #print 'Hough radius', i
423 this_hough
= do_hough(i
)
424 better
= this_hough
> hough
425 hough_radius
[better
] = i
426 hough
[better
] = this_hough
[better
]
428 # Show result of Hough transform in output image
429 #a_output[:,:,:] = hough[:,:,None]*20
431 # Find and classify peaks in Hough transform image
432 thresh
= 0.8 #Hmm 1.5
434 for y
in xrange(height
):
435 for x
in xrange(width
):
436 if hough
[y
, x
] > thresh
:
437 candidates
.append((hough
[y
,x
], y
, x
))
438 candidates
.sort(reverse
= True)
443 file_coords
= open(filename_coords
,'wt')
444 file_coords
.write('Autocount on %dx%d image\n' % (width
,height
))
446 suppress
= numpy
.zeros((height
,width
), 'bool')
447 claim_strength
= numpy
.zeros((height
,width
), 'float64') + 1e30
451 for _
,y
,x
in candidates
:
452 radius
= hough_radius
[y
,x
]
454 # Make sure candidate is not near a previously detected cell,
455 # and does not abut the edge of the image
456 if not suppress
[y
,x
] and \
457 y
>= radius
and x
>= radius
and y
< height
-radius
and x
< width
-radius
:
458 coords
= clipped_coordinates(0, radius
* Cell_inside_radius
, y
, x
, height
, width
)
460 #Only include foreground
461 this_fg
= fg
[coords
[:,0], coords
[:,1]] # [GB] Booleans
462 coords
= coords
[this_fg
]
464 this_strength
= ((coords
[:,0]-y
)**2 + (coords
[:,1]-x
)**2)
465 claim_strength
[coords
[:,0],coords
[:,1]] = numpy
.minimum(
466 claim_strength
[coords
[:,0],coords
[:,1]],
470 cells
.append((y
, x
, coords
, this_strength
))
473 # Suppress detection of cells too close to here
474 # (note: we do this irrespective of whether we treated this location as a hit)
475 coords
= clipped_coordinates(0,radius
* Cell_suppression_radius
, y
,x
,height
,width
)
476 suppress
[coords
[:,0],coords
[:,1]] = True
478 #for oy, ox in offsets(0, radius*Cell_suppression_radius): #suppression_offsets:
481 # if ny >= 0 and nx >= 0 and ny < height and nx < width:
482 # suppress[ny,nx] = True
484 for y
,x
, coords
, this_strength
in cells
:
485 radius
= hough_radius
[y
,x
]
487 if not len(coords
): continue # No fg pixels within radius
489 winners
= claim_strength
[coords
[:,0],coords
[:,1]] >= this_strength
490 coords
= coords
[winners
]
492 if not len(coords
): continue # Did not win any foreground pixels
496 dark_stain_pixels
= 0
504 for ny
, nx
in coords
:
514 if dark_stain
[ny
,nx
]:
515 dark_stain_pixels
+= 1
518 mass_y
= numpy
.average(coords
[:,0])
519 mass_x
= numpy
.average(coords
[:,1])
520 mass_offset
= numpy
.sqrt( (mass_y
-y
)**2 + (mass_x
-x
)**2 )
522 is_infected
= (infected_pixels
>= Infection_pixels_required
and stain_pixels
>= 1)
526 stain_x2
/stain_pixels
- (stain_x
/stain_pixels
)**2 +
527 stain_y2
/stain_pixels
- (stain_y
/stain_pixels
)**2
529 is_infected
= stain_spread
>= Stain_spread_required
532 (radius
> Max_cell_radius
) or
533 (radius
< Min_cell_radius
) or
534 (mass_offset
> Max_offcenter
* radius
) or
535 (dark_stain_pixels
>= fg_pixels
*Max_dark_stain_proportion
)
540 file_coords
.write('%d,%d,Rejected\n' % (x
,y
))
545 file_coords
.write('%d,%d,Infected\n' % (x
,y
))
549 file_coords
.write('%d,%d,Uninfected\n' % (x
,y
))
552 a_output
[coords
[:,0],coords
[:,1]] += [[ random
.random()*64+32 for i
in (0,1,2) ]]
554 # Mark detected cell on output image
555 if Show_computer_marks
:
556 set(a_output
, y
, x
, color
)
557 for i
in xrange(1,3):
558 set(a_output
,y
-i
,x
,color
)
559 set(a_output
,y
+i
,x
,color
)
560 set(a_output
,y
,x
-i
,color
)
561 set(a_output
,y
,x
+i
,color
)
568 # Save diagnostic output image
569 a_output
= numpy
.clip(a_output
,0,255)
570 i2
= Image
.fromarray(a_output
.astype('uint8'))
571 i2
.save(filename_out
)
573 thresh_output
= numpy
.clip(fg
, 0, 255)
574 i3
= Image
.fromarray(thresh_output
.astype('uint8'))
575 i3
.save(filename_out
+ " - thresholded.tiff")
579 def read_ac(filename
):
580 f
= open(filename
,'rU')
585 parts
= line
.rstrip().split(',')
586 if parts
[2] == 'Uninfected':
588 elif parts
[2] == 'Infected':
592 return n_total
, n_infected
595 if __name__
== '__main__':
599 name
, value
= args
.pop(0).split('=')
601 assert hasattr(sys
.modules
[__name__
], name
), 'Parameter named %s does not exist' % name
602 setattr(sys
.modules
[__name__
], name
, value
)
603 print name
, '=', value
606 print Help
% sys
.argv
[0]
609 if os
.path
.isfile(args
[0]):
610 in_dir
, filename
= os
.path
.split(args
[0])
611 filenames
= [ filename
]
616 filenames
= os
.listdir(in_dir
)
618 is_subprocess
= False
620 out_dir
= os
.path
.normpath(in_dir
) + '_output'
621 if not os
.path
.exists(out_dir
):
624 if not is_subprocess
:
625 for filename
in filenames
:
626 root
, ext
= os
.path
.splitext(filename
)
627 ac_filename
= os
.path
.join(out_dir
, root
+'.ac')
628 if os
.path
.exists(ac_filename
): os
.unlink(ac_filename
)
630 coord_filenames
= [ ]
634 for filename
in filenames
:
635 root
, ext
= os
.path
.splitext(filename
)
636 if ext
.lower() not in ('.bmp','.tif','.tiff','.png'): continue
638 filename_in
= os
.path
.join(in_dir
, filename
)
639 filename_out
= os
.path
.join(out_dir
, root
+'.tif')
640 filename_coords
= os
.path
.join(out_dir
, root
+'.ac')
642 if is_subprocess
or N_processors
== 1:
643 process(filename_in
, filename_out
, filename_coords
, root
)
645 while n_running
>= N_processors
:
648 os
.spawnl(os
.P_NOWAIT
, sys
.executable
, sys
.executable
, *(sys
.argv
[:-1]+[filename_in
]))
651 coord_filenames
.append((filename_coords
, root
))
657 if not is_subprocess
:
658 result_file
= open(os
.path
.join(out_dir
,'results.txt'),'wt')
660 for filename_coords
, root
in coord_filenames
:
661 n_total
, n_infected
= read_ac(filename_coords
)
662 print >> result_file
, root
, n_total
, n_infected