MATH590: C2examples, clustering of leukocytes

In this course module, clustering methods are used to classify blood sample images. We investigate whether clustering can recover the historical observations of blood samples that have led to classification of the blood cells revealed by staining and visual analysis of the morphology.

Blood physiology background

A prototypical problem in the biomedical field is high-throughput characterization of tissues. Samples from many subjects are shipped to a processing facility with specialized equipment, and automation of classification processes allows considerable benefits, both economic and for quality of care. 

The problem considered in this module is that of classification of leukocytes, or white blood cells (WBC), a type of blood cell produced in the bone marrow, and also found in lymph tissue. Leukocytes (from the Greek, Leuko = white, Cyte = urn, enclosed container, hence cell), are part of the immune system. Basic characterization of a blood sample includes counts of the number of leukocytes per microliter (the WBC number). Leukocytes are readily identified from a blood smear through the Wright or combined Wright-Giemsa stains. The Wright-Giemsa stain (introduced in 1902), is a combination of a red stain (eosin, Greek Eos = ancient Greek goddess of the dawn) that stains the cytoplasm pink and Methylene blue that stains cell nuclei. Red blood cells (RBC, eryhrocytes, Greek Erythros = red)  contain hemoglobin to carry out oxygen transport and do not have a nucleus. Since WBC do have a nucleus, they can readily be distinguished in a Wright-Giemsa stained blood smear.

ChemicalData["Eosin"]

Further characterization determines the type of observed leukocytes. There are 5 major types :
1. Basophils
2. Eosinophils
3. Lymphocytes
4. Monocytes
5. Neutrophils

Leukemia is a group of cancers of the bone marrow that results in large numbers of abnormal leukocytes, called blasts (in biology, a blast or precursor cell is only partially differentiated and has lost most stem cell properties). There are two types of leukemia:
1. In acute leukemia there is an overproduction of immature blood cells, both white and red that cannot maintain oxygenation and are malignant (i.e., can induce metastasis, the formation of cancers in additional tissues).
2. In chronic leukemia many abnormal leukocytes are produced that leads to general loss of circulatory function including tiredness, fever, risk of infection.

The above observations concisely synthesize decades of medical research into blood composition, cell morphology and pathology. We are interested in whether mathematical analysis of blood image data can reproduce the above body of knowledge.

Blood image database

Read data and carry out visual examination of images

Choose directory of images to process

SetDirectory["/home/student/courses/MATH590/SETdata/BloodCells/JPEGImages"]

List the files in the current directory

f=FileNames[]; Nim=Length[f]

Manipulate[Import[f[[i]]],{i,1,Nim,1}]

Here are two typical images that differentiate between WBC and RBC by staining (WBC appear blue), and also show different morphology of WBCs.

{im1=Import[f[[1]]],im2=Import[f[[9]]]}

The image data can be separated into individual color channels, adjusted to cover the entire intensity scale (contrast enhancement)

