%% Segment Vegetation Regions in Hyperspectral Data Using Spectral Matching
% The spectral matching method compares the spectral signature of each pixel in the hyperspectral data cube with a reference spectral signature for vegetation from an ECOSTRESS spectrum file.
% Read the spectral signature of vegetation from the ECOSTRESS spectral library.

fileroot = matlabshared.supportpkg.getSupportPackageRoot();
filename = fullfile(fileroot,'toolbox','images','supportpackages','hyperspectral','hyperdata',...
               'ECOSTRESSSpectraFiles','vegetation.tree.tsuga.canadensis.vswir.tsca-1-47.ucsb.asd.spectrum.txt');
libData = readEcostressSig(filename);

% Read the hyperspectral data into the workspace.

hcube = hypercube('paviaU.hdr');

% Compute the distance scores of the spectra of the hyperspectral data pixels with respect to the reference spectra.

score = spectralMatch(libData,hcube);

% Display the distance scores. The pixels with low distance scores are stronger matches to the reference spectra and are more likely to belong to the vegetation region.

figure
imagesc(score)
colorbar

% Define a threshold for detecting distance scores that correspond to the vegetation region.

threshold = 0.3;

% Generate a binary image by assigning a intensity value 1 for pixels with score less than a specified threshold. Other regions are assigned the intensity value 0. The maximum intensity regions in the binary image correspond to the vegetation regions in the hyperspectral data cube.

bw = score < threshold;

% Segment the vegetation regions of the hyperspectral data cube by using the indices of the maximum intensity regions in the binary image.

T = reshape(hcube.DataCube,[size(hcube.DataCube,1)*size(hcube.DataCube,2) size(hcube.DataCube,3)]);
Ts = zeros(size(T));
Ts(bw == 1,:) = T( bw==1 ,:);
Ts = reshape(Ts,[size(hcube.DataCube,1) size(hcube.DataCube,2) size(hcube.DataCube,3)]);

% . Create a new hypercube object that contains only the segmented vegetation regions.

segmentedDataCube = hypercube(Ts,hcube.Wavelength);

% Estimate the RGB colour image of the original data cube and the segmented data cube by using the colorize function.

rgbImg = colorize(hcube,'Method','rgb','ContrastStretching',true);
segmentedImg = colorize(segmentedDataCube,'Method','rgb','ContrastStretching',true);

% Overlay the binary image on the RGB version of the original data cube by using the imoverlay function.

B = imoverlay(rgbImg,bw,'Yellow');

% Display the RGB colour images of the original data cube and the segmented data cube along with the overlaid image. The segmented image contains only the vegetation regions that are segmented from the original data cube.

figure
montage({rgbImg segmentedImg B},'Size',[1 3])
title(['Original Image | ' 'Segmented Image | ' 'Overlayed Image'])

% Copyright 2020 The MathWorks, Inc.