%% Deep Learning Classification of Large Multiresolution Images
% This example shows how to train an Inception-v3 deep neural network to classify multiresolution whole slide images (WSIs) that do not fit in memory.
% The example shows how to train an Inception-v3 tumor classification network and also provides a pretrained Inception-v3 network. If you choose to train the Inception-v3 network, use of a CUDA capable NVIDIA™ GPU with compute capability 3.0 or higher is highly recommended. Use of a GPU requires Parallel Computing Toolbox™.
%
%% 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 objects can store and process this type of large multiresolution image. bigimageDatastore 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.
% Download the "lesion_annotations.zip" file. Extract the files to the directory specified by the trainTumorAnnotationDir variable.
% Open the "normal" directory. Download the images to the directory specified by the trainNormalDataDir variable.
% Open the "tumor" directory. Download the images to the directory specified by the trainTumorDataDir variable.
% 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 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"])

% Display the bigimage at a coarse resolution level by using the bigimageshow 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 and Process Big Images Efficiently Using 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 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:
% Create a bigimage object from the .TIF image file.
% Set the spatial referencing of all resolution levels from the image metadata.
% Get the image at the coarse resolution level.
% Convert the coarse image to the L*a*b* color space, then extract the a* channel.
% Create a binary image by thresholding the image using Otsu's method, which minimizes the intraclass variance between the black and white pixels.
% Create a single-resolution bigimage object from the mask and set the spatial referencing of the mask to match the spatial referencing of the input image.
% Write the mask bigimage to memory. Only the bigimage object is in memory. The individual image blocks corresponding to the logical mask image are in a temporary directory. Writing to a directory preserves the custom spatial referencing, which ensures that the normal images and their corresponding mask images have the same spatial referencing.

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)
bigNormalMasks(idx).SpatialReferencing(1)

%
% 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:
% Create a bigimage object from the .TIF image file.
% Set the spatial referencing from the image metadata.
% Read the corresponding lesion annotations in XML files and convert the annotations to polygons (Polygon objects).
% For each image block, use the polygon data to create a mask for the corresponding block. Images with tumor regions can contain some normal regions within them. Use the annotations of normal tissue to exclude those regions.
% Create an output logical mask bigimage object at a coarser resolution level. Write the mask image on a block-by-block basis by using the setBlock function.
% Write the mask bigimage object to a directory in memory. Only the bigimage object is in memory. The individual image blocks corresponding to the logical mask image are in a temporary directory. Writing to a directory preserves the custom spatial referencing, which ensures that the tumor images and their corresponding mask images have the same spatial referencing.

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)
bigTumorMasks(idx).SpatialReferencing(1)

%
% 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 Datastore for Training and Validation
% To extract patches of training data from the bigimage objects, use a bigimageDatastore. This datastore reads patches of bigimage data at a single resolution level. If you specify a mask, then the datastore reads patches only from within the ROI blocks. Note that each mask block can contain a mixture of ROI and background pixels, which can lead to ambiguity about whether the entire block is ROI or background. To remove the ambiguity, the datastore has an InclusionThreshold property. This property enables you to specify the minimum fraction of true pixels within a mask block necessary for the block to be considered ROI.
% Two main issues with the raw training patches can potentially bias the network. The first issue is color imbalance resulting from nonuniform color staining of the tissue. The second issue is class imbalance resulting from an unequal amount of tumor and normal tissue in the data. To correct these issues, you can preprocess and augment the datastore.
% This example shows how to create a bigimageDatastore to extract tumor and normal patches for training the network and how to preprocess and augment the datastores to avoid biasing the network.
%
%% Create Datastore for Normal Images
% 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 using the BlockOffsets property to add a gap between sampled patches. Specify a block 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.
% The datastores dsNormalData and dsNormalDataValidation read image patches from normal images at the finest resolution level for training and validation, respectively.

patchSize = [299,299,3];

normalStrideFactor = 10;
dsNormalData = bigimageDatastore(bigNormalImages(normalTrainIdx),1, ...
    'InclusionThreshold',0.75,'Mask',bigNormalMasks(normalTrainIdx), ...
    'BlockSize',patchSize(1:2),'BlockOffsets',patchSize(1:2)*normalStrideFactor, ...
    'IncompleteBlocks','exclude');

dsNormalDataValidation = bigimageDatastore(bigNormalImages(normalValidationIdx),1, ...
    'InclusionThreshold',0.75,'Mask',bigNormalMasks(normalValidationIdx), ...
    'BlockSize',patchSize(1:2),'IncompleteBlocks','exclude');

% 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")

%
%% Create Datastore for Tumor Images
% Randomly split the tumor images and corresponding masks into two sets. The validation set contains two images and masks. The training set contains the remaining images and masks.

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

%
% Create a bigimageDatastore from the training tumor images and masks. Tumor tissue is more sparse than normal tissue. To increase the sampling density and generate more training patches, specify a smaller block offset. 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.
% The datastores dsTumorData and dsTumorDataValidation read image patches from tumor images at the finest resolution level for training and validation, respectively.

tumorStrideFactor = 4;
dsTumorData = bigimageDatastore(bigTumorImages(tumorTrainIdx),1, ...
    'InclusionThreshold',0.75,'Mask',bigTumorMasks(tumorTrainIdx), ...
    'BlockSize',patchSize(1:2),'BlockOffsets',patchSize(1:2)*tumorStrideFactor, ...
    'IncompleteBlocks','exclude');

dsTumorDataValidation = bigimageDatastore(bigTumorImages(tumorValidationIdx),1, ...
    'InclusionThreshold',0.75,'Mask',bigTumorMasks(tumorValidationIdx), ...
    'BlockSize',patchSize(1:2),'IncompleteBlocks','exclude');

% 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:
% Normalize staining by using the normalizeStaining.m function [4]. Stain normalization is performed using Macenko's method, which separates H&E color channels by color deconvolution using a fixed matrix and then recreates the normalized images with individual corrected mixing. The function returns the normalized image as well as the H&E images.
% Add color jitter by using the jitterColorHSV function. Color jitter varies the color of each patch by perturbing the image contrast, hue, saturation, and brightness. Color jitter is performed in the HSV color space to avoid unwanted color artifacts in the RGB image.
% Apply random combinations of 90 degree rotations and vertical and horizontal reflection. Randomized affine transformations make the network agnostic to the orientation of the input image data.
% Label the patch as 'normal' or 'tumor'.
% 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)

% 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.
% Download the "lesion_annotations.zip" file. Extract all files to the directory specified by the testTumorAnnotationDir variable.
% Open the "images" directory. Download the first two files, "test_001.tif" and "test_002.tif". Move the files to the directory specified by the testDataDir variable.
% 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

%
%% 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 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:
% Calculate the tumor probability score by using the predict function.
% Create a heatmap image patch with pixel values equal to the tumor probability score.
% 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.
% Copyright 2019 The MathWorks, Inc.