Introduction

The only definitive way to diagnose breast cancer is by examining tissue samples collected from biopsy or surgery. The samples are commonly prepared with hematoxylin and eosin (H&E) staining to increase the contrast of structures in the tissue. Traditionally, pathologists examine the tissue on glass slides under a microscope to detect tumor tissue. Diagnosis takes time as pathologists must thoroughly inspect an entire slide at close magnification. Further, pathologists may not notice small tumors. Deep learning methods aim to automate the detection of tumor tissue, saving time and improving the detection rate of small tumors.

Deep learning methods for tumor classification rely on digital pathology, in which whole tissue slides are imaged and digitized. The resulting WSIs have extremely high resolution. WSIs are frequently stored in a multiresolution file to facilitate the display, navigation, and processing of the images.

Reading WSIs is a challenge because the images cannot be loaded as a whole into memory and therefore require out-of-core image processing techniques. Image Processing Toolbox™ bigimage (Image Processing Toolbox) objects can store and process this type of large multiresolution image. bigimageDatastore (Image Processing Toolbox) objects can prepare training patches from a bigimage to feed into a neural network.

This example shows how to train a deep learning network to classify tumors in very large multiresolution images using bigimage and bigimageDatastore. The example presents classification results as heatmaps that depict the probability that local tissue is tumorous. The localization of tumor regions enables medical pathologists to investigate specific regions and quickly identify tumors of any size in an image.

Download Training Data

This example uses WSIs from the Camelyon16 challenge [1]. The data from this challenge contains a total of 400 WSIs of lymph nodes from two independent sources, separated into 270 training images and 130 test images. The WSIs are stored as TIF files in a stripped format with an 11-level pyramid structure.

bigimage objects warn when images contain stripped formats because this format is slow to process. To avoid displaying warnings in the command window, turn off the warning while the example runs and restore the warning state on completion.

warningState = warning('off','images:bigimage:singleStripTiff');
c = onCleanup(@()warning(warningState));

The training data set consists of 159 WSIs of normal lymph nodes and 111 WSIs of lymph nodes with tumor and healthy tissue. Usually, the tumor tissue is a small fraction of the healthy tissue. Ground truth coordinates of the lesion boundaries accompany the tumor images.

The size of each training file is approximately 2 GB. If you do not want to download the training data set or train the network, then go directly to the Download and Preprocess Testing Data section in this example.

Create a directory to store the training data set.

trainingImageDir = fullfile(tempdir,'Camelyon16','training');
if ~exist(trainingImageDir,'dir')
    mkdir(trainingImageDir);
    mkdir(fullfile(trainingImageDir,'normal'));
    mkdir(fullfile(trainingImageDir,'tumor'));
    mkdir(fullfile(trainingImageDir,'lesion_annotations'));
end

trainNormalDataDir = fullfile(trainingImageDir,'normal');
trainTumorDataDir = fullfile(trainingImageDir,'tumor');
trainTumorAnnotationDir = fullfile(trainingImageDir,'lesion_annotations');

To download the training data, go to the Camelyon17 website and click the first "CAMELYON16 data set" link. Open the "training" directory, then follow these steps.

Specify the number of training images. Note that one of the training images of normal tissue, 'normal_144.tif', has metadata that cannot be read by the bigimage object. This example uses the remaining 158 training files.

numNormalFiles = 158;
numTumorFiles = 111; 

Visualize Training Data

To get a better understanding of the training data, display one training image. Because loading the entire image into memory at the finest resolution is not possible, you cannot use traditional image display functions such as imshow. To display and process the image data, use a bigimage (Image Processing Toolbox) object.

Create a bigimage object from a tumor training image.

tumorFileName = fullfile(trainTumorDataDir,'tumor_001.tif');
tumorImage = bigimage(tumorFileName);

Inspect the dimensions of the bigimage at each resolution level. Level 1 has the most pixels and is the finest resolution level. Level 10 has the fewest pixels and is the coarsest resolution level. The aspect ratio is not consistent, which indicates that levels do not all span the same world area.

