%% Segment Satellite Images Using Deep Learning
% This example shows how to perform semantic segmentation of buildings in high-resolution satellite images using a Deeplab v3+ neural network.
% The example also provides a pretrained DeepLab v3+ network with weights initialized from a pretrained Resnet-50 network. For more details, see Pretrained Deep Neural Networks (Deep Learning Toolbox). If you choose to train the network, use of a quad CUDA capable NVIDIA™ GPU with compute capability 3.0 or higher is highly recommended. Use of a GPU requires Parallel Computing Toolbox™.
%
%% Introduction
% Building footprint detection is important for several applications such as disaster relief and urban planning. To analyze satellite images, images must be labeled so that buildings are distinct from other landmarks in the images such as roads, rivers, and agricultural areas. Deep learning methods enable automated labeling of high-resolution satellite images. This example demonstrates two different strategies for segmenting buildings.
% The first approach uses the building boundary labels to define a trinary segmentation with the classes "building boundary", "building interior", and background. However, it is difficult to classify pixels on building boundaries. Misclassified boundary pixels can result in segmented images whose building boundaries are fuzzy or incorrectly located.
% The second approach increases the accuracy of building borders by training using a distance objective. The distance objective works well especially when several buildings are very close to each other or even contiguous. The distance objective relies on a signed distance pixel values indicate the distance to the closest building boundary [6]. Pixels on the building boundary are labeled 0. Pixels inside buildings have a positive value and pixels outside the building have a negative value. The further a pixel is from the boundary of a building, the greater the magnitude of its value. The distance values are quantized to the range [-64 64] and then shifted to all positive values to the range [1 129]. This new 129-class segmentation problem makes the segmentation more fine-grained and building boundaries more accurate.
% Training both types of segmentation networks benefits from using many patches from the same image. Because the images are in TIFF files, you can use a bigimageDatastore to extract several patches from each training image. To use a bigimageDatastore, you must first store information about the images in bigimage objects. A bigimage manages the loading and processing of small image blocks that are out of memory.
% This example shows how to train a deep learning network to perform semantic segmentation of buildings in satellite images using bigimage and bigimageDatastore. The example presents semantic segmentation results as overlay of segmentation pixel labels and corresponding images to show the location of buildings. The display of building pixels allows the GIS practitioners to quickly identify buildings.
%
%% Download Training Data
% This example uses high-resolution satellite images of Las Vegas from the first and second SpaceNet challenges [1].
% The images were captured from the WorldView-3 satellite. This satellite collects multiple types of images, including panchromatic images with a 0.31 m resolution and 8-channel multispectral images at a lower 1.24 m resolution. The SpaceNet data set also includes pansharpened RGB images, in which the resolution of the color channels is enhanced using information from the high-resolution panchromatic images. This example trains the segmentation network using the pansharpened RGB images.
% The SpaceNet data is hosted as an Amazon Web Services (AWS) Public Dataset. To download the data set, create an Amazon AWS™ account and download the Vegas buildings dataset at the SpaceNet AWS portal. For detailed descriptions of data formats and download instructions click on the link here.
% Create a directory to store the training data set.

spaceNetDir = fullfile(tempdir,'SpaceNet');

% Download the training data set, AOI_2_Vegas_Train.tar.gz. Extract the files to the directory specified by the spaceNetDir variable. The extracted directory, 'AOI_2_Vegas_Train' includes six total subdirectories. This example uses two of the subdirectories:
% 'RGB-PanSharpen' contains 3,851 pansharpened RGB images. Each image is a 650-by-650 pixel block representing a 200-by-200 square-meter region of Las Vegas. Images represent both urban and suburban areas. Image files are in the TIF file format with uint16 data type.
% 'summaryData' contains pixel labels for buildings in each RGB image. Labels are stored in a CSV file.
% The size of each image file is approximately 2.5 MB and the entire training data set is about 100 GB.
% Specify the location of the data set.

spaceNetTrainDir = fullfile(spaceNetDir,'AOI_2_Vegas_Train');
trainImageDir = fullfile(spaceNetTrainDir,'RGB-PanSharpen');

% Create a new directory, 'training', to store intermediate images used to train the network.

