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
63 Show_infection
= False
65 Show_hand_marks
= False
66 Show_computer_marks
= True
67 Show_cell_shapes
= True
75 # Tailles pour les images dans '08.10.2015'
76 Resize_width
= 1296 #680
77 Resize_height
= 972 #512
79 # Tailles pour les images dans '22.05.2015' (50%)
86 # Resize_height = 2885
89 # ==========================================
90 # Image enhancement and fg/bg discrimination
92 # Blur slightly to reduce noise, and unsharp-mask
94 # gaussian_blur(image, a) - b * gaussian_blur(image, c)
100 # Larger = enlarge background area, smaller = enlarge foreground area
101 Background_weight
= 0.915 #&median
104 # =============================
105 # Stain and infection detection
107 Dark_stain_level
= 0.839
111 Stain_spread_required
= 3.0
113 Infection_level
= 0.762
114 Infection_bg_blur
= 2.0
115 Infection_pixels_required
= 1
119 # ====================
120 # Cell size parameters
122 # Size of cells, for Hough transform
124 Hough_max_radius
= 32 # 20
126 # Portion of cell used for classification as infected/weird
127 # (may be abutting an infected cell, so be cautious)
128 # (as a proportion of the true radius)
129 Cell_inside_radius
= 0.9
131 # Do not allow another cell center to be detected within this radius
132 Cell_suppression_radius
= 1.05 # 1.25
135 # =======================
136 # Peculiar cell filtering
138 # White blood cells appear to be composed entirely of dark staining, filter them out
139 Max_dark_stain_proportion
= 0.728
141 # Minimum/maximum radius (as determined by Hough transform)
142 Min_cell_radius
= 11 # 9
143 Max_cell_radius
= 32 # 20
145 # Proportion of radius that center of mass can be away from Hough determined center
150 def offsets(inner_radius
, outer_radius
):
152 Generate a list of offsets for pixels within a circular region
155 key
= (inner_radius
,outer_radius
)
156 if key
not in _offset_cache
:
158 int_radius
= int(outer_radius
)+1
159 for y
in xrange(-int_radius
,int_radius
+1):
160 for x
in xrange(-int_radius
,int_radius
+1):
161 if inner_radius
*inner_radius
<= y
*y
+x
*x
<= outer_radius
*outer_radius
:
163 _offset_cache
[key
] = numpy
.array(result
)
164 return _offset_cache
[key
]
167 def clipped_coordinates(inner_radius
,outer_radius
, y
,x
, height
,width
):
169 Coordinates contained within a ring or circle, clipped to rectangular bounds.
172 coords
= offsets(inner_radius
,outer_radius
) + [[y
,x
]]
175 (coords
[:,0] < height
) &
177 (coords
[:,1] < width
)
182 def set(array
, y
, x
, value
):
184 Set an element of an array, if it is in bounds, otherwise silently do nothing
187 if y
>= 0 and x
>= 0 and y
< array
.shape
[0] and x
< array
.shape
[1]:
191 def mark_hand(array
, filename_in
):
193 If there is a file of hand marked cells from Cell Counting Aid, mark them.
196 hand_name
= os
.path
.splitext(filename_in
)[0] + '.txt'
197 if not os
.path
.exists(hand_name
): return
199 f
= open(hand_name
,'rU')
205 parts
= line
.rstrip('\n').split(',')
206 x
,y
,infected
= parts
[:3]
208 x
= int( int(x
)/1005.0*array
.shape
[1] +0.5)
209 y
= int( int(y
)/592.0*array
.shape
[0] +0.5)
210 if infected
== 'True':
214 infection_type
= parts
[3]
215 if infection_type
== '1':
216 color
= (196, 196, 0)
217 elif infection_type
== '3':
220 raise Exception('Unknown infection type')
223 for i
in (-2,-1,0,1,2):
224 set(array
,y
-3,x
+i
,color
)
225 set(array
,y
+3,x
+i
,color
)
226 set(array
,y
+i
,x
-3,color
)
227 set(array
,y
+i
,x
+3,color
)
230 def gaussian_blur(array
, amount
):
232 Gaussian blur each channel of a color image.
234 output
= numpy
.empty(array
.shape
, array
.dtype
)
235 for i
in xrange(array
.shape
[2]):
236 ndimage
.gaussian_filter(array
[:,:,i
],amount
,output
=output
[:,:,i
])
242 #sys.stdout.write(word+' ')
246 def process(filename_in
, filename_out
, filename_coords
, name
):
247 # print 'File: ', filename_in
250 i
= Image
.open(filename_in
)
252 print 'Processing', filename_in
, '->', filename_out
254 # Resize (mostly done for speed increase)
255 i
= i
.resize((Resize_width
,Resize_height
), Image
.ANTIALIAS
) if Resize
else i
257 # Convert to numpy array
258 a
= numpy
.asarray(i
).astype('float64') # [GB] Shape: (512, 680, 3)
262 a_gray
= numpy
.sum(a
, 2) # [GB] Somme de chaque composante RGB (non-normalisé)
264 # Make a copy of the array to doodle various diagnostic markings on
267 # If there are hand-marked cells, show them in the output image
269 mark_hand(a_output
, filename_in
)
272 # Denoise/unsharp mask
273 # [GB] (Concerne les trois channels RGB).
274 a_enhanced
= gaussian_blur(a
, Enhance_a
) - Enhance_b
* gaussian_blur(a
, Enhance_c
)
276 # Split into foreground and background using k-medians
279 # [GB] Chaque canal de couleur est linéarisé dans un tableau.
280 a_raveled
= numpy
.reshape(a_enhanced
,(width
*height
,a_enhanced
.shape
[2]))
282 # [GB] Moyenne de chaque canal.
283 average
= numpy
.average(a_raveled
, axis
=0)
285 # Initial guess of the medians
286 color_bg
= average
* 1.5 # [GB] Plus clair
287 color_fg
= average
* 0.5 # [GB] Plus sombre
289 # Fifty times: (0 to 4).
291 d_bg
= a_raveled
- color_bg
[None,:]
293 d_bg
= numpy
.sum(d_bg
, 1) # [GB] Somme des carrés des couleurs
294 d_fg
= a_raveled
- color_fg
[None,:]
296 d_fg
= numpy
.sum(d_fg
, 1)
297 fg
= d_fg
* Background_weight
< d_bg
298 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
299 color_bg
= numpy
.median(a_raveled
[~fg
,:],axis
=0)
301 fg
= numpy
.reshape(fg
, (height
,width
)) # [GB] Image binaire du premier plan.
303 edge
= ndimage
.maximum_filter(fg
,size
=3) != ndimage
.minimum_filter(fg
,size
=3) # [GB] image binaire des contours du premier plan.
305 # [GB] délinéarisation des tableaux résultant du k-median (niveau de gris).
306 d_fg
= numpy
.sqrt(numpy
.reshape(d_fg
, (height
,width
)))
307 d_bg
= numpy
.sqrt(numpy
.reshape(d_bg
, (height
,width
)))
309 # Staining is taken as foreground pixels that differ markedly from the mean foreground color
310 # (also, they must be darker in the red and green channels)
311 # The background intensity is simply a useful basis for scaling
312 # (should be invariant under changes in brightness)
313 mag_bg
= numpy
.average(color_bg
)
315 # [GB] Les "dark stain" :
316 # * sont d'une intensité plus élevée que la moyenne du bg multiplié par un facteur (Dark_stain_level)
317 # * les valeurs rouges sont plus petites que la mediane du premier-plan pour les rouges
318 # * les valeurs vertes sont plus petites que la médiane du premier-plan pour les verts
319 # * appartiennent au foreground
321 (d_fg
> mag_bg
*Dark_stain_level
) &
322 (a_enhanced
[:,:,0] < color_fg
[0]) & # [GB] Red channel
323 (a_enhanced
[:,:,1] < color_fg
[1]) & # [GB] Green channel
328 fg_float
= fg
.astype('float64')
329 green
= a
[:,:,1] * fg_float
330 green_local_bg
= ndimage
.gaussian_filter(green
, Stain_bg_blur
)
331 green_local_bg_divisor
= ndimage
.gaussian_filter(fg_float
, Stain_bg_blur
)
333 (green
* green_local_bg_divisor
< green_local_bg
* Stain_level
) &
337 # Detect infection as lumpiness in the green channel (green is most affected by staining)
338 # Must also be within a stained region
339 green_local_bg
= ndimage
.gaussian_filter(green
, Infection_bg_blur
)
340 green_local_bg_divisor
= ndimage
.gaussian_filter(fg_float
, Infection_bg_blur
)
342 (green
* green_local_bg_divisor
< green_local_bg
* Infection_level
) &
346 # Show edges and infection on output image
348 a_output
[edge
,:] += 64
350 a_output
[stain
,2] += 128
352 a_output
[infection
,0] += 128
354 #a_output[stain, 1] += 128
360 x_edge
= ndimage
.convolve(a_gray
,
364 y_edge
= ndimage
.convolve(a_gray
,
368 mag
= numpy
.sqrt(x_edge
*x_edge
+y_edge
*y_edge
)
372 # Vectorized for enhanced performance and incomprehensibility:
375 for y
in xrange(height
):
376 for x
in xrange(width
):
382 indexes
= numpy
.arange(len(ys
))
384 def do_hough(radius
):
385 hough
= numpy
.zeros((height
,width
), 'float64')
387 oys
= (ys
+ y_dir
[ys
,xs
]*radius
+ 0.5).astype('int')
388 oxs
= (xs
+ x_dir
[ys
,xs
]*radius
+ 0.5).astype('int')
395 for i
in indexes
[valid
]:
396 hough
[oys
[i
],oxs
[i
]] += 1.
398 #for y in xrange(height):
399 # for x in xrange(width):
401 # oy = y + y_dir[y,x]*radius
402 # ox = x + x_dir[y,x]*radius
405 # if ioy>0 and iox>0 and ioy<height and iox<width:
406 # hough[ioy,iox] += 1.
408 # Blur the transformed image a little
409 # (more blurring allows greater deviation from circularity)
410 # (if you alter this you will also need to alter thresh, below)
412 hough
= ndimage
.gaussian_filter(hough
, 1.0)
416 hough
= numpy
.zeros((height
,width
),'float64')
417 hough_radius
= numpy
.zeros((height
,width
),'int')
418 for i
in xrange(Hough_min_radius
,Hough_max_radius
+1):
419 #print 'Hough radius', i
420 this_hough
= do_hough(i
)
421 better
= this_hough
> hough
422 hough_radius
[better
] = i
423 hough
[better
] = this_hough
[better
]
425 # Show result of Hough transform in output image
426 #a_output[:,:,:] = hough[:,:,None]*20
428 # Find and classify peaks in Hough transform image
431 for y
in xrange(height
):
432 for x
in xrange(width
):
433 if hough
[y
, x
] > thresh
:
434 candidates
.append((hough
[y
,x
], y
, x
))
435 candidates
.sort(reverse
= True)
441 file_coords
= open(filename_coords
,'wt')
442 file_coords
.write('Autocount on %dx%d image\n' % (width
,height
))
444 suppress
= numpy
.zeros((height
,width
), 'bool')
445 claim_strength
= numpy
.zeros((height
,width
), 'float64') + 1e30
449 for _
,y
,x
in candidates
:
450 radius
= hough_radius
[y
,x
]
452 # Make sure candidate is not near a previously detected cell,
453 # and does not abut the edge of the image
454 if not suppress
[y
,x
] and \
455 y
>= radius
and x
>= radius
and y
< height
-radius
and x
< width
-radius
:
456 coords
= clipped_coordinates(0, radius
* Cell_inside_radius
, y
, x
, height
, width
)
458 #Only include foreground
459 this_fg
= fg
[coords
[:,0], coords
[:,1]] # [GB] Booleans
460 coords
= coords
[this_fg
]
462 this_strength
= ((coords
[:,0]-y
)**2 + (coords
[:,1]-x
)**2)
463 claim_strength
[coords
[:,0],coords
[:,1]] = numpy
.minimum(
464 claim_strength
[coords
[:,0],coords
[:,1]],
468 cells
.append((y
, x
, coords
, this_strength
))
471 # Suppress detection of cells too close to here
472 # (note: we do this irrespective of whether we treated this location as a hit)
473 coords
= clipped_coordinates(0,radius
* Cell_suppression_radius
, y
,x
,height
,width
)
474 suppress
[coords
[:,0],coords
[:,1]] = True
476 #for oy, ox in offsets(0, radius*Cell_suppression_radius): #suppression_offsets:
479 # if ny >= 0 and nx >= 0 and ny < height and nx < width:
480 # suppress[ny,nx] = True
482 for y
,x
, coords
, this_strength
in cells
:
483 radius
= hough_radius
[y
,x
]
485 if not len(coords
): continue # No fg pixels within radius
487 winners
= claim_strength
[coords
[:,0],coords
[:,1]] >= this_strength
488 coords
= coords
[winners
]
490 if not len(coords
): continue # Did not win any foreground pixels
494 dark_stain_pixels
= 0
502 for ny
, nx
in coords
:
512 if dark_stain
[ny
,nx
]:
513 dark_stain_pixels
+= 1
516 mass_y
= numpy
.average(coords
[:,0])
517 mass_x
= numpy
.average(coords
[:,1])
518 mass_offset
= numpy
.sqrt( (mass_y
-y
)**2 + (mass_x
-x
)**2 )
520 is_infected
= (infected_pixels
>= Infection_pixels_required
and stain_pixels
>= 1)
524 stain_x2
/stain_pixels
- (stain_x
/stain_pixels
)**2 +
525 stain_y2
/stain_pixels
- (stain_y
/stain_pixels
)**2
527 is_infected
= stain_spread
>= Stain_spread_required
530 (radius
> Max_cell_radius
) or
531 (radius
< Min_cell_radius
) or
532 (mass_offset
> Max_offcenter
* radius
) or
533 (dark_stain_pixels
>= fg_pixels
*Max_dark_stain_proportion
)
538 file_coords
.write('%d,%d,Rejected\n' % (x
,y
))
543 file_coords
.write('%d,%d,Infected\n' % (x
,y
))
547 file_coords
.write('%d,%d,Uninfected\n' % (x
,y
))
550 a_output
[coords
[:,0],coords
[:,1]] += [[ random
.random()*64+32 for i
in (0,1,2) ]]
552 # Mark detected cell on output image
553 if Show_computer_marks
:
554 set(a_output
, y
, x
, color
)
555 for i
in xrange(1,3):
556 set(a_output
,y
-i
,x
,color
)
557 set(a_output
,y
+i
,x
,color
)
558 set(a_output
,y
,x
-i
,color
)
559 set(a_output
,y
,x
+i
,color
)
566 # Save diagnostic output image
567 a_output
= numpy
.clip(a_output
,0,255)
568 i2
= Image
.fromarray(a_output
.astype('uint8'))
569 i2
.save(filename_out
)
573 def read_ac(filename
):
574 f
= open(filename
,'rU')
579 parts
= line
.rstrip().split(',')
580 if parts
[2] == 'Uninfected':
582 elif parts
[2] == 'Infected':
586 return n_total
, n_infected
589 if __name__
== '__main__':
593 name
, value
= args
.pop(0).split('=')
595 assert hasattr(sys
.modules
[__name__
], name
), 'Parameter named %s does not exist' % name
596 setattr(sys
.modules
[__name__
], name
, value
)
597 print name
, '=', value
600 print Help
% sys
.argv
[0]
603 if os
.path
.isfile(args
[0]):
604 in_dir
, filename
= os
.path
.split(args
[0])
605 filenames
= [ filename
]
610 filenames
= os
.listdir(in_dir
)
612 is_subprocess
= False
614 out_dir
= os
.path
.normpath(in_dir
) + '_output'
615 if not os
.path
.exists(out_dir
):
618 if not is_subprocess
:
619 for filename
in filenames
:
620 root
, ext
= os
.path
.splitext(filename
)
621 ac_filename
= os
.path
.join(out_dir
, root
+'.ac')
622 if os
.path
.exists(ac_filename
): os
.unlink(ac_filename
)
624 coord_filenames
= [ ]
628 for filename
in filenames
:
629 root
, ext
= os
.path
.splitext(filename
)
630 if ext
.lower() not in ('.bmp','.tif','.tiff','.png'): continue
632 filename_in
= os
.path
.join(in_dir
, filename
)
633 filename_out
= os
.path
.join(out_dir
, root
+'.tif')
634 filename_coords
= os
.path
.join(out_dir
, root
+'.ac')
636 if is_subprocess
or N_processors
== 1:
637 process(filename_in
, filename_out
, filename_coords
, root
)
639 while n_running
>= N_processors
:
642 os
.spawnl(os
.P_NOWAIT
, sys
.executable
, sys
.executable
, *(sys
.argv
[:-1]+[filename_in
]))
645 coord_filenames
.append((filename_coords
, root
))
651 if not is_subprocess
:
652 result_file
= open(os
.path
.join(out_dir
,'results.txt'),'wt')
654 for filename_coords
, root
in coord_filenames
:
655 n_total
, n_infected
= read_ac(filename_coords
)
656 print >> result_file
, root
, n_total
, n_infected