levelSizeInfo = table((1:length(tumorImage.LevelSizes))', ...
    tumorImage.LevelSizes(:,1), ...
    tumorImage.LevelSizes(:,2), ...
    tumorImage.LevelSizes(:,1)./tumorImage.LevelSizes(:,2), ...
    'VariableNames',["Resolution Level" "Image Width" "Image Height" "Aspect Ratio"])

levelSizeInfo=11×4 table
    Resolution Level    Image Width    Image Height    Aspect Ratio
    ________________    ___________    ____________    ____________

            1           2.2118e+05        97792           2.2618   
            2           1.1059e+05        49152             2.25   
            3                55296        24576             2.25   
            4                27648        12288             2.25   
            5                13824         6144             2.25   
            6                 7168         3072           2.3333   
            7                 3584         1536           2.3333   
            8                 2048         1024                2   
            9                 1024          512                2   
           10                  512          512                1   
           11                 1577         3629          0.43455   

Display the bigimage at a coarse resolution level by using the bigimageshow (Image Processing Toolbox) function. Return a handle to the bigimageshow object. You can use the handle to adjust the display. The image contains a lot of empty white space. The tissue occupies only a small portion of the image.

h = bigimageshow(tumorImage,'ResolutionLevel',7);

Zoom in on one part of the image by setting the horizontal and vertical spatial extents with respect to the finest resolution level. The image looks blurry because this resolution level is very coarse.

xlim([29471,29763]);
ylim([117450,118110]);

To see more details, change the resolution level to a finer level.

h.ResolutionLevel = 1;

Create Masks

You can reduce the amount of computation by processing only regions of interest (ROIs). Use a mask to define ROIs. A mask is a logical image in which true pixels represent the ROI.

To further reduce the amount of computation, create masks at a coarse resolution level that can be processed entirely in memory instead of on a block-by-block basis. If the spatial referencing of the coarse resolution level matches the spatial referencing of finer resolution levels, then locations at the coarse level correspond to locations in finer levels. In this case, you can use the coarse mask to select which blocks to process at the finer levels. For more information, see Set Spatial Referencing for Big Images (Image Processing Toolbox) and Process Big Images Efficiently Using Mask (Image Processing Toolbox).

Specify a resolution level to use for creating the mask. This example uses resolution level 7, which is coarse and fits in memory.

resolutionLevel = 7;

Create Masks for Normal Images

In normal images, the ROI consists of healthy tissue. The color of healthy tissue is distinct from the color of the background, so use color thresholding to segment the image and create an ROI. The L*a*b* color space provides the best color separation for segmentation. Convert the image to the L*a*b* color space, then threshold the a* channel to create the tissue mask.

You can use the helper function createMaskForNormalTissue to create masks using color thresholding. This helper function is attached to the example as a supporting file.

The helper function performs these operations on each training image of normal tissue:

trainNormalMaskDir = fullfile(trainNormalDataDir,['normal_mask_level' num2str(resolutionLevel)]);
createMaskForNormalTissue(trainNormalDataDir,trainNormalMaskDir,resolutionLevel)

Now that both normal images and masks are on disk, create bigimage objects to manage the data by using the helper function createBigImageAndMaskArrays. This function creates an array of bigimage objects from the normal images and a corresponding array of bigimage objects from the normal mask images. The helper function is attached to the example as a supporting file.

[bigNormalImages,bigNormalMasks] = createBigImageAndMaskArrays(trainNormalDataDir,trainNormalMaskDir);

Select a sample normal image and mask. Confirm that the spatial world extents of the mask match the extents of the image at the finest resolution level. The spatial world extents are specified by the XWorldLimits and YWorldLimits properties.

idx = 2;
bigNormalImages(idx).SpatialReferencing(1)