preprocessDataDir = fullfile(spaceNetDir,'preprocessedDataset');
trainingDir = preprocessDataDir;

if ~exist(trainingDir,'dir')
    mkdir(trainingDir);
    mkdir(fullfile(trainingDir,'maskImage'));
    mkdir(fullfile(trainingDir,'distanceLabelImage'));
    mkdir(fullfile(trainingDir,'buildingLabelImage'));
end

trainMaskDir = fullfile(trainingDir,'maskImage');
trainDistanceLabelDir = fullfile(trainingDir,'distanceLabelImage');
trainBuildingLabelDir = fullfile(trainingDir,'buildingLabelImage');

% The SpaceNet test data for Las Vegas does not include ground truth building labels so this example does not use that test data. Instead, this example splits the training data into three independent sets that are used for training, validation, and testing. Specify the number of each type of image. To speed up the training process, use a portion of the data set.

numTrain = 800;
numVal = 300;
numTest = 10;
numImages = numTrain + numVal + numTest;

%
%% Create Masks of Satellite Images
% Some satellite images contain regions that are not useful for training, such as black borders or large landmarks without buildings. Exclude these regions from training by using masks.
% The RGB color of satellite imagery is distinct from the color of the black 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 image mask.
% To create masks using color thresholding, use the helper function createSpaceNetImageMasks. This helper function is attached to the example as a supporting file. The function performs these operations on each training image that has building polygons. The function does not create masks for images with no corresponding building polygons.
% Read each TIF image file.
% Convert the image to the L*a*b* color space, then extract the a* channel. Black borders and dark regions such as shadows have a large a* value and appear bright. Buildings appear bright against vegetation and other landmarks because small values of a* correspond to green hues whereas large values of a* correspond to red hues.
% Create a binary image by thresholding the image using Otsu's method, which minimizes the intraclass variance between the black and white pixels. In the binary image, borders and buildings are white in the binary image and vegetation is black.
% Perform morphological opening on the binary image. Use a large disk-shaped structuring element of 100 pixels. Buildings and shadows that have a small size vanish during the erosion phase and do not return during the dilation phase. Large borders continue to exist after the erosion, and return to their original size after the dilation.
% Specify the mask as the inverse of the opened binary image. In the mask image, the white border is black (false) and the desired ROI consisting of buildings and landmarks are white (true).
% Write the mask image to memory as a TIF file.

createSpaceNetImageMasks(spaceNetTrainDir,trainMaskDir);

%
% Proper images for training, validation, and testing have associated masks. These steps create an image datastore that contains only proper, or valid, files.
% First, specify the set of all image file names by using an ImageDatastore.

imds = imageDatastore(trainImageDir,'FileExtensions',{'.tif'});

% Shuffle the datastore so that the datastore reads image files in a random order.

rng default
imds = shuffle(imds);

% Get the indices of valid images in the datastore. The datastore has more valid images than are needed for training, validation, and testing, so only return the needed number of images.

validIndices = getValidSpaceNetIndices(imds,trainMaskDir,numImages);

% Create a new ImageDatastore that contains only files specified by valid indices.

imdsValidImages = subset(imds,validIndices);

%
% Create arrays of bigimage objects that manage the valid image and mask data.

bigArrayIdx = 1;
while hasdata(imdsValidImages)

    [~,info] = read(imdsValidImages);
    imageFullFileName = info.Filename;    
    [~,name] = fileparts(imageFullFileName);
    maskImageFileName = strcat(name,'_mask');
    
    bigImages(bigArrayIdx) = bigimage(imageFullFileName); 
    bigMasks(bigArrayIdx) = bigimage(fullfile(trainMaskDir,[maskImageFileName '.tif']));
    
    bigArrayIdx = bigArrayIdx+1;
end

%
% Display a sample image with the mask overlay. This example displays the twelfth image because this image includes  Because the image and mask are small enough to fit in memory, you can bring the image and mask into memory by using the getFullLevel function.

displayImageIdx = 12;
RGB = getFullLevel(bigImages(displayImageIdx));
mask = getFullLevel(bigMasks(displayImageIdx));

