{"id":10597,"date":"2023-11-29T08:57:11","date_gmt":"2023-11-29T13:57:11","guid":{"rendered":"https:\/\/blogs.mathworks.com\/student-lounge\/?p=10597"},"modified":"2023-11-29T09:23:58","modified_gmt":"2023-11-29T14:23:58","slug":"kelp-wanted-challenge-starter-code","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/student-lounge\/2023\/11\/29\/kelp-wanted-challenge-starter-code\/","title":{"rendered":"Kelp Wanted Challenge Starter Code"},"content":{"rendered":"
\n

Getting Started with MATLAB<\/h2>\n
We at MathWorks<\/a>, in collaboration with DrivenData<\/a>, are excited to bring you this challenge<\/a>! The goal is to develop an algorithm that can use provided satellite imagery to predict where kelp is present and where it is not. Kelp is a type of seaweed or algae that often grows in clusters known as kelp forests<\/a>, which provide shelter and stability for many coastal ecosystems. The presence and growth of kelp is an important measurement for evaluating the health of these ecosystems, so the ability to easily monitor kelp forests could be a huge step forward in coastal climate science. In this blog, we will explore the data using the Hyperspectral Viewer app, preprocess the dataset, then create, evaluate, and use a basic<\/span> semantic segmentation model to solve this challenge.<\/div>\n
To request your complimentary MATLAB license<\/span> and access additional learning resources, check out this website<\/a>!<\/div>\n

Table of Contents:<\/h3>\n
    \n
  1. Explore and Understand the Data<\/a><\/li>\n
  2. Import the Data<\/a><\/li>\n
  3. Preprocess the Data<\/a><\/li>\n
  4. Design and Train a Neural Network<\/a><\/li>\n
  5. Evaluate the Model<\/a><\/li>\n
  6. Create Submissions<\/a><\/li>\n<\/ol>\n

    Explore and Understand the Data<\/h2>\n
    Instructions for accessing and downloading the competition data can be found here<\/a>.<\/div>\n

    The Input: Satellite Images<\/h3>\n
    The input data is a set of augmented satellite images that have seven layers or “bands”, so you can think of it as 7 separate images all stacked on top of each other, as shown below<\/div>\n
    \"\"<\/div>\n
    Each band is looking at the same exact patch of earth, but they each contain different measurements. The first 5 bands contain measurements taken at different wavelengths of the light spectrum, and the last two are supplementary metrics to better understand the environment. The following list shows what each of the seven bands measures:<\/div>\n
      \n
    1. Short-wave infrared (SWIR)<\/li>\n
    2. Near infrared (NIR)<\/li>\n
    3. Red<\/li>\n
    4. Green<\/li>\n
    5. Blue<\/li>\n
    6. Cloud Mask (binary – is there cloud or not)<\/li>\n
    7. Digital Elevation Model (meters above sea-level)<\/li>\n<\/ol>\n
      Typically, most standard images just measure the red, green, and blue values, but by including additional measurements, hyperspectral images can enable us to identify objects and patterns that may not be easily seen with the naked eye, such as underwater kelp.<\/div>\n
      \"kelp-example.png\"<\/div>\n
      [Left: A true color image of an example tile using the RGB bands. Center: A false color image using the SWIR, NIR, and Red bands. Right: The false color image with the labeled kelp mask overlayed in cyan.]<\/div>\n
      Let’s read in a sample image and label for tile ID AA498489, which we will explore to gain a better understanding of the data.<\/div>\n
      \n
      \n
      firstImage = imread(‘train_features\/AA498489_satellite.tif’<\/span>);<\/span><\/div>\n<\/div>\n
      \n
      firstLabel = imread(‘train_labels\/AA498489_kelp.tif’<\/span>);<\/span><\/div>\n<\/div>\n<\/div>\n

      The Spectral Bands (1-5)<\/h4>\n
      Let’s start by exploring the first five layers. The rescale<\/span><\/a> function adjusts the values of the bands so that they can be visualized as grayscale images, and the montage<\/span><\/a> function displays each band next to each other.<\/div>\n
      \n
      \n
      montage(rescale(firstImage(:, :, 1:5)));<\/span><\/div>\n
      \n
      \n
      <\/div>\n
      <\/div>\n
      <\/div>\n
      <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
      Here we can see that there’s some land masses present, and that the SWIR and NIR bands have higher values than the red, green, and blue bands when looking at this patch of earth, as they are brighter. This doesn’t tell us much about the data, but gives us an idea of what we are looking at.<\/div>\n

      Hyperspectral Viewer<\/h4>\n
      You can use the Hyperspectral Viewer app<\/a> to further explore the first five layers. Note that this requires the Image Processing Toolbox\u2122 Hyperspectral Imaging Library<\/a>, which can be installed through the Add-On Explorer<\/a>, and is not supported in MATLAB Online. The center wavelengths shown below are approximated from this resource<\/a>.<\/div>\n
      \n
      \n
      firstImSatellite = firstImage(:, :, 1:5);<\/span><\/div>\n<\/div>\n
      \n
      centerWavelengths = [1650, 860, 650, 550, 470]; % in nanometers<\/span><\/span><\/div>\n<\/div>\n
      \n
      hcube = hypercube(firstImSatellite, centerWavelengths);<\/span><\/div>\n<\/div>\n
      \n
      hyperspectralViewer(hcube);<\/span><\/div>\n<\/div>\n<\/div>\n
      When the app opens, you’ll have the ability to view single bands on the left pane and various band combinations on the right. Note that the bands are shown in order of wavelength, not<\/span> in the order they are loaded, so in the app the bands are in reverse order. Band 1 = Blue, Band 5 = SWIR.<\/div>\n
      On the left pane, you can scroll through and view each band one at a time. You can also manually adjust the contrast to make it easier to see or to make it representative of a different spectrum than the default.<\/div>\n
      \"ExploreBands.gif\"<\/div>\n
      On the right, you’ll have the ability to see False Color, RGB, and CIR images. RGB images are just standard color images, and show the earth as we would see it from a typical camera. False Color and CIR images convert the measurements from the SWIR and NIR bands, which are not visible from the human eye, to colors that we can see. You can manually adjust the bands to create custom images as well.<\/div>\n
      In this pane, you also have the ability to create spectral plots for a single pixel, which shows what value that pixel holds for each band. Since this image has land, sea, and coast, I’ll create spectral plots for a pixel in each of these areas to see how they differ.<\/div>\n
      \"BandCombos.gif\"<\/div>\n
      This app also provides the ability to plot and interact with various spectral indices that calculate different measurements related to vegetation, which could provide helpful additional information when looking for kelp. Learn more about these spectral indices by checking out this documentation link<\/a>.<\/div>\n
      \"Indices.gif\"<\/div>\n
      If you have some plots that you’d like to work with further, you can export any of these to the MATLAB workspace. I’ll use the RGB image in a moment, so let’s export it.<\/div>\n
      \"Export.gif\"<\/div>\n

      The Physical Property Bands<\/h4>\n
      The other two layers of the input images are not based on the light spectrum, but on physical properties. The cloud mask can be visualized as a black-and-white image, where black means there was no cloud present and white means there was cloud blocking that part of the image.<\/div>\n
      \n
      \n
      cloudMask = firstImage(:, :, 6);<\/span><\/div>\n<\/div>\n
      \n
      imshow(double(cloudMask));<\/span><\/div>\n
      \n
      \n
      <\/div>\n
      <\/div>\n
      <\/div>\n
      <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
      This image is almost all black, so there was very little cloud blocking the satellite, but there are a few white pixels as highlighted in the image below.<\/div>\n
      \"\"<\/div>\n
      The elevation mask can be visualized using the imagesc<\/span><\/a> function, which will colorize different parts of the image based on how high above sea level each pixel is. As one might expect, the highest elevation in our image correlates to the large land mass.<\/div>\n
      \n
      \n
      elevationModel = firstImage(:, :, 7);<\/span><\/div>\n<\/div>\n
      \n
      imagesc(elevationModel);<\/span><\/div>\n<\/div>\n
      \n
      colormap(turbo);<\/span><\/div>\n<\/div>\n
      \n
      colorbar;<\/span><\/div>\n
      \n
      \n
      <\/div>\n
      <\/div>\n
      <\/div>\n
      <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n

      The Output: A Binary Mask<\/h3>\n
      The corresponding label for this satellite image is a binary mask, similar to the cloud mask. It is 350×350 – the same height and width of the satellite images – and each pixel is labeled as either 1 (kelp detected) or 0 (no kelp detected).<\/div>\n
      \n
      \n
      imshow(double(firstLabel))<\/span><\/div>\n
      \n
      \n
      <\/div>\n
      <\/div>\n
      <\/div>\n
      <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
      You can add these labels over the RGB satellite image we exported earlier to see where the kelp is in relation to the land masses.<\/div>\n
      \n
      \n
      labeledIm = labeloverlay(rgb, firstLabel);<\/span><\/div>\n<\/div>\n
      \n
      imshow(labeledIm);<\/span><\/div>\n
      \n
      \n
      <\/div>\n
      <\/div>\n
      <\/div>\n
      <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n

      Import the Data<\/h2>\n
      To start working with all of the data in MATLAB, you can use an imageDatastore<\/span><\/a> and pixelLabelDatastore<\/span><\/a>. pixelLabelDatastore<\/span> expects uint8 data, but the labels are currenlty int8, so I’ve created a custom read function (readLabelData<\/span><\/a>) to convert the label data to the correct format.<\/div>\n
      \n
      \n
      trainImagesPath = ‘.\/train_features’<\/span>;<\/span><\/div>\n<\/div>\n
      \n
      trainLabelsPath = ‘.\/train_labels’<\/span>;<\/span><\/div>\n<\/div>\n
      \n
      <\/div>\n<\/div>\n
      \n
      allTrainIms = imageDatastore(trainImagesPath);<\/span><\/div>\n<\/div>\n
      \n
      <\/div>\n<\/div>\n
      \n
      classNames = [“nokelp”<\/span>, “kelp”<\/span>];<\/span><\/div>\n<\/div>\n
      \n
      pixelLabelIDs = [0, 1];<\/span><\/div>\n<\/div>\n
      \n
      allTrainLabels = pixelLabelDatastore(trainLabelsPath, classNames, pixelLabelIDs, ReadFcn=@readLabelData);<\/span><\/div>\n<\/div>\n<\/div>\n
      Now we can divide the data into a training, validation, and testing data. The training set will be used to train our model, the validation set will be used to check in on training and make sure the model is not overfitting, and the testing set will be used after the model is trained to see how well it generalizes to new data.<\/div>\n
      \n
      \n
      numObservations = numel(allTrainIms.Files);<\/span><\/div>\n<\/div>\n
      \n
      <\/div>\n<\/div>\n
      \n
      numTrain = round(0.7 * numObservations);<\/span><\/div>\n<\/div>\n
      \n
      numVal = round(0.15 * numObservations);<\/span><\/div>\n<\/div>\n
      \n
      <\/div>\n<\/div>\n
      \n
      trainIms = subset(allTrainIms, 1:numTrain);<\/span><\/div>\n<\/div>\n
      \n
      trainLabels = subset(allTrainLabels, 1:numTrain);<\/span><\/div>\n<\/div>\n
      \n
      <\/div>\n<\/div>\n
      \n
      valIms = subset(allTrainIms, (numTrain + 1):(numTrain + numVal));<\/span><\/div>\n<\/div>\n
      \n
      valLabels = subset(allTrainLabels, (numTrain + 1):(numTrain + numVal));<\/span><\/div>\n<\/div>\n
      \n
      <\/div>\n<\/div>\n
      \n
      testIms = subset(allTrainIms, (numTrain + numVal + 1):numObservations);<\/span><\/div>\n<\/div>\n
      \n
      testLabels = subset(allTrainLabels, (numTrain + numVal + 1):numObservations);<\/span><\/div>\n<\/div>\n<\/div>\n

      Preprocess The Data<\/h2>\n

      Clean up the sample image<\/h3>\n
      Now that we have a better understanding of our data, we can preprocess it! In this section, I will show some ways you can:<\/div>\n
        \n
      1. Resize the data<\/li>\n
      2. Normalize the data<\/li>\n
      3. Augment the data<\/li>\n<\/ol>\n
        While ideally each image in the dataset will be the same size, data is messy, and this isn’t always the case. I’ll use imresize<\/span><\/a> to ensure the height and width of each image is correct.<\/div>\n
        \n
        \n
        inputSize = [350 350 8];<\/span><\/div>\n<\/div>\n
        \n
        firstImage = imresize(firstImage, inputSize(1:2));<\/span><\/div>\n<\/div>\n<\/div>\n
        Each band has a different minimum and maximum, so while a 1 may be low for some bands it could be a high value for other bands. Let’s go through each layer (except for the cloud mask) and rescale it so that the minimum values are 0 and the maximum values are 1. There are many ways to normalize your data, so I suggest testing out other algorithms.<\/div>\n
        \n
        \n
        normalizedImage = zeros(inputSize); % preallocate for speed<\/span><\/span><\/div>\n<\/div>\n
        \n
        continuousBands = [1 2 3 4 5 7];<\/span><\/div>\n<\/div>\n
        \n
        for <\/span>band = continuousBands<\/span><\/div>\n<\/div>\n
        \n
        normalizedImage(:, :, band) = rescale(firstImage(:, :, band));<\/span><\/div>\n<\/div>\n
        \n
        end<\/span><\/span><\/div>\n<\/div>\n
        \n
        normalizedImage(:, :, 6) = firstImage(:, :, 6);<\/span><\/div>\n<\/div>\n<\/div>\n
        You can also use the provided data to create more data! This is called feature extraction<\/span>. Since I know that kelp is often found along coasts, I’ll use an edge detection algorithm to show the edges that exist in the image, which will often include coastlines.<\/div>\n
        \n
        \n
        normalizedImage(:, :, 8) = edge(firstImage(:, :, 4), “sobel”<\/span>);<\/span><\/div>\n<\/div>\n<\/div>\n
        Now we can view our preprocessed data!<\/div>\n
        \n
        \n
        montage(normalizedImage)<\/span><\/div>\n
        \n
        \n
        <\/div>\n
        <\/div>\n
        <\/div>\n
        <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n

        Apply Preprocessing to the Entire Dataset<\/h3>\n
        To make sure these preprocessing steps are applied to every image in the dataset, you can use the transform<\/span><\/a> function. This allows you to apply a function of your choice to each image as it is read, so I have defined a function cleanSatelliteData<\/span><\/a> (shown at the end of the blog) that applies these steps to every image.<\/div>\n
        \n
        \n
        trainImsProcessed = transform(trainIms, @cleanSatelliteData);<\/span><\/div>\n<\/div>\n
        \n
        valImsProcessed = transform(valIms, @cleanSatelliteData);<\/span><\/div>\n<\/div>\n<\/div>\n
        Then we combine the input and output datastores so that each satellite image can easily be associated with it’s expected output.<\/div>\n
        \n
        \n
        trainData = combine(trainImsProcessed, trainLabels);<\/span><\/div>\n<\/div>\n
        \n
        valData = combine(valImsProcessed, valLabels);<\/span><\/div>\n<\/div>\n<\/div>\n
        If you preview the resulting datastore, the satellite images are now 350x350x8 instead of 350x350x7 since we added a band in the transformation function.<\/div>\n
        \n
        \n
        firstSample = preview(trainData)<\/span><\/div>\n
        \n
        \n
        \n
        firstSample = 1\u00d72 cell <\/span><\/div>\n
        \n\n\n\n\n\n
        <\/th>\n1<\/th>\n2<\/th>\n<\/tr>\n<\/thead>\n
        1<\/th>\n350\u00d7350\u00d78\u00a0double<\/td>\n350\u00d7350\u00a0categorical<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n
        <\/div>\n
        <\/div>\n
        <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n

        Design and Train a Neural Network<\/h2>\n

        Create the network layers<\/h3>\n
        Once the data is ready, it’s time to create a neural network.I’m going to create a simple network for semantic segmentation using the segnetLayers<\/span><\/a> function.<\/div>\n
        \n
        \n
        numClasses = 2;<\/span><\/div>\n<\/div>\n
        \n
        lgraph = segnetLayers(inputSize, numClasses, 5);<\/span><\/div>\n<\/div>\n<\/div>\n

        Balance the Classes<\/h3>\n
        In the sample “firstImage”, there were a lot of pixels with the 0 label, meaning no kelp was detected. Ideally, we would have equal amounts of “kelp” and “nokelp” labels so that the network would learn each equally, but most images probably don’t show 50% or more kelp. To see the exact distribution of class labels in the dataset, use countEachLabel<\/span><\/a>, which counts the number of pixels by class label.<\/div>\n
        \n
        \n
        labelCounts = countEachLabel(trainLabels)<\/span><\/div>\n
        \n
        \n
        \n
        labelCounts = 2\u00d73 table <\/span><\/div>\n
        \n\n\n\n\n\n\n
        <\/th>\nName<\/th>\nPixelCount<\/th>\nImagePixelCount<\/th>\n<\/tr>\n<\/thead>\n
        1<\/th>\n‘nokelp’<\/td>\nundefined<\/td>\nundefined<\/td>\n<\/tr>\n
        2<\/th>\n‘kelp’<\/td>\nundefined<\/td>\nundefined<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n
        <\/div>\n
        <\/div>\n
        <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        ‘PixelCount’ shows how many total pixels contained that class, and ‘ImagePixelCount’ shows the total number of pixels in all images that contained that class. This shows that not only are there way <\/span>more “nokelp” labels than “kelp” labels, but also that there are images that don’t contain any “kelp” labels. If not handled correctly, this imbalance can be detrimental to the learning process because the learning is biased in favor of “nokelp”. To improve training, you can use class weights to balance the classes<\/a>. Class weights define the relative importance of each class to the training process, and by default is set to 1 for each class. By assigning class weights that are inversely proportional to the frequency of each class (i.e., giving the “kelp” class a higher weight than “nokelp”), we reduce the chance of the network having a strong bias towards more common classes.<\/div>\n
        Use the pixel label counts from above to calculate the median frequency class weights:<\/div>\n
        \n
        \n
        imageFreq = labelCounts.PixelCount .\/ labelCounts.ImagePixelCount;<\/span><\/div>\n<\/div>\n
        \n
        classWeights = median(imageFreq) .\/ imageFreq<\/span><\/div>\n<\/div>\n<\/div>\n
        You can then pass the class weights to the network by creating a new pixelClassificationLayer<\/span><\/a> and replacing the default one.<\/div>\n
        \n
        \n
        pxLayer = pixelClassificationLayer(‘Name’<\/span>,‘labels’<\/span>,‘Classes’<\/span>,labelCounts.Name,‘ClassWeights’<\/span>,classWeights);<\/span><\/div>\n<\/div>\n
        \n
        lgraph = replaceLayer(lgraph,“pixelLabels”<\/span>,pxLayer);<\/span><\/div>\n<\/div>\n<\/div>\n

        Train the Network<\/h3>\n
        Specify the settings you want to use for training with the trainingOptions<\/span><\/a> function, and train the network!<\/div>\n
        \n
        \n
        tOps = trainingOptions(“sgdm”<\/span>, InitialLearnRate=0.001, …<\/span><\/span><\/div>\n<\/div>\n
        \n
        MiniBatchSize=32, …<\/span><\/span><\/div>\n<\/div>\n
        \n
        MaxEpochs=5, …<\/span><\/span><\/div>\n<\/div>\n
        \n
        ValidationData=valData);<\/span><\/div>\n<\/div>\n
        \n
        <\/div>\n<\/div>\n
        \n
        trainedNet = trainNetwork(trainData, lgraph, tOps);<\/span><\/div>\n<\/div>\n<\/div>\n
        \"\"<\/div>\n
        This is an example of training a neural network from the command line, but if you want to explore your neural networks visually or go through the deep learning steps interactively, check out the Deep Network Designer app documentation<\/a> and starter video<\/a>!<\/div>\n

        Evaluate the Model<\/h2>\n
        To test the quality of your model before submission, you need to process your testing data (which we created earlier<\/a>) the same way you processed your training data<\/div>\n
        \n
        \n
        testIms = transform(testIms, @cleanSatelliteData);<\/span><\/div>\n<\/div>\n<\/div>\n
        We need to create a folder to contain the predictions<\/div>\n
        \n
        \n
        if <\/span>~exist(‘evaluationTest’<\/span>, ‘dir’<\/span>)<\/span><\/div>\n<\/div>\n
        \n
        mkdir evaluationTest<\/span>;<\/span><\/div>\n<\/div>\n
        \n
        end<\/span><\/span><\/div>\n<\/div>\n<\/div>\n
        Then we make predictions on the test data!<\/div>\n
        \n
        \n
        allPreds = semanticseg(testIms,trainedNet,…<\/span><\/span><\/div>\n<\/div>\n
        \n
        MiniBatchSize=32, …<\/span><\/span><\/div>\n<\/div>\n
        \n
        WriteLocation=“evaluationTest”<\/span>);<\/span><\/div>\n
        \n
        \n
        Running semantic segmentation network
        \n————————————-
        \n* Processed 846 images.<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        Once we have a set of predictions, we can use the evaluateSemanticSegmentation<\/span><\/a> function to compare the predictions with the actual labels and get a sense of how well the model will perform on new data.<\/div>\n
        \n
        \n
        metrics = evaluateSemanticSegmentation(allPreds,testLabels);<\/span><\/div>\n
        \n
        \n
        Evaluating semantic segmentation results
        \n—————————————-
        \n* Selected metrics: global accuracy, class accuracy, IoU, weighted IoU, BF score.
        \n* Processed 846 images.
        \n* Finalizing… Done.
        \n* Data set metrics:GlobalAccuracy<\/strong> MeanAccuracy<\/strong> MeanIoU<\/strong> WeightedIoU<\/strong> MeanBFScore<\/strong>
        \n______________<\/strong> ____________<\/strong> _______<\/strong> ___________<\/strong> ___________<\/strong>0.94677 0.52232 0.47932 0.94021 0.15665<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        To understand how often the network predicted each class correctly and incorrectly, we can extract the confusion matrix. In a confusion matrix:<\/div>\n