Main Content

Hyperspectral Image Analysis Using Maximum Abundance Classification

This example shows how to identify different regions in a hyperspectral image by performing maximum abundance classification (MAC). An abundance map characterizes the distribution of an endmember across a hyperspectral image. Each pixel in the image is either a pure pixel or a mixed pixel. The set of abundance values obtained for each pixel represents the percentage of each endmembers present in that pixel. In this example, you will classify the pixels in a hyperspectral image by finding the maximum abundance value for each pixel and assigning it to the associated endmember class.

This example requires the Hyperspectral Imaging Library for Image Processing Toolbox™. You can install the Hyperspectral Imaging Library for Image Processing Toolbox from Add-On Explorer. For more information about installing add-ons, see Get and Manage Add-Ons. The Hyperspectral Imaging Library for Image Processing Toolbox requires desktop MATLAB®, as MATLAB® Online™ and MATLAB® Mobile™ do not support the library.

This example uses a data sample from the Pavia University dataset as test data. The test data contains nine endmembers that represent these ground truth classes: Asphalt, Meadows, Gravel, Trees, Painted metal sheets, Bare soil, Bitumen, Self blocking bricks, and Shadows.

Load and Visualize Data

Load the .mat file containing the test data into the workspace. The .mat file contains an array paviaU, representing the hyperspectral data cube and a matrix signatures, representing the nine endmember signatures taken from the hyperspectral data. The data cube has 103 spectral bands with wavelengths ranging from 430 nm to 860 nm. The geometric resolution is 1.3 meters and the spatial resolution of each band image is 610-by-340.

load('paviaU.mat');
img = paviaU;
sig = signatures;

Compute the central wavelength for each spectral band by evenly spacing the wavelength range across the number of spectral bands.

wavelengthRange = [430 860];
numBands = 103;
wavelength = linspace(wavelengthRange(1),wavelengthRange(2),numBands);

Create a hypercube object using the hyperspectral data cube and the central wavelengths. Then estimate an RGB image from the hyperspectral data. Set the ContrastStretching parameter value to true in order to improve the contrast of the RGB output. Visualize the RGB image.

hcube = hypercube(img,wavelength);
rgbImg = colorize(hcube,'Method','RGB','ContrastStretching',true);
figure
imshow(rgbImg)

Figure contains an axes object. The hidden axes object contains an object of type image.

The test data contains the endmember signatures of nine ground truth classes. Each column of sig contain the endmember signature of a ground truth class. Create a table that lists the class name for each endmember and the corresponding column of sig.

num = 1:size(sig,2);
endmemberCol = num2str(num');
classNames = {'Asphalt';'Meadows';'Gravel';'Trees';'Painted metal sheets';'Bare soil';...
              'Bitumen';'Self blocking bricks';'Shadows'};
table(endmemberCol,classNames,'VariableName',{'Column of sig';'Endmember Class Name'})
ans=9×2 table
    Column of sig      Endmember Class Name  
    _____________    ________________________

          1          {'Asphalt'             }
          2          {'Meadows'             }
          3          {'Gravel'              }
          4          {'Trees'               }
          5          {'Painted metal sheets'}
          6          {'Bare soil'           }
          7          {'Bitumen'             }
          8          {'Self blocking bricks'}
          9          {'Shadows'             }

Plot the endmember signatures.

figure
plot(sig)
xlabel('Band Number')
ylabel('Data Values')
ylim([400 2700])
title('Endmember Signatures')
legend(classNames,'Location','NorthWest')

Figure contains an axes object. The axes object with title Endmember Signatures, xlabel Band Number, ylabel Data Values contains 9 objects of type line. These objects represent Asphalt, Meadows, Gravel, Trees, Painted metal sheets, Bare soil, Bitumen, Self blocking bricks, Shadows.

Estimate Abundance Maps

Create abundance maps the endmembers by using the estimateAbundanceLS function and select the method as full constrained least squares (FCLS). The function outputs the abundance maps as a 3-D array with the spatial dimensions as the input data. Each channel is the abundance map of the endmember from the corresponding column of signatures. In this example, the spatial dimension of the input data is 610-by-340 and the number of endmembers is 9. So, the size of the output abundance map is 610-by-340-by-9.

abundanceMap = estimateAbundanceLS(hcube,sig,'Method','fcls');

Display the abundance maps.

fig = figure('Position',[0 0 1100 900]);
n = ceil(sqrt(size(abundanceMap,3)));
for cnt = 1:size(abundanceMap,3)
    subplot(n,n,cnt)
    imagesc(abundanceMap(:,:,cnt))
    title(['Abundance of ' classNames{cnt}])
    hold on
end
hold off

Figure contains 9 axes objects. Axes object 1 with title Abundance of Asphalt contains an object of type image. Axes object 2 with title Abundance of Meadows contains an object of type image. Axes object 3 with title Abundance of Gravel contains an object of type image. Axes object 4 with title Abundance of Trees contains an object of type image. Axes object 5 with title Abundance of Painted metal sheets contains an object of type image. Axes object 6 with title Abundance of Bare soil contains an object of type image. Axes object 7 with title Abundance of Bitumen contains an object of type image. Axes object 8 with title Abundance of Self blocking bricks contains an object of type image. Axes object 9 with title Abundance of Shadows contains an object of type image.

Perform Maximum Abundance Classification

Find the channel number of the largest abundance value for each pixel. The channel number returned for each pixel corresponds to the column in sig that contains the endmember signature associated with the maximum abundance value of that pixel. Display a color coded image of the pixels classified by maximum abundance value.

[~,matchIdx] = max(abundanceMap,[],3);
figure
imagesc(matchIdx)
colormap(jet(numel(classNames)))
colorbar('TickLabels',classNames)

Figure contains an axes object. The axes object contains an object of type image.

Segment the classified regions and overlay each of them on the RGB image estimated from the hyperspectral data cube.

segmentImg = zeros(size(matchIdx));
overlayImg = zeros(size(abundanceMap,1),size(abundanceMap,2),3,size(abundanceMap,3));
for i = 1:size(abundanceMap,3)
    segmentImg(matchIdx==i) = 1;
    overlayImg(:,:,:,i) = imoverlay(rgbImg,segmentImg);
    segmentImg = zeros(size(matchIdx));
end

Display the classified and the overlaid hyperspectral image regions along with their class names. From the images, you can see that the asphalt, trees, bare soil, and brick regions have been accurately classified.

figure('Position',[0 0 1100 900]);
n = ceil(sqrt(size(abundanceMap,3)));
for cnt = 1:size(abundanceMap,3)
    subplot(n,n,cnt);
    imagesc(uint8(overlayImg(:,:,:,cnt)));
    title(['Regions Classified as ' classNames{cnt}])
    hold on
end
hold off

Figure contains 9 axes objects. Axes object 1 with title Regions Classified as Asphalt contains an object of type image. Axes object 2 with title Regions Classified as Meadows contains an object of type image. Axes object 3 with title Regions Classified as Gravel contains an object of type image. Axes object 4 with title Regions Classified as Trees contains an object of type image. Axes object 5 with title Regions Classified as Painted metal sheets contains an object of type image. Axes object 6 with title Regions Classified as Bare soil contains an object of type image. Axes object 7 with title Regions Classified as Bitumen contains an object of type image. Axes object 8 with title Regions Classified as Self blocking bricks contains an object of type image. Axes object 9 with title Regions Classified as Shadows contains an object of type image.

See Also

| |

Related Topics