% The SpaceNet images are stored as uint16 TIFF images. To display, rescale the intensity values input image as these values are stored in a very small range.

RGB = rescale(RGB);

% display the corresponding mask image with partial transparency. The mask highlights the regions containing building images.

RGBmasked = labeloverlay(RGB,mask,'IncludedLabels',[0 1],'Transparency',0.8);
imshow(RGBmasked)

%
%% Create Ground Truth Label Images
% Create ground truth pixel label images using the building footprint corresponding to each image. These are stored as a CSV file with these fields: ImageId, BuildingId, PolygonWKT_Pix and PolygonWKT_Geo. The PolygonWKT_Pix field defines the building polygon intrinsic coordinates in Well-Known Text (WKT) format [5].
% To create label images from the polygon coordinates, use the helper function createSpaceNetLabelImages. This helper function is attached to the example as a supporting file.
% The function performs these preprocessing operations to extract the polygon coordinates. The function does not create label images for images with no corresponding building polygons.
% To load the polygon coordinates from the CSV file, create a TabularTextDatastore then read all of the data into a table.
% Read each TIF image file and the corresponding polygon labels. The labels are a cell array of strings. Each string denotes polygon coordinates for a distinct building in the image.
% Parse the strings and extract the coordinates into n-by-2 matrices of data type double. Round the coordinates to the nearest integer value within the bounds of the image.
% The function performs these additional steps to create the building label image.
% Create grids covering all (x,y) image coordinates.
% Query if each point is within or on the edge of polygons by using the inpolygon function. The result is a binary mask in which pixels within or on the boundary of a polygon are white (true) and pixels outside a polygon are black (false).
% To identify pixels within a polygon, erode the binary segmentation with a small disk shaped structuring element with a radius of three pixels. The result is a binary mask in which pixels within the boundary of a polygon are white (true) and other pixels are black (false).
% To identify pixels on the boundary of a polygon, subtract the eroded image from the binary segmentation. The result is a binary mask in which pixels on the boundary of a polygon are white (true) and other pixels are black (false).
% Combine the three classes into a single numeric image by taking a linear combination of the segmentation masks of the three pixel types.
% Quantize the image to three pixel values: 1, 2, and 3.
% The function performs these additional steps to create a distance label image, starting with the binary mask of the boundary pixels.
% For each pixel in the image, calculate the Euclidean distance between that pixel and the nearest nonzero pixel by using the bwdist function. The result is a matrix of nonnegative numbers of data type single.
% Create a signed image by flipping the sign of pixels outside a polygon.
% Quantize the signed image to 129 discrete integer values in the range [-64, 63] by using the imquantize function.
% Finally, the function writes the building label images and distance label images to memory as a TIF files.

createSpaceNetLabelImages(spaceNetTrainDir,trainBuildingLabelDir,trainDistanceLabelDir)

% Specify class names for the building label images. The three classes signify building exteriors, building boundaries, and building interiors, respectively.

buildingClassNames = ["notBuilding","boundary","building"];

% Specify class names for the distance label images. The 129 classes represent quantized distances from building borders with an additive offset of 65. The offset shifts the original range of distance values from  [-64, 64] to all positive values in the range [1, 129]. Border pixels have the class 65. Pixels in building interiors have class names in the range [66, 129]. Pixels outside buildings have class names in the range [1, 64].

distanceClassNames = split(string(num2str(1:129)));

%
% Now that the building label images and distance label images are on disk, manage the data using bigimage objects.
% Specify the set of image file names by using the same ImageDatastore as for creating the mask images. To read the same set of indices in the same order, reset the datastore.

reset(imdsValidImages);

% Create an array of categorical bigimage objects from the building label images and distance label images.

