%% Detect Anomalous Pixels in Hyperspectral Data Using RX Detector
% Detect anomalous pixels in hyperspectral data by computing the RX score for each pixel in a hyperspectral data cube. Then compute the threshold for detecting true anomalous pixels by using cumulative probability distribution of RX score values.
% Read hyperspectral data containing anomalous pixels into the workspace.

hcube = hypercube('indian_pines.dat');

% Find anomalous pixels in the input hyperspectral data by using the RX detector. The detector searches for pixels with a high intensity difference within a homogeneous region.

rxScore = anomalyRX(hcube);

% Reduce the dynamic range of the RX score values by rescaling them to the range [0, 255].

rxScore = im2uint8(rescale(rxScore));

%  Display the RX score map. Pixels with a high RX score are likely anomalous pixels.

figure
imagesc(rxScore)
colorbar

% Compute and plot the cumulative probability distribution of RX score values.

count = imhist(rxScore);
pdf = count/prod(size(rxScore,[1 2]));
cdf = cumsum(pdf(:));
figure
plot(cdf)
xlabel('RX Score')
ylabel('Cumulative Probability Values')

% Set the confidence coefficient value to 0.998. Select the first RX score with cumulative probability distribution value greater than the confidence coefficient as the threshold. This threshold represents the RX score above which a pixel is an anomaly with 99.8 percent confidence.

confCoefficient = 0.998;
rxThreshold = find(cdf > confCoefficient,1);

% Apply thresholding to detect anomalous pixels with RX score greater than the computed threshold. The result is a binary image in which the anomalous pixels are assigned the intensity value 1 and other pixels are assigned 0.

bw = rxScore > rxThreshold;

%  Derive the RGB version of the data cube by using the colorize function. Overlay the binary image of anomalous pixels on the RGB image.

rgbImg  = colorize(hcube,'Method','rgb');
B = imoverlay(rgbImg,bw);

% Display both the binary image and the overlaid image.

fig = figure('Position',[0 0 800 400]);
axes1 = axes('Parent',fig,'Position',[0 0.1 0.5 0.8]);
imagesc(bw,'Parent',axes1);
title('Detected Anomalous Pixels')
axis off
colormap gray
axes2 = axes('Parent',fig,'Position',[0.5 0.1 0.5 0.8]);
imagesc(B,'Parent',axes2)
title('Overlaid Image');
axis off  

% Copyright 2020 The MathWorks, Inc.