ans = 
  imref2d with properties:

           XWorldLimits: [0 97792]
           YWorldLimits: [0 221184]
              ImageSize: [221184 97792]
    PixelExtentInWorldX: 1
    PixelExtentInWorldY: 1
    ImageExtentInWorldX: 97792
    ImageExtentInWorldY: 221184
       XIntrinsicLimits: [0.5000 9.7793e+04]
       YIntrinsicLimits: [0.5000 2.2118e+05]

bigNormalMasks(idx).SpatialReferencing(1)

ans = 
  imref2d with properties:

           XWorldLimits: [0 97792]
           YWorldLimits: [0 221184]
              ImageSize: [3456 1527]
    PixelExtentInWorldX: 64.0419
    PixelExtentInWorldY: 64
    ImageExtentInWorldX: 97792
    ImageExtentInWorldY: 221184
       XIntrinsicLimits: [0.5000 1.5275e+03]
       YIntrinsicLimits: [0.5000 3.4565e+03]

Verify that the mask contains the correct ROIs and spatial referencing. Display the sample image by using the bigimageshow function. Get the axes containing the display.

figure
hBigNormal = bigimageshow(bigNormalImages(idx));
hNormalAxes = hBigNormal.Parent;

Create a new axes on top of the displayed bigimage. In the new axes, display the corresponding mask image with partial transparency. The mask highlights the regions containing normal tissue.

hMaskAxes = axes;
hBigMask = bigimageshow(bigNormalMasks(idx),'Parent',hMaskAxes, ...
    'Interpolation','nearest','AlphaData',0.5);
hMaskAxes.Visible = 'off';

Link the axes of the image with the axes of the mask. When you zoom and pan, both axes are updated identically.

linkaxes([hNormalAxes,hMaskAxes]);

Zoom in on one part of the image by setting the horizontal and vertical spatial extents. The mask overlaps the normal tissue correctly.

xlim([45000 80000]);
ylim([130000 165000]);

Create Masks for Tumor Images

In tumor images, the ROI consists of tumor tissue. The color of tumor tissue is similar to the color of healthy tissue, so you cannot use color segmentation techniques. Instead, create ROIs by using the ground truth coordinates of the lesion boundaries that accompany the tumor images.

You can use the helper function createMaskForTumorTissue to create a mask using ROIs. This helper function is attached to the example as a supporting file.

The helper function performs these operations on each training image of tumor tissue:

trainTumorMaskDir = fullfile(trainTumorDataDir,['tumor_mask_level' num2str(resolutionLevel)]);
createMaskForTumorTissue(trainTumorDataDir,trainTumorAnnotationDir, ...
    trainTumorMaskDir,resolutionLevel);

Now that both tumor images and masks are on disk, create bigimage objects to manage the data by using the helper function createBigImageAndMaskArrays. This function creates an array of bigimage objects from the tumor images and a corresponding array of bigimage objects from the tumor mask images. The helper function is attached to the example as a supporting file.

[bigTumorImages,bigTumorMasks] = createBigImageAndMaskArrays(trainTumorDataDir,trainTumorMaskDir);

Select a sample tumor image and mask. Confirm that the spatial world extents of the mask match the extents of the image at the finest resolution level. The spatial world extents are specified by the XWorldLimits and YWorldLimits properties.

idx = 5;
bigTumorImages(idx).SpatialReferencing(1)

ans = 
  imref2d with properties:

           XWorldLimits: [0 97792]
           YWorldLimits: [0 219648]
              ImageSize: [219648 97792]
    PixelExtentInWorldX: 1
    PixelExtentInWorldY: 1
    ImageExtentInWorldX: 97792
    ImageExtentInWorldY: 219648
       XIntrinsicLimits: [0.5000 9.7793e+04]
       YIntrinsicLimits: [0.5000 2.1965e+05]

bigTumorMasks(idx).SpatialReferencing(1)

