detect and count cell nuclei in whole slide images -凯发k8网页登录

download sample image from camelyon17 data set

this example uses one wsi from the challenge. navigate to one of the online repositories specified on the page, then download the file patient_000_node_1.tif from the directory training/center_0/patient_000. update the filename variable to point to the full path of the downloaded image.

filename = "patient_000_node_1.tif";

create and preprocess blocked image

create a object from the sample image. this image has nine resolution levels. the finest resolution level, which is the first level, has a size of 197226-by-96651 pixels. the coarsest resolution level, which is the last level, has a size of 770-by-377 pixels and easily fits in memory.

bim = blockedimage(filename);
disp(bim.numlevels)

     9

disp(bim.size)

      197226       96651           3
       98613       48325           3
       49306       24162           3
       24653       12081           3
       12326        6040           3
        6163        3020           3
        3081        1510           3
        1540         755           3
         770         377           3

adjust spatial extents of each level

the tif image metadata contains information about the world coordinates. the exact format of the metadata varies across different sources. the file in this example stores the mapping between pixel subscripts and world coordinates using the xresolution, yresolution, and resolutionunit fields. the indicates that xresolution and yresolution are the number of pixels per resolutionunit in the corresponding directions.

fileinfo = imfinfo(bim.source);
info = table((1:9)',[fileinfo.xresolution]',[fileinfo.yresolution]',{fileinfo.resolutionunit}', ...
    variablenames=["level","xresolution","yresolution","resolutionunit"]);
disp(info)

    level    xresolution    yresolution    resolutionunit
    _____    ___________    ___________    ______________
      1         41136          41136       {'centimeter'}
      2         20568          20568       {'centimeter'}
      3         10284          10284       {'centimeter'}
      4        5142.1         5142.1       {'centimeter'}
      5          2571           2571       {'centimeter'}
      6        1285.5         1285.5       {'centimeter'}
      7        642.76         642.76       {'centimeter'}
      8        321.38         321.38       {'centimeter'}
      9        160.69         160.69       {'centimeter'}

calcuate the size of a single pixel, or the pixel extent, in world coordinates at the finest resolution level. the pixel extent is identical in the x and the y directions because the values of xresolution and yresolution are equal. then, calculate the corresponding full image size in world coordinates.

level = 1;
pixelextentinworld = 1/fileinfo(level).xresolution;
imagesizeinworld = bim.size(level,1:2)*pixelextentinworld;

adjust the spatial extents of each level such that the levels cover the same physical area in the real world. to ensure that subscripts map to pixel centers, offset the starting coordinate and the ending coordinate by a half pixel.

bim.worldstart = [pixelextentinworld/2 pixelextentinworld/2 0.5];
bim.worldend = bim.worldstart   [imagesizeinworld 3];

display the image using the function. bigimageshow displays the image in world coordinates, with each axis measured in centimeters.

bigimageshow(bim)
title("wsi of lymph node")

develop basic nuclei detection algorithm

this example develops the nuclei detection algorithm using a representative region of interest (roi), then applies the algorithm on all regions of the image that contain tissue.

identify and display a representative roi in the image.

xregion = [0.53 0.55];
yregion = [2.064 2.08];
xlim(xregion)
ylim(yregion)
title("region of interest")

convert the world coordinates of the roi to pixel subscripts. then, extract the roi data at the finest resolution level using the getregion function.

substart = world2sub(bim,[yregion(1) xregion(1)]);
subend   = world2sub(bim,[yregion(2) xregion(2)]);
imregion = getregion(bim,substart,subend,level=1);

detect nuclei using color thresholding

the nuclei have a distinct stain color that enables segmentation using color thresholding. this example provides a helper function, createnucleimask, that segments nuclei using color thresholding in the rgb color space. the thresholds have been selected, and the code for the helper function has been generated, using the colorthresholder app. the helper function is attached to the example as a supporting file.

if you want to segment the nuclei regions using a different color space or select different thresholds, use the colorthresholder app. open the app and load the image using this command:

colorthresholder(imregion)

after you select a set of thresholds in the color space of your choice, generate the segmentation function by selecting export from the toolstrip and then selecting export function. save the generated function as createnucleimask.m.

apply the segmentation to the roi, and display the mask over the roi as a falsecolor composite image.

nucleimask = createnucleimask(imregion);
imshowpair(imregion,nucleimask)

select nuclei using image region properties

explore properties of the regions in the mask using the imageregionanalyzer app. open the app by using this command:

imageregionanalyzer(nucleimask)

clean up regions with holes inside the nuclei by selecting fill holes on the toolstrip.

sort the regions by area. regions with the largest areas typically correspond to overlapping nuclei. becauses your goal is to segment single cells, you can ignore the largest regions at the expense of potentially underestimating the total number of nuclei.