bigArrayIdx = 1;
while hasdata(imdsValidImages)

    [~,info] = read(imdsValidImages);
    imageFullFileName = info.Filename;    
    [~,name] = fileparts(imageFullFileName);
    buildingLabelFileName = strcat(name,'_label');
    distanceLabelFileName = strcat(name,'_distance'); 
    
    bigBuildingLabels(bigArrayIdx) = bigimage(fullfile(trainBuildingLabelDir,[buildingLabelFileName '.tif']), ...
        'Classes',buildingClassNames,'PixelLabelIDs',[1 2 3],'UnloadedValue',"notBuilding",'UndefinedID',0);
    bigDistanceLabels(bigArrayIdx) = bigimage(fullfile(trainDistanceLabelDir,[distanceLabelFileName '.tif']), ...
        'Classes',distanceClassNames,'PixelLabelIDs',1:129,'UnloadedValue','1','UndefinedID',0);
    
    bigArrayIdx = bigArrayIdx+1;
end

%
% Display a sample label image and distance label image.
% The building label images are represented as categorical big images with 3 classes. To display the building label image, convert the data type to double and rescale the data to the range [0, 1] by using the apply function. The "notBuilding", "boundary", and "building" classes display as black, white, and gray, respectively.

displayImageIdx = 1;

bldg = getFullLevel(bigBuildingLabels(displayImageIdx));
bldg = rescale(double(bldg));
imshow(bldg)

% The distance label images are categorical big images with 129 classes. To display the distance label image as a grayscale image, convert the data data type to double and rescale the data to the range [0, 1] by using the apply function

dis = getFullLevel(bigDistanceLabels(displayImageIdx));
dis = rescale(double(dis));
imshow(dis)

%
%% Create Big Image Datastores 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.
% Specify a patch size of 250-by-250 pixels to match the input size of the network.

patchSize = [256,256,3];
strideFactor = 1;

% Create a bigimageDatastore object to read the satellite images. To exclude patches from the black border regions, specify the 'Masks' property as the corresponding mask from bigTrainMasks.

trainImageIndices = 1:numTrain;
bigimdsTrain = bigimageDatastore(bigImages(trainImageIndices),1, ...
    'InclusionThreshold',1,'Mask',bigMasks(trainImageIndices), ...
    'BlockSize',patchSize(1:2),'BlockOffsets',patchSize(1:2)*strideFactor, ...
    'IncompleteBlocks','exclude');

% The SpaceNet images are uint16 images with intensity values occupying only a small portion of the range. Preprocess the images by normalizing the images based on the range of the intensity values.

bimdsTrainNormalized = transform(bigimdsTrain,@(x)mat2gray(x{:}));

% Preview patches from the datastore containing training images.

imagesToPreview = zeros([patchSize 10]);
for n = 1:10
    im = read(bimdsTrainNormalized);
    imagesToPreview(:,:,:,n) = im;
end
montage(imagesToPreview,'Size',[2 5],'BorderSize',10,'BackgroundColor','w');

%
% Create two bigimageDatastore objects to read the building label images and distance label images. Use the same masks as for the RGB images.

bigimdsBuildingTrain = bigimageDatastore(bigBuildingLabels(trainImageIndices),1, ...
    'InclusionThreshold',1,'Mask',bigMasks(trainImageIndices), ...
    'BlockSize',patchSize(1:2),'BlockOffsets',patchSize(1:2)*strideFactor, ...
    'IncompleteBlocks','exclude');
bigimdsDistanceTrain = bigimageDatastore(bigDistanceLabels(trainImageIndices),1, ...
    'InclusionThreshold',1,'Mask',bigMasks(trainImageIndices), ...
    'BlockSize',patchSize(1:2),'BlockOffsets',patchSize(1:2)*strideFactor, ...
    'IncompleteBlocks','exclude');

% Combine the corresponding image and label bigimageDatastore objects. The combined datastores are used to train two networks.

dsTrainBuilding = combine(bimdsTrainNormalized,bigimdsBuildingTrain);
dsTrainDistance = combine(bimdsTrainNormalized,bigimdsDistanceTrain);

%
% Create combined datastores for validation using a similar processing pipeline as for the training data.

validationImageIndices = numTrain+1:numTrain+numVal;

bigimdsVal = bigimageDatastore(bigImages(validationImageIndices),1, ...
    'InclusionThreshold',1,'Mask',bigMasks(validationImageIndices), ...
    'BlockSize',patchSize(1:2),'BlockOffsets',patchSize(1:2)*strideFactor, ...
    'IncompleteBlocks','exclude');