ans = 
  imref2d with properties:

           XWorldLimits: [0 97792]
           YWorldLimits: [0 219648]
              ImageSize: [3432 1527]
    PixelExtentInWorldX: 64.0419
    PixelExtentInWorldY: 64
    ImageExtentInWorldX: 97792
    ImageExtentInWorldY: 219648
       XIntrinsicLimits: [0.5000 1.5275e+03]
       YIntrinsicLimits: [0.5000 3.4325e+03]

Verify that the mask contains the correct ROIs and spatial referencing. Display the sample image by using the bigimageshow function. Get the axes containing the display.

figure
hBigTumor = bigimageshow(bigTumorImages(idx));
hTumorAxes = hBigTumor.Parent;

Create a new axes on top of the displayed big image. In the new axes, display the corresponding mask image with partial transparency. The mask highlights the regions containing normal tissue.

hMaskAxes = axes;
hBigMask = bigimageshow(bigTumorMasks(idx),'Parent',hMaskAxes, ...
    'Interpolation','nearest','AlphaData',0.5);
hMaskAxes.Visible = 'off';

Link the axes of the image with the axes of the mask. When you zoom and pan, both axes are updated identically.

linkaxes([hTumorAxes,hMaskAxes]);

Zoom in on one part of the image by setting the horizontal and vertical spatial extents. The mask overlaps the tumor tissue correctly.

xlim([45000 65000]);
ylim([130000 150000]);

Create Big Image Datastores for Training and Validation

To extract patches of training data from the bigimage objects, use a bigimageDatastore (Image Processing Toolbox). This datastore reads patches of bigimage data at a single resolution level.

Color imbalance and class imbalance in the raw training patches can potentially bias the network. Color imbalance results from nonuniform color staining of the tissue. Class imbalance results from an unequal amount of tumor and normal tissue in the data. To correct these imbalances, you can preprocess and augment the datastore.

This example shows how to create a bigimageDatastore that extracts tumor and normal patches for training the network. The example also shows how to preprocess and augment datastores to avoid biasing the network.

Select Location of Normal Tissue Patches to Read

Randomly split the normal images and corresponding masks into two sets. The validation set contains two randomly selected images and corresponding masks. The training set contains the remaining images and masks.

normalValidationIdx = randi(numNormalFiles,[1 2]);
normalTrainIdx = setdiff(1:numNormalFiles,normalValidationIdx);

The patch size is small compared to the size of features in the image. By default, a bigimageDatastore extracts patches with no overlap and no gap, which generates a huge quantity of training patches. You can reduce the amount of training data by specifying a subset of patches. Specify the coordinates of patches using the selectBlockLocations (Image Processing Toolbox) function. Add a gap between the sampled training patches using the BlockOffsets name-value argument. Specify an offset that is larger than the patch size. Increase the inclusion threshold from the default value of 0.5 so that the network trains on relatively homogenous patches.

patchSize = [299,299,3];
normalStrideFactor = 10;
blsNormalData = selectBlockLocations(bigNormalImages(normalTrainIdx), ...
    "BlockSize",patchSize(1:2),"BlockOffsets",patchSize(1:2)*normalStrideFactor, ...
    "Masks",bigNormalMasks(normalTrainIdx),"InclusionThreshold",0.75,"ExcludeIncompleteBlocks",true);

Select the location of validation patches to read. Because there are fewer validation images, you do not need to add a gap between patches.

blsNormalDataValidation = selectBlockLocations(bigNormalImages(normalValidationIdx), ...
    "BlockSize",patchSize(1:2), ...
    "Masks",bigNormalMasks(normalValidationIdx),"InclusionThreshold",0.75,"ExcludeIncompleteBlocks",true);

Create Datastores for Normal Images

Create datastores dsNormalData and dsNormalDataValidation that read image patches from normal images at the finest resolution level for training and validation, respectively. Specify the coordinates of patches using the BlockLocationSet name-value pair argument.

dsNormalData = bigimageDatastore(bigNormalImages(normalTrainIdx), ...
    "BlockLocationSet",blsNormalData);