{imR,imG,imB}=Map[ImageAdjust[#]&,ColorSeparate[im1]]

It can also be reduced to a grayscale image, again adjusted to cover the entire intensity scale

imBW = ImageAdjust[ColorNegate[ColorConvert[im1,"Grayscale"]]]

Blood cell classification

A first question for mathematical data analysis is whether clustering techniques might reproduce the success of Wright staining. This is investigated by removing the blue color information (i.e., considering the situation before Wright staining) and processing the color image data.

Undoing the effect of Wright-Giemsa staining

To remove the Wright-Giemsa staining in the image data base, use a high-intensity threshold on the blue channel to identify the stained region, and smooth the resulting mask to account for experimental noise

imWrightGiemsa=Binarize[GaussianFilter[FillingTransform[Binarize[imB,0.7]],4]]

Mask each color channel to remove the Wright-Giemsa stain

{imRremovestain,imGremovestain,imBremovestain}=Map[ImageMultiply[#,ColorNegate[imWrightGiemsa]]&,{imR,imG,imB}]

A gray level representative of the cytosol background is placed in the masked region. To find the appropriate gray level, mask the average light absorption by  the RBC

imCyto=Binarize[imBremovestain,0.25]

These are the pixels within the cytosol

{imRcyto,imGcyto,imBcyto}={ImageMultiply[imR,imCyto],ImageMultiply[imG,imCyto],ImageMultiply[imB,imCyto]}

Define a function to create a list of non-zero pixel values

SelectCyto[im_,mask_]:=Select[Flatten[ImageData[ImageMultiply[im,mask]]],#>0.&];

Create lists of values that correspond to cytosol pixels

{imDataR,imDataG,imDataB}= Map[SelectCyto[#,imCyto]&, {imR,imG,imB}];

Find the average and standard deviation of the gray value in the cytosol pixels in each color channel

{avGrayR,sdGrayR}={Mean[imDataR],StandardDeviation[imDataR]}

{avGrayG,sdGrayG}={Mean[imDataG],StandardDeviation[imDataG]}

{avGrayB,sdGrayB}={Mean[imDataB],StandardDeviation[imDataB]}

Create a speckle gray pattern using a Gaussian distribution, and smooth out the variation

{mx,my}=ImageDimensions[imB]; nPixels=mx*my;
bR=Image[ArrayReshape[Table[RandomVariate[NormalDistribution[avGrayR,sdGrayR]],nPixels],{my,mx}]];
bG=Image[ArrayReshape[Table[RandomVariate[NormalDistribution[avGrayG,sdGrayG]],nPixels],{my,mx}]];
bB=Image[ArrayReshape[Table[RandomVariate[NormalDistribution[avGrayB,sdGrayB]],nPixels],{my,mx}]];
{backgroundR,backgroundG,backgroundB}=Map[GaussianFilter[#,4]&,{bR,bG,bB}]

Place the speckle gray pattern in the stained region in each color channel

imRnostain=ImageAdd[imRremovestain, ImageMultiply[imWrightGiemsa,backgroundR]];
imGnostain=ImageAdd[imGremovestain, ImageMultiply[imWrightGiemsa,backgroundG]];
imBnostain=ImageAdd[imBremovestain, ImageMultiply[imWrightGiemsa,backgroundB]];
{imRnostain,imGnostain,imBnostain}

A color image without the blue stain can now be recreated

im =ColorCombine[{imRnostain,imGnostain,imBnostain},"RGB"]

Note that the threshold values in the above were arbitrary, and a considerable number of steps were needed to undo the effect of the Wright-Giemsa staining. Combine all these steps into a single function that we can then Manipulate to determine appropriate threshold levels. The following function gathers all the individual steps taken above. Two modifications are made:
1) The individual color channels are not contrast-adjusted individually since this would modify the color balance in the processed image;
2) An overall blue brightness value for the image used in the testing is computed and used to scale the blue level in other images

im1=Import[f[[1]]];
im1B = ColorSeparate[im1][[3]]; 
im1Blevel =Mean[Flatten[ImageData[im1B]]];

SelectCyto[im_,mask_]:=Select[Flatten[ImageData[ImageMultiply[im,mask]]],#>0.&];
RemoveWrightGiemsa[im_,tB_]:=Module[{mx,my,nPixels,imR,imG,imB,imBdata,imBlevel,t,imWrightGiemsa,imRremovestain,imGremovestain,imBremovestain,imCyto,imRcyto,imGcyto,imBcyto,imDataR,imDataG,imDataB,avGrayR,sdGrayR,avGrayG,sdGrayG,avGrayB,sdGrayB,bR,bG,bB,backgroundR,backgroundG,backgroundB,imRnostain,imGnostain,imBnostain,imNoStain },
{mx,my}=ImageDimensions[im]; nPixels=mx*my;
{imR,imG,imB}=ColorSeparate[im];
imBdata=ImageData[imB];
imBlevel =Mean[Flatten[imBdata]];
imB = Image[imBdata*im1Blevel/imBlevel];
t =tB*im1Blevel/imBlevel;
imWrightGiemsa=Binarize[GaussianFilter[FillingTransform[Binarize[imB,t]],4]];
{imRremovestain,imGremovestain,imBremovestain}=Map[ImageMultiply[#,ColorNegate[imWrightGiemsa]]&,{imR,imG,imB}];
imCyto=Binarize[imBremovestain,0.5t];
{imRcyto,imGcyto,imBcyto}={ImageMultiply[imR,imCyto],ImageMultiply[imG,imCyto],ImageMultiply[imB,imCyto]};
{imDataR,imDataG,imDataB}= Map[SelectCyto[#,imCyto]&, {imR,imG,imB}];
{avGrayR,sdGrayR}={Mean[imDataR],StandardDeviation[imDataR]};
{avGrayG,sdGrayG}={Mean[imDataG],StandardDeviation[imDataG]};
{avGrayB,sdGrayB}={Mean[imDataB],StandardDeviation[imDataB]};bR=Image[ArrayReshape[Table[RandomVariate[NormalDistribution[avGrayR,sdGrayR]],nPixels],{my,mx}]];
bG=Image[ArrayReshape[Table[RandomVariate[NormalDistribution[avGrayG,sdGrayG]],nPixels],{my,mx}]];
bB=Image[ArrayReshape[Table[RandomVariate[NormalDistribution[avGrayB,sdGrayB]],nPixels],{my,mx}]];
{backgroundR,backgroundG,backgroundB}=Map[GaussianFilter[#,8]&,{bR,bG,bB}];
imRnostain=ImageAdd[imRremovestain, ImageMultiply[imWrightGiemsa,backgroundR]];
imGnostain=ImageAdd[imGremovestain, ImageMultiply[imWrightGiemsa,backgroundG]];
imBnostain=ImageAdd[imBremovestain, ImageMultiply[imWrightGiemsa,backgroundB]];
{imRnostain,imGnostain,imBnostain};
imNoStain =ColorCombine[{imRnostain,imGnostain,imBnostain},"RGB"];
Return[imNoStain ]
];

With the above encapsulation of all the image processing steps needed to remove the effect of the Wright-Giemsa stain, the appropriate threshold values can be determined interactively

Manipulate[RemoveWrightGiemsa[im1,tB],{{tB,0.81},0.75,0.9}]

For image im1 the appropriate threshold value that removes most of the blue stain but maintains an observable separation of the leukocyte to the surrounding fluid seems to be tB=0.81. Check this value on another image

Manipulate[RemoveWrightGiemsa[im2,tB],{{tB,0.81},0.75,0.9}]

Though not all the blue stain is eliminated, the threshold seems effective, and will be used henceforth;

{im1NoB,im2NoB}={RemoveWrightGiemsa[im1,0.81],RemoveWrightGiemsa[im2,0.81]}

Background separation

Define a function that interactively investigates the masking of the image to isolate the blood cells from the background. Again, this is first presented as a sequence of individual image processing steps, and then encapsulated into a single function to be used in subsequent analyses.

First, generate gray-scale images for testing

{im1BW,im2BW} = Map[ImageAdjust[ColorNegate[ColorConvert[#,"Grayscale"]]]&,{im1NoB,im2NoB}]

Interactively determine the appropriate threshold to obtain masking of the cytosol background

Manipulate[ImageMultiply[Binarize[FillingTransform[im1BW],t],im1BW],{{t,0.44},0.4,0.6}]

Adopt a threshold of 0.44 and define a mask to isolate the blood cells

{mask1,mask2} =Map[ Binarize[FillingTransform[#],.44]&,{im1BW,im2BW}]

The masks exhibit numerous small regions related to noisy data. Assume that at most MaxCells could be present in the image, and select only those groups with enough pixels to contain at least one cell. The steps to do this are:

1. Compute the gradient of the mask to identify edges

{Gmask1,Gmask2}=Map[Binarize[GradientFilter[#,1]]&,{mask1,mask2}]

2. Fill in closed regions

{FGmask1,FGmask2}=Map[FillingTransform[#]&,{Gmask1,Gmask2}]

3. Assume the image can contain at most MaxCells, and compute from this the minimum number of pixels needed to represent a cell

MaxCells = 100; MinPixels =Floor[Apply[Times,ImageDimensions[mask1]]/MaxCells];

4. Define a function that provides a mask of the image noise

ImageNoiseMask[im_]:=ColorNegate[Binarize[ColorConvert[Colorize[SelectComponents[im, #Count< MinPixels&]],"GrayScale"]]]

{Nmask1,Nmask2}=Map[ImageNoiseMask[#]&,{FGmask1,FGmask2}]

5. Update the cell background mask

{mask1,mask2}={ImageMultiply[mask1,Nmask1],ImageMultiply[mask2,Nmask2]}

6. There still exist numerous small, convoluted regions that cannot correspond to cells. To eliminate these consider the concept of an image skeleton, i.e., the lines that are equidistant to the boundary

Smask1=SkeletonTransform[mask1]

Smask2=SkeletonTransform[mask2]

7. Now, prune the skeleton by removing outermost branches of the tree structure

minLength=Floor[Sqrt[MinPixels]/2];
PSmask1=Pruning[Smask1,minLength]

PSmask2=Pruning[Smask2,minLength]

8. Reconstruct a smoothed mask shape

{SPSmask1,SPSmask2}=Map[InverseDistanceTransform[#]&,{PSmask1,PSmask2}]

The cells are the result of applying the mask to the B&W image

{cells1,cells2}={ImageMultiply[SPSmask1,im1NoB],ImageMultiply[SPSmask2,im2NoB]};

cells1

cells2

Notice that some cells can be isolated, e.g., for cells2:

MorphologicalComponents[cells2] // Colorize

However, consistent isolation of individual cells is not possible due to overlapping cells. This is an indication of the limitations of methods based solely on clustering of data values. In addition to data value clustering, topological analysis of the distribution of pixel brightness values values will be necessary. In particular, image data will need to be partitioned into convex subregions. This shall be carried out in the TOP module after introducing the necessary theoretical concepts.

Isolating leukocytes by Wright-Giemsa staining

The additional information provided by chemical staining and the relative paucity of leukocytes (i.e., WBC are well separated in the blood smear data) allows straightforward identification of white blood cells.

Repeat the isolation of the individual color channels, and adjust for maximum contrast.

{imR1,imG1,imB1}=Map[ImageAdjust[#]&,ColorSeparate[im1]]

The green channel immediately furnishes a mask to eliminate the background fluid

mask1=ColorNegate[Binarize[imG1]]

Eliminate the background from the blue channel data

cellsB1=ImageMultiply[mask1,imB1]

Generate a mask for the leukocytes

Lmask1=FillingTransform[Binarize[cellsB1,0.6]]

Filter out small components from the mask

MaxCells = 100; MinPixels =Floor[Apply[Times,ImageDimensions[Lmask1]]/MaxCells];

FLmask1=Binarize[ColorConvert[Colorize[SelectComponents[Lmask1,#Count>MinPixels&]],"Grayscale"]]

Define the smallest possible image that contains just leukocytes

wbc=ImageCrop[ImageMultiply[FLmask1,im1]]

Encapsulate all these operations into a single function

Leukocytes[im_,t_]:=Module[{imR,imG,imB,mask,MaxCells,MinPixels,cellsB,Lmask,FLmask,wbc},
{imR,imG,imB}=Map[ImageAdjust[#]&,ColorSeparate[im]];
mask=ColorNegate[Binarize[imG]];
MaxCells = 100; MinPixels =Floor[Apply[Times,ImageDimensions[mask]]/MaxCells];
cellsB=ImageMultiply[mask,imB];
Lmask=FillingTransform[Binarize[cellsB,t]];
FLmask=Binarize[ColorConvert[Colorize[SelectComponents[Lmask,#Count>MinPixels&]],"Grayscale"]];
wbc=ImageCrop[ImageMultiply[FLmask,im]];
Return[wbc];
]

Leukocytes[im2,0.6]

Manipulate[Leukocytes[Import[f[[k]]],t],{k,1,Length[f],1},{{t,0.6},0.5,0.7}]