bigimdsValNormalized = transform(bigimdsVal,@(x)mat2gray(x{:}));

bigimdsBuildingVal = bigimageDatastore(bigBuildingLabels(validationImageIndices),1, ...
    'InclusionThreshold',1,'Mask',bigMasks(validationImageIndices), ...
    'BlockSize',patchSize(1:2),'BlockOffsets',patchSize(1:2)*strideFactor, ...
    'IncompleteBlocks','exclude');
bigimdsDistanceVal = bigimageDatastore(bigDistanceLabels(validationImageIndices),1, ...
    'InclusionThreshold',1,'Mask',bigMasks(validationImageIndices), ...
    'BlockSize',patchSize(1:2),'BlockOffsets',patchSize(1:2)*strideFactor, ...
    'IncompleteBlocks','exclude');
dsValBuilding = combine(bigimdsValNormalized,bigimdsBuildingVal);
dsValDistance = combine(bigimdsValNormalized,bigimdsDistanceVal);

%
%% Set Up Deeplabv3 Network Layers
% This example uses the Deeplab v3+ network, a convolutional neural network for semantic image segmentation. The network uses encoder-decoder architecture, dilated convolutions, and skip connections to segment images. The network is based on ResNet-50 which is a convolutional neural network that is trained on more than a million images from the ImageNet database [4]. The network is 50 layers deep and can classify images into 1,000 object categories such as keyboard, mouse, pencil, and many animals. As a result, the network has learned rich feature representations for a wide range of images. The image input size is 250-by-250 with 3 channels.
% Use the deeplabv3plusLayers function to create a DeepLab v3+ network based on ResNet-50. This requires the Deep Learning Toolbox™ Model for ResNet-50 Network support package for the pretrained weights. If this support package is not installed, then the function provides a download link.

% TODO remove this line in the HTML edit
addpath(fullfile(matlabroot,'toolbox','nnet','supportpackages','resnet50'))

% Create a layer graph for building labels. Specify the second argument as the number of classes in the trinary building segmentation, 3.

lgraphBuilding = deeplabv3plusLayers(patchSize,3,'resnet50');

% Create a layer graph for distance labels. Specify the second argument as the number of classes in the distance segmentation, 129.

lgraphDistance = deeplabv3plusLayers(patchSize,129,'resnet50');

%
%% Define Training Options
% Train the network using the adam optimization solver. This solver automatically adjusts the learning rate and momentum for faster convergence. Specify other hyperparameter settings by using the trainingOptions function.
% 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.

l2reg = 0.0005;
maxEpochs = 50;

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

options = trainingOptions('adam', ...
    'InitialLearnRate',3e-4, ...
    'SquaredGradientDecayFactor',0.99, ...
    'MaxEpochs',maxEpochs, ...
    'MiniBatchSize',32, ...
    'Shuffle','every-epoch', ...
    'Plots','training-progress', ...
    'ExecutionEnvironment','auto', ...
    'CheckpointPath',checkpointsDir);

% Except for the validation data, use identical training parameters for the building segmentation network and distance segmentation network.

optionsBuilding = options;
optionsBuilding.ValidationData = dsValBuilding;

optionsDistance = options;
optionsDistance.ValidationData = dsValDistance;

%
%% Train Networks
% After configuring the training options and the data source, train two DeepLab v3+ networks by using the trainNetwork function.
% By default, this example downloads two pretrained versions of the DeepLab v3+ network for this data set by using the helper function downloadTrainedSpaceNetBuilding and downloadTrainedSpaceNetDistance, respectively. The helper functions are attached to the example as supporting files. The pretrained networks enable you to run the entire example without waiting for training to complete.
% Train the building segmentation network using the 3-class building label images.
% Note: Training takes 10 hours on a quad NVIDIA™ Titan X and can take even longer depending on your GPU hardware.

doTraining = false;
if doTraining
    trainedNetBuilding = trainNetwork(dsTrainBuilding,lgraph,optionsBuilding);
    trainedNetBuildingStr = regexprep(['trainedNetBuilding_' datestr(now, 'dd-mmm-yyyy-HH-MM-SS')],'-','_');
    evalin('base',strcat(trainedNetBuildingStr, ' = trainedNetBuilding;'));
    save(['deeplabv3plusResNet50SpaceNetBuilding-' datestr(now, 'dd-mmm-yyyy-HH-MM-SS') '.mat'],trainedNetBuildingStr,'options')