dsNormalDataValidation = bigimageDatastore(bigNormalImages(normalValidationIdx), ...
    "BlockLocationSet",blsNormalDataValidation);

Preview patches from the datastore containing normal training images.

imagesToPreview = zeros([patchSize 10],'uint8');
for n = 1:10
    im = read(dsNormalData);
    imagesToPreview(:,:,:,n) = im{1};
end
figure
montage(imagesToPreview,'Size',[2 5],'BorderSize',10,'BackgroundColor','k');
title("Training Patches of Normal Tissue")

Select Location of Tumor Tissue Patches to Read

Randomly split the tumor images and corresponding masks into two sets. The validation set contains two randomly selected images and corresponding masks. The training set contains the remaining images and masks.

tumorValidationIdx = randi(numTumorFiles,[1 2]);
tumorTrainIdx = setdiff(1:numTumorFiles,tumorValidationIdx);

Specify the coordinates of patches to read using the selectBlockLocations function. Tumor tissue is more sparse than normal tissue, so increase the sampling density by specifying a smaller block offset than for normal tissue. Note that if you want to train using fewer training images, then you might need to increase the size of the training set by decreasing the block offset even further.

tumorStrideFactor = 4;
blsTumorData = selectBlockLocations(bigTumorImages(tumorTrainIdx), ...
    "BlockSize",patchSize(1:2),"BlockOffsets",patchSize(1:2)*tumorStrideFactor, ...
    "Masks",bigTumorMasks(tumorTrainIdx),"InclusionThreshold",0.75,"ExcludeIncompleteBlocks",true);

Select the location of validation patches to read. Because there are fewer validation images, you do not need to add a gap between patches.

blsTumorDataValidation = selectBlockLocations(bigTumorImages(tumorValidationIdx), ...
    "BlockSize",patchSize(1:2), ...
    "Masks",bigTumorMasks(tumorValidationIdx),"InclusionThreshold",0.75,"ExcludeIncompleteBlocks",true);

Create Datastores for Tumor Images

Create a bigimageDatastore from the training tumor images and masks. The datastores dsTumorData and dsTumorDataValidation read image patches from tumor images at the finest resolution level for training and validation, respectively.

dsTumorData = bigimageDatastore(bigTumorImages(tumorTrainIdx), ...
    "BlockLocationSet",blsTumorData);
dsTumorDataValidation = bigimageDatastore(bigTumorImages(tumorValidationIdx), ...
    "BlockLocationSet",blsTumorDataValidation);

Preview patches from the datastore containing tumor training images.

imagesToPreview = zeros([patchSize 10],'uint8');
for n = 1:10
    im = read(dsTumorData);
    imagesToPreview(:,:,:,n) = im{1};
end
montage(imagesToPreview,'Size',[2 5],'BorderSize',10,'BackgroundColor','k');
title("Training Patches of Tumor Tissue")

Normalize Colors and Augment Training Data

The training images have different color distributions because the data set came from different sources and color staining the tissue does not result in identically stained images. Additional preprocessing is necessary to avoid biasing the network.

To prevent color variability, this example preprocesses the data with standard stain normalization techniques. Apply stain normalization and augmentation by using the transform function with custom preprocessing operations specified by the helper function augmentAndLabel. This function is attached to the example as a supporting file.

The augmentAndLabel function performs these operations:

Each image patch generates five augmented and labeled patches: the stain-normalized patch, the stain-normalized patch with color jitter, the stain-normalized patch with color jitter and random affine transformation, the hematoxylin image with random affine transformation, and the eosin image with random affine transformation.

Create datastores that transform the normal training and validation images and label the generated patches as 'normal'.

dsLabelledNormalData = transform(dsNormalData, ...
    @(x,info)augmentAndLabel(x,info,'normal'),'IncludeInfo',true);
dsLabelledNormalDataValidation = transform(dsNormalDataValidation, ...
    @(x,info)augmentAndLabel(x,info,'normal'),'IncludeInfo',true);

