%% Measure Vegetation Cover in Hyperspectral Data Using NDVI Image
% Read hyperspectral data into the workspace.

hcube = hypercube('indian_pines.dat');

% Compute the NDVI value for each pixel in the data cube.

ndviImg = ndvi(hcube);

% Estimate a contrast-stretched RGB image from the original data cube by using the colorize function.

rgbImg = colorize(hcube,'Method','RGB','ContrastStretching',true);

% Display the original and the NDVI image.

fig = figure('Position',[0 0 1200 600]);
axes1 = axes('Parent',fig,'Position',[0 0.1 0.4 0.8]);
imshow(rgbImg,'Parent',axes1)
title('RGB Image of Data Cube')
axes2 = axes('Parent',fig,'Position',[0.45 0.1 0.4 0.8]);
imagesc(ndviImg,'Parent',axes2)
colorbar
title('NDVI Image')

% Vegetation regions typically have NDVI values from 0.2 and 0.8. NDVI values less than or equal to 0.2 indicate the absence of vegetation. Perform thresholding of NDVI image to segment the vegetation regions. Specify the threshold value.

threshold = 0.2;

% Generate a binary image with a intensity value 1 for pixels with a score greater than or equal to the specified threshold. All other pixels have a value 0. The regions in the binary image with a value of 1 correspond to the vegetation regions in the data cube with NDVI values greater than the threshold.

bw = ndviImg > threshold;

% Overlay the binary image on to the RGB image and display the overlaid image.

overlayImg = imoverlay(rgbImg,bw,[0 1 0]);
figure
imagesc(overlayImg)
title('Vegetation Region Overlaid on RGB Image')

% Compute the vegetation cover based on the total number of pixels in a spectral band and the number of pixels with an NDVI value greater than 0.2.

numVeg = find(bw == 1);
imgSize = size(hcube.DataCube,1)*size(hcube.DataCube,2);
vegetationCover = length(numVeg)/imgSize

% Copyright 2020 The MathWorks, Inc.