scroll down the list until you find the largest region that corresponds to a single nucleus. specify maxarea as the area of this region, and specify equivdiameter as the equivalent diameter, represented by the equivdiameter property in the region properties pane. calculate the equivalent radius of the region as half the value of equivdiameter.

maxarea = 405;
equivdiameter = 22.7082;
equivradius = round(equivdiameter/2);

regions with very small area typically correspond to partial nuclei. continue scrolling until you find the smallest region that corresponds to a whole nucleus. set minarea as the area of this region. because your goal is to segment whole cells, you can ignore regions with areas smaller than minarea.

minarea = 229;

filter the image so it contains only the regions with areas between minarea and maxarea. from the toolstrip, select filter. the binary image updates in real time to reflect the current filtering parameters.

you can optionally change the range of areas or add additional region filtering constraints. for example, you can require that nuclei be approximately round by excluding regions with large perimeters. specify the filter parameters in the filter regions dialog box.

this example provides a helper function, calculatenucleiproperties, that calculates the area and centroid of regions within the range of areas [229, 405]. the code for the helper function has been generated using the imageregionanalyzer app. the helper function is attached to the example as a supporting file.

if you specify different filtering constraints, then you must generate a new filtering function by selecting export from the toolstrip and then selecting export function. save the generated function as calculatenucleiproperties.m. after you generate the function, you must open the function and edit the last line of code so that the function returns the region centroids. add 'centroid' to the cell array of properties measured by the regionprops function. for example:

properties = regionprops(bw_out,{'centroid','area','perimeter'});

calculate the centroids for each likely detection by using the calculatenucleiproperties helper function. overlay the centroids as yellow points over the roi.

[~,stats] = calculatenucleiproperties(nucleimask);
centroidcoords = vertcat(stats.centroid);
figure
imshow(imregion)
hold on
plot(centroidcoords(:,1),centroidcoords(:,2),"y.",markersize=10)
hold off

apply nuclei detection algorithm

to process wsi data efficiently at fine resolution levels, apply the nuclei detection algorithm on blocks consisting of tissue and ignore blocks without tissue. you can identify blocks with tissue using a coarse tissue mask. for more information, see .

create tissue mask

the tissue is brighter than the background, which enables segmentation using color thresholding. this example provides a helper function, createtissuemask, that segments tissue using color thresholding in the l*a*b* color space. the thresholds have been selected, and the code for the helper function has been generated, using the colorthresholder app. the helper function is attached to the example as a supporting file.

if you want to segment the tissue using a different color space or select different thresholds, use the colorthresholder app. open a low-resolution version of the image in the app using these commands:

imcoarse = gather(bim);
colorthresholder(imcoarse)

after you select a set of thresholds, generate the segmentation function by selecting export from the toolstrip and then selecting export function. save the generated function as createtissuemask.m.

create a mask of the regions of interest that contain tissue by using the function. in this code, the apply function runs the createtissuemask function on each block of the original blocked image bim and assembles the tissue masks into a single blockedimage with the correct world coordinates. this example creates the mask at the seventh resolution level, which fits in memory and follows the contours of the tissue better than a mask at the coarsest resolution level.

masklevel = 7;
tissuemask = apply(bim,@(bs)createtissuemask(bs.data),level=masklevel);

overlay the mask on the original wsi image using the bigimageshow function.

hbim = bigimageshow(bim);
showlabels(hbim,tissuemask,alphadata=tissuemask,alphamap=[0.5 0])
title("tissue mask; background has green tint")

select blocks with tissue

select the set of blocks on which to run the nuclei detection algorithm by using the selectblocklocations function. include only blocks that have at least one pixel inside the mask by specifying the inclusionthreshold name-value argument as 0. specify a large block size to perform the nuclei detection more efficiently. the upper limit of the block size depends on the amount of system ram. if you perform the detection using parallel processing, you may need to reduce the block size.

blocksize = [2048 2048 3];
bls = selectblocklocations(bim,blocksize=blocksize, ...
    masks=tissuemask,inclusionthreshold=0,levels=1);

apply algorithm on blocks with tissue

perform the nuclei detection on the selected blocks using the function. in this code, the apply function runs the countnuclei helper function (defined at the end of the example) on each block of the blocked image bim and returns a structure with information about the area and centroid of each nucleus in the block. to ensure the function detects and counts nuclei on the edges of a block only once, add a border size equal to the equivalent radius value of the largest nucleus region. to accelerate the processing time, process the blocks in parallel.

if(exist("bndetections","dir"))
    rmdir("bndetections","s")
end
bndetections = apply(bim,@(bs)countnuclei(bs), ...
    blocklocationset=bls,bordersize=[equivradius equivradius], ...
    adapter=images.blocked.matblocks,outputlocation="bndetections", ...
    useparallel=true);

starting parallel pool (parpool) using the 'processes' profile ...
connected to parallel pool with 6 workers.