else
    trainedSpaceNetBuildingURL = 'https://www.mathworks.com/supportfiles/vision/data/deeplabv3plusResNet50SpaceNetBuilding.mat';
    netDir = fullfile(tempdir,'SpaceNet');
    downloadTrainedSpaceNetBuildingNet(trainedSpaceNetBuildingURL,netDir);
    load(fullfile(netDir,'deeplabv3plusResNet50SpaceNetBuilding.mat'));
end

% Train the distance segmentation network using the 129-class distance label images
% Note: Training takes 10 hours on a quad NVIDIA™ Titan X and can take even longer depending on your GPU hardware.

doTraining = false;
if doTraining
    trainedNetDistance = trainNetwork(dsTrainDistance,lgraph,optionsDistance);
    trainedNetDistanceStr = regexprep(['trainedNetDistance_' datestr(now, 'dd-mmm-yyyy-HH-MM-SS')],'-','_');
    evalin('base',strcat(trainedNetDistanceStr, ' = trainedNetDistance;'));
    save(['deeplabv3plusResNet50SpaceNetDistance-' datestr(now, 'dd-mmm-yyyy-HH-MM-SS') '.mat'],trainedNetDistance,'options')
else
    trainedSpaceNetDistanceURL = 'https://www.mathworks.com/supportfiles/vision/data/deeplabv3plusResNet50SpaceNetDistance.mat';
    netDir = fullfile(tempdir,'SpaceNet');
    downloadTrainedSpaceNetDistanceNet(trainedSpaceNetDistanceURL,netDir);
    load(fullfile(netDir,'deeplabv3plusResNet50SpaceNetDistance.mat'));
end

%
%% Perform Semantic Segmentation of Test Data
% A GPU is highly recommended for segmentation of the test images (requires Parallel Computing Toolbox™).
% For each test image use apply method of bigimage object to run semantic segmentation. The apply method runs the algorithm on all the blocks in each image. The block-wise results for each test image are assembled back into one output segmented image. The segmented results are overlayed on top of the test image. Only the regions that have the mask set to true are processed. The images must be rescaled so that they display correctly. The mask transparency is set such that the test image is visible underneath the segmented results.
% Display the overlay segmented results for 2 test images using the network trained using the building label images. The function createSegMap evaluates semantic segmentation results for the network. After segmentation, the helper function get pixels corresponding to building by thresholding pixel values greater than the boundary label to . The segmented values are expected to be in the range [1 3]. Any pixel with a value greater than 1 is considered to be the building boundary or within the building.

testImageIndices = numTrain+numVal+1:numTrain+numVal+numTest;
bigTestImages = bigImages(testImageIndices);
bigTestMasks = bigMasks(testImageIndices);
bigTestBuildingsGT = bigBuildingLabels(testImageIndices);
bigTestDistancesGT = bigDistanceLabels(testImageIndices);

boundaryLabelBuilding = 1;
boundaryLabelDistance = 64;

for idx = 1:numTest
    
    bigTestImage = bigTestImages(idx);
    bigMaskImage = bigTestMasks(idx);

    bigSegBuildings(idx) = apply(bigTestImage,1,@(x)createSegMap(x,trainedNetBuilding,boundaryLabelBuilding), ...
        'Mask',bigMaskImage,'InclusionThreshold',0, ...
        'BlockSize',patchSize(1:2),'PadPartialBlocks',true);
    
    bigSegDistances(idx) = apply(bigTestImage,1,@(x)createSegMap(x,trainedNetDistance,boundaryLabelDistance), ...
        'Mask',bigMaskImage,'InclusionThreshold',0, ...
        'BlockSize',patchSize(1:2),'PadPartialBlocks',true);

end

