%% Identify Vegetation Regions using NDVI and Interactive Threshold Selection
% Normalized difference vegetation index(NDVI) is a one of the most widely used vegetation index used to determine the density of the vegetation. NDVI is computed from near-infrared and visible light reflected by vegetation. Healthy vegetation absorbs more visible light than near-infrared light and Unhealthy or sparse vegetation reflects more visible light than infrared light.
% NDVI image can be calculated as
% \textrm{NDVI}=\frac{\left(\textrm{NIR}-R\right)}{\left(\textrm{NIR}+R\right)}
% Where NIR is a near-infrared channel and R is Red channel. NDVI pixel value ranges from -1 to 1. NDVI values closer to 1 is more likely to be healthy vegetation, close to 0 is likely to be unhealthy vegetation and below 0 is probably not at all vegetation.
% This example shows how to interactively change threshold value for NDVI image. Interactive Normalized Difference Vegetation Index allows you change threshold value for NDVI image to find the vegetation region in hyperspectral datacube.
% Finding vegetation in hyperspectral datacube comprises three steps:
% Step 1: Read Hyperspectral datacube
% Step 2: Compute NDVI image
% Step 3: Apply a simple threshold to get a logical map for vegetation region
% Read hyperspectral datacube into workspace
% Read hyperspectral data file in ENVI format. Here we used pavia dataset which contains nine classes asphalt, meadows, gravel, trees, painted metal sheets, bare soil, bitumen, self-blocking bricks and shadows. NDVI can be used to detect vegetation in it.

hcube = hypercube('paviaU.dat','paviaU.hdr');

% Compute NDVI image
% Function ndvi finds NIR and R bands present in hyperspectral data and computes NDVI image.

ndviImage = ndvi(hcube);

% Set initial threshold value  and overlay the thresholded NDVI image on RGB image

initialThreshold = 0.5;
ndviThreshIm = ndviImage > initialThreshold;
rgb = colorize(hcube, 'method','rgb','ContrastStretching',true);
overlayImage = imoverlay(rgb, ndviThreshIm, 'green');

% Create an Interactive Figure Window
% Create an interactive figure window to display the NDVI image and overlaid RGB image side by side. In the window, you can change the threshold value and find the vegetation region corresponding to the threshold.

f = figure;

% Create a panel in the current figure to interactively adjust the threshold values.
dataPanel = uipanel(f,'Position',[0.01 0.01 0.99 0.2],'Title',...
    'Adjust Threshold Value','FontSize',10);

% Add an user control interface for adjusting the slider.
% Set the default slider value to 0.5.
b = uicontrol(dataPanel,'Style','slider','Position',[75,35,419,23],...
    'value',0.5, 'max',1,'SliderStep',[0.1 0.2]);

% Add an user control interface for displaying the threshold values.
uicontrol(dataPanel,'Style','text','Position',[500,33,45,23], ...
    'String',[num2str(0.5),'/','1']);
uicontrol(dataPanel,'Style','text','Position',[240,5,100,23], ...
    'String','NDVI Index');

% Create a panel in the current figure to display the image.
viewPanel = uipanel(f,'Position',[0.01 0.25 0.99 0.75], ...
    'Title','Interactive Normalized Difference Vegetation Index','FontSize',10);
ax = axes('Parent',viewPanel);

% Display the NDVI image and overlaid image side by side in the interactive figure window.

montage({ndviImage, overlayImage})
title('NDVI Image & Overlay Image' )

% Change the threshold value interactively
% Change the threshold value interactively and it dynamically updates the logical map overlaid on RGB image using callback function ndviThreshold. Assign a function handle that references the  ndviThreshold function to the Callback property to the control object.

b.Callback = @(es,ed) ndviThreshold(b,dataPanel,ndviImage,rgb);

%
% Create Callback Function to update logical mask according to threshold

function ndviThreshold(hObject,dataPanel,ndviImage,rgb)
th = round(hObject.Value,2);

% Apply simple thresholding to get logical map
ndviThreshIm = ndviImage > th;

% Add an user control interface for displaying the threshold value.
uicontrol(dataPanel,'Style','text','Position',[500,33,45,23], ...
    'String',[num2str(th),'/','1']);

% Overlay the thresholded ndvi image into rgb image. 
overlayImage = imoverlay(rgb, ndviThreshIm, 'green');
montage({ndviImage, overlayImage});
title('NDVI Image & Overlay Image' )
end

% Changing the threshold interactively updates the logical map. Higher threshold shows the logical map with healthy vegetation. Lower threshold shows the logical map with healthy vegetation as well as sparse vegetation.
%%
%% References
% [1] Rouse, J. W., Haas, R. H., Schell, J. A., and Deering, D. W., Monitoring vegetation systems in the great plains with ERTS, in Third ERTS Symposium, NASA SP-351, NASA, Washington, DC, Vol. 1, p. 309-317. 1973.
% [2] Haboudane, D., Miller, J.R., Pattey, E., Zarco-Tejada, P.J. and Strachan, I.B., 2004. Hyperspectral vegetation indices and novel algorithms for predicting green LAI of crop canopies: Modeling and validation in the context of precision agriculture. Remote sensing of environment, 90(3), pp.337-352.
% Copyright 2020 The MathWorks, Inc.