Create datastores that transform the tumor training and validation images and label the generated patches as 'tumor'.

dsLabelledTumorData = transform(dsTumorData, ...
    @(x,info)augmentAndLabel(x,info,'tumor'),'IncludeInfo',true);
dsLabelledTumorDataValidation = transform(dsTumorDataValidation, ...
    @(x,info)augmentAndLabel(x,info,'tumor'),'IncludeInfo',true);

Balance Tumor and Normal Classes

The amount of cancer tissue in the tumor images is very small compared to the amount of normal tissue. Additional preprocessing is necessary to avoid training the network on class-imbalanced data containing a large amount of normal tissue and a very small amount of tumor tissue.

To prevent class imbalance, this example defines a custom datastore called a randomSamplingDatastore that randomly selects normal and tumor training patches in a balanced way. The script to define this custom datastore is attached to the example as a supporting file. For more information, see Develop Custom Datastore.

Create the custom randomSamplingDatastore from the normal and tumor training datastores. The random sampling datastore dsTrain provides mini-batches of training data to the network at each iteration of the epoch.

dsTrain = randomSamplingDatastore(dsLabelledTumorData,dsLabelledNormalData);

To limit the number of patches used during validation, this example defines a custom datastore called a validationDatastore that returns five validation patches from each class. The script to define this custom datastore is attached to the example as a supporting file.

Create the custom validationDatastore from the normal and tumor validation datastores.

numValidationPatchesPerClass = 5;
dsValidation = validationDatastore(dsLabelledTumorDataValidation, ...
    dsLabelledNormalDataValidation,numValidationPatchesPerClass);

Set Up Inception-v3 Network Layers

This example uses the Inception-v3 network, a convolutional neural network that is trained on more than a million images from the ImageNet database [3]. The network is 48 layers deep and can classify images into 1,000 object categories, such as keyboard, mouse, pencil, and many animals. The network expects an image input size of 299-by-299 with 3 channels.

The inceptionv3 function returns a pretrained Inception-v3 network. Inception-v3 requires the Deep Learning Toolbox™ Model for Inception-v3 Network support package. If this support package is not installed, then the function provides a download link.

net = inceptionv3;

Replace Final Layers

The convolutional layers of the network extract image features that the last learnable layer and the final classification layer use to classify the input image. These two layers contain information on how to combine the features that the network extracts into class probabilities, a loss value, and predicted labels. To retrain a pretrained network to classify new images, replace these two layers with new layers adapted to the new data set. For more information, see Train Deep Learning Network to Classify New Images.

Extract the layer graph from the trained network.

lgraph = layerGraph(net);

Find the names of the two layers to replace by using the supporting function findLayersToReplace. This function is attached to the example as a supporting file. In Inception-v3, these two layers are named 'predictions' and 'ClassificationLayer_predictions'.

[learnableLayer,classLayer] = findLayersToReplace(lgraph)

learnableLayer = 
  FullyConnectedLayer with properties:

          Name: 'predictions'

   Hyperparameters
     InputSize: 2048
    OutputSize: 1000

   Learnable Parameters
       Weights: [1000×2048 single]
          Bias: [1000×1 single]

  Show all properties

classLayer = 
  ClassificationOutputLayer with properties:

            Name: 'ClassificationLayer_predictions'
         Classes: [1000×1 categorical]
      OutputSize: 1000

   Hyperparameters
    LossFunction: 'crossentropyex'

The goal of this example is to perform binary segmentation between two classes, tumor and nontumor regions. Create a new fully connected layer for two classes. Replace the original final fully connected layer with the new layer.

numClasses = 2;
newLearnableLayer = fullyConnectedLayer(numClasses,'Name','predictions');
lgraph = replaceLayer(lgraph,learnableLayer.Name,newLearnableLayer);

Create a new classification layer for two classes. Replace the original final classification layer with the new layer.