display heatmap of detected nuclei

create a heatmap that shows the distribution of detected nuclei across the tissue sample. each pixel in the heatmap contains the count of detected nuclei in a square of size 0.01-by-0.01 square centimeters in world coordinates.

determine the number of pixels that this real world area contains:

numpixelsinsquare = world2sub(bim,[0.01 0.01]);

calculate the size of the heatmap, in pixels, at the finest level. then initialize a heatmap consisting of all zeros.

bimsize = bim.size(1,1:2);
heatmapsize = round(bimsize./numpixelsinsquare);
heatmap = zeros(heatmapsize);

load the blocked image of detected nuclei into memory. get the list of the (row, column) pixel indices of the centroids of all detected nuclei.

ndetections = gather(bndetections);
centroidspixelindices = vertcat(ndetections.ncentroidsrc);

map the nuclei centroid coordinates to pixels in the heatmap. first, normalize the (row, column) centroid indices relative to the size of the blocked image, in pixels. then, multiply the normalized coordinates by the size of the heatmap to yield the (row, column) pixel indices of the centroids in the heatmap.

centroidsrelative = centroidspixelindices./bimsize;
centroidsinheatmap = ceil(centroidsrelative.*heatmapsize);

loop through the nuclei and, for each nucleus, increment the value of the corresponding heatmap pixel by one.

for ind = 1:size(centroidsinheatmap)
    r = centroidsinheatmap(ind,1);
    c = centroidsinheatmap(ind,2);
    heatmap(r,c) = heatmap(r,c) 1;
end

display the heatmap.

maxcount = max(heatmap,[],"all");
figure
imshow(heatmap,hot(maxcount))
hcb = colorbar;
hcb.label.string = "nuclei density (count/0.01 cm square)";

convert the heatmap into a blocked image with the correct world coordinates.

bheatmap = blockedimage(heatmap, ...
    worldstart=bim.worldstart(1,1:2),worldend=bim.worldend(1,1:2));

display the original wsi data next to the heatmap.

tiledlayout(1,2)
nexttile
hl = bigimageshow(bim);
title("wsi image")
nexttile
hr = bigimageshow(bheatmap,cdatamapping="direct");
colormap(hr.parent,hot(maxcount))
title("nuclei count heatmap")

link the axes of the two images to enable comparison.

linkaxes([hl.parent hr.parent])

zoom in on an interesting region.

xlim(hl.parent,[0.56514 0.85377])
ylim(hl.parent,[2.6624 3.2514])

inspect the detected nuclei across the entire wsi using the helper function explorenucleidetectionresults, which is attached to the example as a supporting file. the helper function displays the wsi with a draggable overview rectangle (shown in blue) next to a detail view of the portion of the image within the overview rectangle. you can move and resize the overview window interactively, and the detail view updates to display the new selected region and its centroids.

roi = [0.62333 2.0725 0.005 0.005]; 
h = explorenucleidetectionresults(bim,bndetections,roi);

supporting function

the countnuclei helper function performs these operations to count nuclei within a block of blocked image data:

to implement a more sophisticated nuclei detection algorithm, replace the calls to the createnucleimask and calculatenucleiproperties functions with calls to your custom functions. preserve the rest of the countnuclei helper function for use with the apply function.

function nucleiinfo = countnuclei(bs)
im = bs.data;
% create a binary mask of regions where nuclei are likely to be present
mask = createnucleimask(im);
% calculate the area and centroid of likely nuclei
[~,stats] = calculatenucleiproperties(mask);
% initialize result structure 
nucleiinfo.topleftrc = bs.start(1:2);
nucleiinfo.botrightrc = bs.end(1:2);
nucleiinfo.narea = [];
nucleiinfo.ncentroidsrc = [];
if ~isempty(stats)
    % go from (x,y) coordinates to (row, column) pixel subscripts
    % relative to the block
    nucleiinfo.ncentroidsrc = fliplr(vertcat(stats.centroid));
    nucleiinfo.narea = vertcat(stats.area);
    % trim out centroids on the borders to prevent double-counting
    onborder = nucleiinfo.ncentroidsrc(:,1) < bs.bordersize(1) ...
        | nucleiinfo.ncentroidsrc(:,1) > size(im,1)-2*bs.bordersize(1) ...
        | nucleiinfo.ncentroidsrc(:,2) < bs.bordersize(2) ...
        | nucleiinfo.ncentroidsrc(:,2) > size(im,2)-2*bs.bordersize(2);
    nucleiinfo.ncentroidsrc(onborder,:) = [];
    nucleiinfo.narea(onborder) = [];
    % convert the coordinates of block centroids from (row, column)
    % subscripts relative to the block to (row, column) subscripts relative
    % to the global image
    nucleiinfo.ncentroidsrc = nucleiinfo.ncentroidsrc   nucleiinfo.topleftrc;
end
end