3 # Copyright (C) 2010 Monash University
4 # This program is free software; you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation; either version 2 of the License, or
7 # (at your option) any later version.
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
14 # You should have received a copy of the GNU General Public License
15 # along with this program; if not, write to the Free Software
16 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 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
63 Show_infection
= False
65 Show_hand_marks
= True
66 Show_computer_marks
= True
67 Show_cell_shapes
= False
77 # ==========================================
78 # Image enhancement and fg/bg discrimination
80 # Blur slightly to reduce noise, and unsharp-mask
82 # gaussian_blur(image, a) - b * gaussian_blur(image, c)
88 # Larger = enlarge background area, smaller = enlarge foreground area
89 Background_weight
= 0.915 #&median
92 # =============================
93 # Stain and infection detection
95 Dark_stain_level
= 0.839
99 Stain_spread_required
= 3.0
101 Infection_level
= 0.762
102 Infection_bg_blur
= 2.0
103 Infection_pixels_required
= 1
107 # ====================
108 # Cell size parameters
110 # Size of cells, for Hough transform
112 Hough_max_radius
= 20
114 # Portion of cell used for classification as infected/weird
115 # (may be abutting an infected cell, so be cautious)
116 # (as a proportion of the true radius)
117 Cell_inside_radius
= 0.9
119 # Do not allow another cell center to be detected within this radius
120 Cell_suppression_radius
= 1.25
123 # =======================
124 # Peculiar cell filtering
126 # White blood cells appear to be composed entirely of dark staining, filter them out
127 Max_dark_stain_proportion
= 0.728
129 # Minimum/maximum radius (as determined by Hough transform)
133 # Proportion of radius that center of mass can be away from Hough determined center
142 def offsets(inner_radius
, outer_radius
):
144 Generate a list of offsets for pixels within a circular region
147 key
= (inner_radius
,outer_radius
)
148 if key
not in _offset_cache
:
150 int_radius
= int(outer_radius
)+1
151 for y
in xrange(-int_radius
,int_radius
+1):
152 for x
in xrange(-int_radius
,int_radius
+1):
153 if inner_radius
*inner_radius
<= y
*y
+x
*x
<= outer_radius
*outer_radius
:
155 _offset_cache
[key
] = numpy
.array(result
)
156 return _offset_cache
[key
]
159 def clipped_coordinates(inner_radius
,outer_radius
, y
,x
, height
,width
):
161 Coordinates contained within a ring or circle, clipped to rectangular bounds.
164 coords
= offsets(inner_radius
,outer_radius
) + [[y
,x
]]
167 (coords
[:,0] < height
) &
169 (coords
[:,1] < width
)
174 def set(array
, y
, x
, value
):
176 Set an element of an array, if it is in bounds, otherwise silently do nothing
179 if y
>= 0 and x
>= 0 and y
< array
.shape
[0] and x
< array
.shape
[1]:
183 def mark_hand(array
, filename_in
):
185 If there is a file of hand marked cells from Cell Counting Aid, mark them.
188 hand_name
= os
.path
.splitext(filename_in
)[0] + '.txt'
189 if not os
.path
.exists(hand_name
): return
191 f
= open(hand_name
,'rU')
197 parts
= line
.rstrip('\n').split(',')
198 x
,y
,infected
= parts
[:3]
200 x
= int( int(x
)/1005.0*array
.shape
[1] +0.5)
201 y
= int( int(y
)/592.0*array
.shape
[0] +0.5)
202 if infected
== 'True':
206 infection_type
= parts
[3]
207 if infection_type
== '1':
208 color
= (196, 196, 0)
209 elif infection_type
== '3':
212 raise Exception('Unknown infection type')
215 for i
in (-2,-1,0,1,2):
216 set(array
,y
-3,x
+i
,color
)
217 set(array
,y
+3,x
+i
,color
)
218 set(array
,y
+i
,x
-3,color
)
219 set(array
,y
+i
,x
+3,color
)
222 def gaussian_blur(array
, amount
):
224 Gaussian blur each channel of a color image.
226 output
= numpy
.empty(array
.shape
, array
.dtype
)
227 for i
in xrange(array
.shape
[2]):
228 ndimage
.gaussian_filter(array
[:,:,i
],amount
,output
=output
[:,:,i
])
234 #sys.stdout.write(word+' ')
238 def process(filename_in
, filename_out
, filename_coords
, name
):
239 # print 'File: ', filename_in
242 i
= Image
.open(filename_in
)
244 print 'Processing', filename_in
, '->', filename_out
246 # Resize (mostly done for speed increase)
247 i
= i
.resize((Resize_width
,Resize_height
), Image
.ANTIALIAS
)
249 # Convert to numpy array
250 a
= numpy
.asarray(i
).astype('float64')
254 a_gray
= numpy
.sum(a
, 2)
256 # Make a copy of the array to doodle various diagnostic markings on
259 # If there are hand-marked cells, show them in the output image
261 mark_hand(a_output
, filename_in
)
264 # Denoise/unsharp mask
265 a_enhanced
= gaussian_blur(a
, Enhance_a
) - Enhance_b
*gaussian_blur(a
,Enhance_c
)
268 # Split into foreground and background using k-medians
271 a_raveled
= numpy
.reshape(a_enhanced
,(width
*height
,a_enhanced
.shape
[2]))
272 average
= numpy
.average(a_raveled
, axis
=0)
274 color_bg
= average
*1.5
275 color_fg
= average
*0.5
278 d_bg
= a_raveled
-color_bg
[None,:]
280 d_bg
= numpy
.sum(d_bg
,1)
281 d_fg
= a_raveled
-color_fg
[None,:]
283 d_fg
= numpy
.sum(d_fg
,1)
284 fg
= d_fg
*Background_weight
< d_bg
285 color_fg
= numpy
.median(a_raveled
[fg
,:],axis
=0)
286 color_bg
= numpy
.median(a_raveled
[~fg
,:],axis
=0)
288 fg
= numpy
.reshape(fg
, (height
,width
))
290 edge
= ndimage
.maximum_filter(fg
,size
=3) != ndimage
.minimum_filter(fg
,size
=3)
292 d_fg
= numpy
.sqrt(numpy
.reshape(d_fg
, (height
,width
)))
293 d_bg
= numpy
.sqrt(numpy
.reshape(d_bg
, (height
,width
)))
295 # Staining is taken as foreground pixels that differ markedly from the mean foreground color
296 # (also, they must be darker in the red and green channels)
297 # The background intensity is simply a useful basis for scaling
298 # (should be invariant under changes in brightness)
299 mag_bg
= numpy
.average(color_bg
)
302 (d_fg
> mag_bg
*Dark_stain_level
) &
303 (a_enhanced
[:,:,0] < color_fg
[0]) &
304 (a_enhanced
[:,:,1] < color_fg
[1]) &
309 fg_float
= fg
.astype('float64')
310 green
= a
[:,:,1] * fg_float
311 green_local_bg
= ndimage
.gaussian_filter(green
, Stain_bg_blur
)
312 green_local_bg_divisor
= ndimage
.gaussian_filter(fg_float
, Stain_bg_blur
)
314 (green
* green_local_bg_divisor
< green_local_bg
* Stain_level
) &
318 # Detect infection as lumpiness in the green channel (green is most affected by staining)
319 # Must also be within a stained region
320 green_local_bg
= ndimage
.gaussian_filter(green
, Infection_bg_blur
)
321 green_local_bg_divisor
= ndimage
.gaussian_filter(fg_float
, Infection_bg_blur
)
323 (green
* green_local_bg_divisor
< green_local_bg
* Infection_level
) &
327 # Show edges and infection on output image
329 a_output
[edge
,:] += 64
331 a_output
[stain
,2] += 128
333 a_output
[infection
,0] += 128
339 x_edge
= ndimage
.convolve(a_gray
,
343 y_edge
= ndimage
.convolve(a_gray
,
347 mag
= numpy
.sqrt(x_edge
*x_edge
+y_edge
*y_edge
)
351 # Vectorized for enhanced performance and incomprehensibility:
354 for y
in xrange(height
):
355 for x
in xrange(width
):
361 indexes
= numpy
.arange(len(ys
))
363 def do_hough(radius
):
364 hough
= numpy
.zeros((height
,width
), 'float64')
366 oys
= (ys
+ y_dir
[ys
,xs
]*radius
+ 0.5).astype('int')
367 oxs
= (xs
+ x_dir
[ys
,xs
]*radius
+ 0.5).astype('int')
374 for i
in indexes
[valid
]:
375 hough
[oys
[i
],oxs
[i
]] += 1.
377 #for y in xrange(height):
378 # for x in xrange(width):
380 # oy = y + y_dir[y,x]*radius
381 # ox = x + x_dir[y,x]*radius
384 # if ioy>0 and iox>0 and ioy<height and iox<width:
385 # hough[ioy,iox] += 1.
387 # Blur the transformed image a little
388 # (more blurring allows greater deviation from circularity)
389 # (if you alter this you will also need to alter thresh, below)
391 hough
= ndimage
.gaussian_filter(hough
, 1.0)
395 hough
= numpy
.zeros((height
,width
),'float64')
396 hough_radius
= numpy
.zeros((height
,width
),'int')
397 for i
in xrange(Hough_min_radius
,Hough_max_radius
+1):
398 #print 'Hough radius', i
399 this_hough
= do_hough(i
)
400 better
= this_hough
> hough
401 hough_radius
[better
] = i
402 hough
[better
] = this_hough
[better
]
404 # Show result of Hough transform in output image
405 #a_output[:,:,:] = hough[:,:,None]*20
407 # Find and classify peaks in Hough transform image
410 for y
in xrange(height
):
411 for x
in xrange(width
):
412 if hough
[y
,x
] > thresh
:
413 candidates
.append((hough
[y
,x
],y
,x
))
414 candidates
.sort(reverse
=True)
418 file_coords
= open(filename_coords
,'wt')
419 file_coords
.write('Autocount on %dx%d image\n' % (width
,height
))
421 suppress
= numpy
.zeros((height
,width
), 'bool')
422 claim_strength
= numpy
.zeros((height
,width
), 'float64') + 1e30
423 cells
= [ ] # [(y,x)]
427 for _
,y
,x
in candidates
:
428 radius
= hough_radius
[y
,x
]
430 # Make sure candidate is not near a previously detected cell,
431 # and does not abut the edge of the image
432 if not suppress
[y
,x
] and \
433 y
>= radius
and x
>= radius
and y
< height
-radius
and x
< width
-radius
:
434 coords
= clipped_coordinates(0,radius
*Cell_inside_radius
, y
,x
,height
,width
)
436 #Only include foreground
437 this_fg
= fg
[coords
[:,0],coords
[:,1]]
438 coords
= coords
[this_fg
]
440 this_strength
= ((coords
[:,0]-y
)**2 + (coords
[:,1]-x
)**2)
441 claim_strength
[coords
[:,0],coords
[:,1]] = numpy
.minimum(
442 claim_strength
[coords
[:,0],coords
[:,1]],
446 cells
.append((y
,x
, coords
,this_strength
))
449 # Suppress detection of cells too close to here
450 # (note: we do this irrespective of whether we treated this location as a hit)
452 coords
= clipped_coordinates(0,radius
*Cell_suppression_radius
, y
,x
,height
,width
)
453 suppress
[coords
[:,0],coords
[:,1]] = True
455 #for oy, ox in offsets(0, radius*Cell_suppression_radius): #suppression_offsets:
458 # if ny >= 0 and nx >= 0 and ny < height and nx < width:
459 # suppress[ny,nx] = True
461 for y
,x
, coords
, this_strength
in cells
:
462 radius
= hough_radius
[y
,x
]
464 if not len(coords
): continue # No fg pixels within radius
466 winners
= claim_strength
[coords
[:,0],coords
[:,1]] >= this_strength
467 coords
= coords
[winners
]
469 if not len(coords
): continue # Did not win any foreground pixels
473 dark_stain_pixels
= 0
481 for ny
, nx
in coords
:
491 if dark_stain
[ny
,nx
]:
492 dark_stain_pixels
+= 1
495 mass_y
= numpy
.average(coords
[:,0])
496 mass_x
= numpy
.average(coords
[:,1])
497 mass_offset
= numpy
.sqrt( (mass_y
-y
)**2 + (mass_x
-x
)**2 )
499 is_infected
= (infected_pixels
>= Infection_pixels_required
and stain_pixels
>= 1)
503 stain_x2
/stain_pixels
- (stain_x
/stain_pixels
)**2 +
504 stain_y2
/stain_pixels
- (stain_y
/stain_pixels
)**2
506 is_infected
= stain_spread
>= Stain_spread_required
509 (radius
> Max_cell_radius
) or
510 (radius
< Min_cell_radius
) or
511 (mass_offset
> Max_offcenter
* radius
) or
512 (dark_stain_pixels
>= fg_pixels
*Max_dark_stain_proportion
)
517 file_coords
.write('%d,%d,Rejected\n' % (x
,y
))
522 file_coords
.write('%d,%d,Infected\n' % (x
,y
))
526 file_coords
.write('%d,%d,Uninfected\n' % (x
,y
))
528 # Mark detected cell on output image
529 if Show_computer_marks
:
530 set(a_output
,y
,x
,color
)
531 for i
in xrange(1,3):
532 set(a_output
,y
-i
,x
,color
)
533 set(a_output
,y
+i
,x
,color
)
534 set(a_output
,y
,x
-i
,color
)
535 set(a_output
,y
,x
+i
,color
)
538 a_output
[coords
[:,0],coords
[:,1]] += [[ random
.random()*64+32 for i
in (0,1,2) ]]
544 # Save diagnostic output image
545 a_output
= numpy
.clip(a_output
,0,255)
546 i2
= Image
.fromarray(a_output
.astype('uint8'))
547 i2
.save(filename_out
)
551 def read_ac(filename
):
552 f
= open(filename
,'rU')
557 parts
= line
.rstrip().split(',')
558 if parts
[2] == 'Uninfected':
560 elif parts
[2] == 'Infected':
564 return n_total
, n_infected
567 if __name__
== '__main__':
571 name
, value
= args
.pop(0).split('=')
573 assert hasattr(sys
.modules
[__name__
], name
), 'Parameter named %s does not exist' % name
574 setattr(sys
.modules
[__name__
], name
, value
)
575 print name
, '=', value
578 print Help
% sys
.argv
[0]
581 if os
.path
.isfile(args
[0]):
582 in_dir
, filename
= os
.path
.split(args
[0])
583 filenames
= [ filename
]
588 filenames
= os
.listdir(in_dir
)
590 is_subprocess
= False
592 out_dir
= os
.path
.normpath(in_dir
) + '_output'
593 if not os
.path
.exists(out_dir
):
596 if not is_subprocess
:
597 for filename
in filenames
:
598 root
, ext
= os
.path
.splitext(filename
)
599 ac_filename
= os
.path
.join(out_dir
, root
+'.ac')
600 if os
.path
.exists(ac_filename
): os
.unlink(ac_filename
)
602 coord_filenames
= [ ]
606 for filename
in filenames
:
607 root
, ext
= os
.path
.splitext(filename
)
608 if ext
.lower() not in ('.bmp','.tif','.tiff','.png'): continue
610 filename_in
= os
.path
.join(in_dir
, filename
)
611 filename_out
= os
.path
.join(out_dir
, root
+'.tif')
612 filename_coords
= os
.path
.join(out_dir
, root
+'.ac')
614 if is_subprocess
or N_processors
== 1:
615 process(filename_in
, filename_out
, filename_coords
, root
)
617 while n_running
>= N_processors
:
620 os
.spawnl(os
.P_NOWAIT
, sys
.executable
, sys
.executable
, *(sys
.argv
[:-1]+[filename_in
]))
623 coord_filenames
.append((filename_coords
, root
))
629 if not is_subprocess
:
630 result_file
= open(os
.path
.join(out_dir
,'results.txt'),'wt')
632 for filename_coords
, root
in coord_filenames
:
633 n_total
, n_infected
= read_ac(filename_coords
)
634 print >> result_file
, root
, n_total
, n_infected