newClassLayer = classificationLayer('Name','ClassificationLayer_predictions');
lgraph = replaceLayer(lgraph,classLayer.Name,newClassLayer);

Specify Training Options

Train the network using the rmsprop optimization solver. This solver automatically adjusts the learning rate and momentum for faster convergence. Specify other hyperparameter settings by using the trainingOptions function. Reduce MaxEpochs to a small number because the large amount of training data enables the network to reach convergence sooner. If you have Parallel Computing Toolbox™ and the hardware resources available for multi-GPU or parallel training, then accelerate training by specifying the ExecutionEnvironment name-value pair argument as 'multi-gpu' or 'parallel', respectively.

checkpointsDir = fullfile(trainingImageDir,'checkpoints');
if ~exist(checkpointsDir,'dir')
    mkdir(checkpointsDir);
end

options = trainingOptions('rmsprop', ...
    'InitialLearnRate',1e-5, ...
    'SquaredGradientDecayFactor',0.99, ...
    'MaxEpochs',3, ...
    'MiniBatchSize',32, ...
    'Plots','training-progress', ...
    'CheckpointPath',checkpointsDir, ... 
    'ValidationData',dsValidation, ...
    'ExecutionEnvironment','auto', ...
    'Shuffle','every-epoch');

Download and Preprocess Test Data

The Camelyon16 test data set consists of 130 WSIs. These images have both normal and tumor tissue. This example uses two test images from the Camelyon16 test data. The size of each file is approximately 2 GB.

Create a directory to store the test data.

testingImageDir = fullfile(tempdir,'Camelyon16','testing');
if ~exist(testingImageDir,'dir')
    mkdir(testingImageDir);
    mkdir(fullfile(testingImageDir,'images'));
    mkdir(fullfile(testingImageDir,'lesion_annotations'));
end

testDataDir = fullfile(testingImageDir,'images');
testTumorAnnotationDir = fullfile(testingImageDir,'lesion_annotations');

To download the test data, go to the Camelyon17 website and click the first "CAMELYON16 data set" link. Open the "testing" directory, then follow these steps.

Specify the number of test images.

numTestFiles = 2;

Create Masks for Test Images

The test images contain a mix of normal and tumor images. To reduce the amount of computation during classification, define the ROIs by creating masks.

Specify a resolution level to use for creating the mask. This example uses resolution level 7, which is coarse and fits in memory.

resolutionLevel = 7;

Create masks for regions containing tissue. You can use the helper function createMaskForNormalTissue to create masks using color thresholding. This helper function is attached to the example as a supporting file. For more information about this helper function, see Create Masks for Normal Images.

testTissueMaskDir = fullfile(testDataDir,['test_tissuemask_level' num2str(resolutionLevel)]);
createMaskForNormalTissue(testDataDir,testTissueMaskDir,resolutionLevel);

Create masks for images that contain tumor tissue. Skip images that do not contain tumor tissue. You can use the helper function createMaskForTumorTissue to create masks using ROI objects. This helper function is attached to the example as a supporting file. For more information about this helper function, see Create Masks for Tumor Images.

testTumorMaskDir = fullfile(testDataDir,['test_tumormask_level' num2str(resolutionLevel)]);
createMaskForTumorTissue(testDataDir,testTumorAnnotationDir,testTumorMaskDir,resolutionLevel);

Train Network

After configuring the training options and the data source, train the Inception-v3 network by using the trainNetwork function. By default, this example downloads a pretrained version of the Inception-v3 network for this data set by using the helper function downloadTrainedCamelyonNet. The helper function is attached to the example as a supporting file. The pretrained network enables you to run the entire example without waiting for training to complete.

Note: Training takes 20 hours on an NVIDIA™ Titan X and can take even longer depending on your GPU hardware.

doTraining = false;
if doTraining
    trainedNet = trainNetwork(dsTrain,lgraph,options);
    save(['trainedCamelyonNet-' datestr(now, 'dd-mmm-yyyy-HH-MM-SS') '.mat'],'trainedNet');