%
% Display the results of the segmentation.
% We notice that the results are not as well defined for some buildings. The borders are a bit fuzzy and bleed into other regions. So let's evaluate the  semantic segmentation results for network trained using distance labels to see if it produces sharper results.
% Now, display the overlay segmented results for 2 test images using the network trained using the distance label images. The function createSegMap performs semantic segmentation using the network trained using distance label images. Then, it  thresholds pixel values greater than 64 to get only pixels corresponding to building. The segmented values are expected to be in the range [1 129]. Any pixel with a distance value greater than 64 is considered to be the building boundary or within the building.
% The class number represents the signed distance of each pixel in the range [-64 64]. We shift the distance to all positive numbers (by adding 65 to all values). Therefore, values 1 to 64 represent outside building pixels, 65 represents bounday pixel and 66 to 129 represents inside building pixels.
% Notice how the results from the network trained using distance label images do a better job of segmenting the building pixels by providing a sharper, accurate boundary.

for idx = 1:2
    
    bigTestImage = getFullLevel(bigTestImages(idx));
    bigSegBuilding = getFullLevel(bigSegBuildings(idx)); 
    bigSegDistance = getFullLevel(bigSegDistances(idx));

    I = rescale(bigTestImage);
    B = rescale(bigSegDistance);
    D = rescale(bigSegBuilding);

    fusedImageBuilding = imfuse(I,B,'blend');
    fusedImageDistance = imfuse(I,D,'blend');

    figure
    montage({fusedImageBuilding,fusedImageDistance},'BorderSize',10,'BackgroundColor','w');
    title(['Test Image ',num2str(idx),': Building-Based Segmentation (Left) Versus Distance-Based Segmentation (Right)'])
end

%
%% Quantify Segmentation Accuracy
% The visual evaluation shows that the network trained using distance label images performs better. Now, compute Jaccard scores for the test images to see how the two networks perform.

jaccardScoreBuilding = zeros(numTest,1);
jaccardScoreDistance = zeros(numTest,1);

for idx = 1:numTest
    
    bigSegBuilding = getFullLevel(bigSegBuildings(idx));
    bigTestBuilding = getFullLevel(bigTestBuildingsGT(idx));
    gtBuilding = double(bigTestBuilding)>boundaryLabelBuilding;
    
    bigSegDistance = getFullLevel(bigSegDistances(idx));
    bigTestDistance =  getFullLevel(bigTestDistancesGT(idx));
    gtDistance = double(bigTestDistance)>boundaryLabelDistance;

    jaccardScoreBuilding(idx) = jaccard(bigSegBuilding,gtBuilding);
    jaccardScoreDistance(idx) = jaccard(bigSegDistance,gtDistance);
    
end

% Display the average Jaccard score for all test images. The network trained using distance labels performs better.

disp(['Average Jaccard score of building-labeled segmentation: ',num2str(mean(jaccardScoreBuilding))])
disp(['Average Jaccard score of distance-labeled segmentation: ',num2str(mean(jaccardScoreDistance))])

%
%% Supporting Functions
%% Evaluate Semantic Segmentation Results
% The function createSegMap evaluates semantic segmentation results for the segmentation network.

function segmentedOutput = createSegMap(im,trainedNet,boundaryLabel)
    im = mat2gray(im);
    segmentedOutput = semanticseg(im,trainedNet);
    segmentedOutput = double(segmentedOutput)>boundaryLabel;
end

%% References
% [1] SpaceNet on Amazon Web Services (AWS). “Datasets.” The SpaceNet Catalog.  Last modified April 30, 2018. Accessed on January 6, 2020. https://spacenetchallenge.github.io/datasets/datasetHomePage.html.
% [2] A Van Etten, D Lindenbaum, TM Bacastow "Spacenet: A remote sensing dataset and challenge series." arXiv preprint arXiv:1807.01232, 2018.
% [3] The SpaceNet Dataset by SpaceNet Partners is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
% [4] ImageNet. http://www.image-net.org.
% https://spacenetchallenge.github.io/Challenges/Challenge-1.html
% https://spacenetchallenge.github.io/Challenges/Challenge-2.html
% https://spacenetchallenge.github.io/AOI_Lists/AOI_2_Vegas.html
% Copyright 2019 The MathWorks, Inc.