else
    trainedCamelyonNet_url = 'https://www.mathworks.com/supportfiles/vision/data/trainedCamelyonNet.mat';
    netDir = fullfile(tempdir,'Camelyon16');
    downloadTrainedCamelyonNet(trainedCamelyonNet_url,netDir);
    load(fullfile(netDir,'trainedCamelyonNet.mat'));
end

Downloading pretrained Inception-v3 network for Cameylon16 data.
This can take several minutes to download...
Done.

Classify Test Data and Create Heatmaps

A GPU is highly recommended for classification of the test images (requires Parallel Computing Toolbox™).

Each test image has two masks, one indicating normal tissue and one indicating tumor tissue. Create bigimage objects to manage the test data and masks by using the helper function createBigImageAndMaskArrays. The helper function is attached to the example as a supporting file.

[bigTestImages,bigTestTissueMasks] = createBigImageAndMaskArrays(testDataDir,testTissueMaskDir);
[~,bigTestTumorMasks] = createBigImageAndMaskArrays(testDataDir,testTumorMaskDir);

Use the trained Inception-v3 network to identify tumor patches in the test images, bigTestImages. Classify the test images on a block-by-block basis by using the apply (Image Processing Toolbox) function with a custom processing pipeline specified by the helper function createHeatMap. This helper function is attached to the example as a supporting file. To reduce the amount of computation required, specify the tissue mask bigTestTissueMask so that the apply function processes only patches that contain tissue. If you have Parallel Computing Toolbox™ and a GPU, then you can evaluate blocks in parallel by specifying the 'UseParallel' name-value pair argument as true.

The createHeatMap function performs these operations on each tissue block:

The apply function stitches the heatmap of each block into a single heatmap for the test image. The heatmap shows where the network detects regions containing tumors.

To visualize the heatmap, overlay the heatmap on the original image and set the transparency 'AlphaData' property as the tissue mask. The overlay shows how well the tumor is localized in the image. Regions with a high probability of being tumors are displayed with red pixels. Regions with a low probability of being tumors are displayed as blue pixels.

patchSize = [299,299,3];

for idx = 1:numTestFiles
    bigTestHeatMaps(idx) = apply(bigTestImages(idx),1,@(x)createHeatMap(x,trainedNet), ...
        'Mask',bigTestTissueMasks(idx),'InclusionThreshold',0, ...
        'BlockSize',patchSize(1:2),'UseParallel',false);
    
    figure
    hBigTest = bigimageshow(bigTestImages(idx));
    hTestAxes = hBigTest.Parent;
    hTestAxes.Visible = 'off';
        
    hMaskAxes = axes;
    hBigMask = bigimageshow(bigTestHeatMaps(idx),'Parent',hMaskAxes, ...
        "Interpolation","nearest","AlphaData",bigTestTissueMasks(idx));
    colormap(jet(255));
    hMaskAxes.Visible = 'off';
    
    linkaxes([hTestAxes,hMaskAxes]);
    title(['Tumor Heatmap of Test Image ',num2str(idx)])
end

References

[1] Ehteshami B. B., et al. "Diagnostic Assessment of Deep Learning Algorithms for Detection of Lymph Node Metastases in Women With Breast Cancer." Journal of the American Medical Association. Vol. 318, No. 22, 2017, pp. 2199–2210. doi:10.1001/jama.2017.14585

[2] Szegedy, C., V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna. "Rethinking the Inception Architecture for Computer Vision." In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2818–2826. Las Vegas, NV: IEEE, 2016.

[3] ImageNet. http://www.image-net.org.

[4] Macenko, M., et al. "A Method for Normalizing Histology Slides for Quantitative Analysis." In 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, 1107–1110. Boston, MA: IEEE, 2009.

[5] xml2struct https://www.mathworks.com/matlabcentral/fileexchange/28518-xml2struct.