スペクトルライブラリを使用したエンドメンバー物質の同定

この例では次を使用します。

この例では、ハイパースペクトルイメージ内に存在するエンドメンバー物質のクラスを同定する方法を示します。エンドメンバーは、1 つの表面物質に属するピクセルの反射率特性を示すピュアスペクトルシグネチャです。既存のエンドメンバー抽出または同定アルゴリズムでは、ハイパースペクトルイメージ内のピュアピクセルを抽出または同定します。ただし、これらの手法でエンドメンバースペクトルが属する物質名またはクラスは同定されません。この例では、エンドメンバーシグネチャを抽出してから、スペクトルマッチングを使用してハイパースペクトルイメージ内のエンドメンバー物質のクラスを分類または同定します。

この例では、Image Processing Toolbox™ Hyperspectral Imaging Library が必要となります。Image Processing Toolbox Hyperspectral Imaging Library はアドオンエクスプローラーからインストールできます。アドオンのインストールの詳細については、アドオンの取得と管理を参照してください。Image Processing Toolbox Hyperspectral Imaging Library は MATLAB® Online™ および MATLAB® Mobile™ ではサポートされないため、デスクトップの MATLAB® が必要となります。

この例では、エンドメンバー物質の同定のため、1) ECOSTRESS スペクトルライブラリのスペクトルシグネチャを基準スペクトルとして使用し、2) Jasper Ridge データセットのデータサンプルをテストデータとして使用します。

ECOSTRESS スペクトルライブラリからの基準データの読み取り

ECOSTRESS ライブラリファイルが格納されているファイルの絶対パスを追加して、ライブラリから読み取るファイルの名前を指定します。

fileroot = matlabshared.supportpkg.getSupportPackageRoot();
addpath(fullfile(fileroot,'toolbox','images','supportpackages','hyperspectral','hyperdata','ECOSTRESSSpectraFiles'));
filenames = ["water.seawater.none.liquid.tir.seafoam.jhu.becknic.spectrum.txt",...
             "water.tapwater.none.liquid.all.tapwater.jhu.becknic.spectrum.txt",...
             "water.ice.none.solid.all.ice_dat_.jhu.becknic.spectrum.txt",...
             "vegetation.tree.eucalyptus.maculata.vswir.jpl087.jpl.asd.spectrum.txt",...
             "soil.utisol.hapludult.none.all.87p707.jhu.becknic.spectrum.txt",...
             "soil.mollisol.cryoboroll.none.all.85p4663.jhu.becknic.spectrum.txt",...
             "manmade.road.tar.solid.all.0099uuutar.jhu.becknic.spectrum.txt",...   
             "manmade.concrete.pavingconcrete.solid.all.0092uuu_cnc.jhu.becknic.spectrum.txt"];
lib = readEcostressSig(filenames);

lib データを表示してその値を検査します。このデータは、クラス、サブクラス、波長、および反射率関係の情報を指定する構造体変数です。

lib

lib=1×8 struct array with fields:
    Name
    Type
    Class
    SubClass
    ParticleSize
    Genus
    Species
    SampleNo
    Owner
    WavelengthRange
    Origin
    CollectionDate
    Description
    Measurement
    FirstColumn
    SecondColumn
    WavelengthUnit
    DataUnit
    FirstXValue
    LastXValue
    NumberOfXValues
    AdditionalInformation
    Wavelength
    Reflectance
      ⋮

ECOSTRESS スペクトルライブラリから読み取ったスペクトルシグネチャをプロットします。

figure
hold on
for idx = 1:numel(lib)
    plot(lib(idx).Wavelength,lib(idx).Reflectance,'LineWidth',2);
end
axis tight
box on
xlabel('Wavelength (\mum)');
ylabel('Reflectance (%)');
classNames = {lib.Class};
legend(classNames,'Location','northeast')
title('Reference Spectra from ECOSTRESS Library');
hold off

テストデータの読み取り

関数 hypercube を使用して Jasper Ridge データセットからテストデータを読み取ります。この関数は、データキューブおよびテストデータから読み取ったメタデータ情報を格納する hypercube オブジェクトを返します。テストデータには、198 個のスペクトルバンドと、399.4 nm ～ 2457 nm の範囲の波長が含まれています。スペクトル分解能は最大 9.9 nm、各バンドイメージの空間分解能は 100 x 100 です。テストデータには、道路、土壌、水、および樹木を含む 4 つの潜在エンドメンバーが格納されています。

hcube = hypercube('jasperRidge2_R198.hdr');

エンドメンバースペクトルの抽出

テストデータに存在するスペクトルが特徴的なエンドメンバーの総数を計算するには、関数 countEndmembersHFC を使用します。この関数は、Harsanyi–Farrand–Chang (HFC) 法を使用してエンドメンバーの数を求めます。誤検知を防止するために偽警報の確率 (PFA) を低い値に設定します。

numEndmembers = countEndmembersHFC(hcube,'PFA',10^-27);

N-FINDR 法を使用してテストデータのエンドメンバーを抽出します。

endMembers = nfindr(hcube,numEndmembers);

hypercube オブジェクト hcube から波長の値を読み取ります。抽出したエンドメンバーシグネチャをプロットします。このテストデータは、4 つのエンドメンバー物質で構成されています。これらの物質のクラス名はスペクトルマッチングで同定できます。

figure
plot(hcube.Wavelength,endMembers,'LineWidth',2)
axis tight
xlabel('Wavelength (nm)')
ylabel('Data Values')
title('Endmembers Extracted using N-FINDR')
num = 1:numEndmembers;
legendName = strcat('Endmember',{' '},num2str(num'));
legend(legendName)

エンドメンバー物質の同定

エンドメンバー物質の名前を同定するには、関数 spectralMatch を使用します。この関数は、ライブラリファイルと分類対象のエンドメンバースペクトルの間のスペクトルの類似度を計算します。マッチングスコアの計算にはスペクトル情報発散 (SID) 法を選択します。通常、SID スコアの値が低いほど、テストスペクトルと基準スペクトルのマッチングが良好であることを意味します。次に、最適に一致する基準スペクトルのクラスに入るようにテストスペクトルを分類します。

たとえば、3 番目と 4 番目のエンドメンバー物質のクラスを同定するには、ライブラリシグネチャとそれぞれのエンドメンバースペクトルの間でスペクトルの類似度を求めます。最小 SID スコア値のインデックスは、スペクトルライブラリ内のクラス名を指定します。3 番目のエンドメンバースペクトルは海水と同定され、4 番目のエンドメンバースペクトルは樹木と同定されます。

wavelength = hcube.Wavelength;
detection = cell(1,1);
cnt = 1;
queryEndmember = [3 4];
for num = 1:numel(queryEndmember)
    spectra = endMembers(:,queryEndmember(num));
    scoreValues = spectralMatch(lib,spectra,wavelength,'Method','sid');
    [~, matchIdx] = min(scoreValues);
    detection{cnt} = lib(matchIdx).Class;
    disp(strcat('Endmember spectrum ',{' '},num2str(queryEndmember(num)),' is identified as ',{' '},detection{cnt}))
    cnt=cnt+1;
end

Endmember spectrum 3 is identified as Sea Water
Endmember spectrum 4 is identified as Tree

テストデータ内のエンドメンバー領域のセグメント化

同定結果を視覚的に検査するには、テストデータ内のエンドメンバー物質に固有のイメージ領域の位置を推定してセグメント化します。関数 sid を使用して、ピクセルスペクトルと抽出したエンドメンバースペクトルの間のピクセル単位のスペクトル類似度を計算します。次に、しきい値処理を実行し、テストデータの対象エンドメンバー領域をセグメント化して、セグメント化されたイメージを生成します。しきい値として 15 を選択し、最適に一致するピクセルを選択します。

可視化のため、関数 colorize を使用してテストデータの RGB バージョンを生成してから、セグメント化されたイメージをテストイメージの上に重ね合わせます。

threshold = 15;
rgbImg = colorize(hcube,'method','rgb','ContrastStretching',true);
overlayImg = rgbImg;
labelColor = {'Blue','Green'};
segmentedImg = cell(size(hcube.DataCube,1),size(hcube.DataCube,2),numel(queryEndmember));
for num = 1:numel(queryEndmember)
    scoreMap = sid(hcube,endMembers(:,queryEndmember(num)));
    segmentedImg{num} = scoreMap <= threshold;
    overlayImg = imoverlay(overlayImg,segmentedImg{num},labelColor{num});   
end

結果の表示

セグメント化されたイメージ、およびテストデータの海水エンドメンバー領域と樹木エンドメンバー領域を強調表示するオーバーレイイメージを表示して、同定結果を視覚的に検査します。

figure('Position',[0 0 900 400])
plotdim = [0.02 0.2 0.3 0.7;0.35 0.2 0.3 0.7];
for num = 1:numel(queryEndmember)
    subplot('Position',plotdim(num,:))
    imagesc(segmentedImg{num})
    title(strcat('Segmented Endmember region :',{' '},detection{num}));
    colormap([0 0 0;1 1 1])
    axis off
end

figure('Position',[0 0 900 400])
subplot('Position',[0 0.2 0.3 0.7])
imagesc(rgbImg)
title('RGB Transformation of Test Data');
axis off
subplot('Position',[0.35 0.2 0.3 0.7])
imagesc(overlayImg)
title('Overlay Segmented Regions')
hold on
dim = [0.66 0.6 0.3 0.3];
annotation('textbox',dim,'String','Sea Water','Color',[1 1 1],'BackgroundColor',[0 0 1],'FitBoxToText','on');
dim = [0.66 0.5 0.3 0.3];
annotation('textbox',dim,'String','Tree','BackgroundColor',[0 1 0],'FitBoxToText','on');
hold off
axis off

参考文献

[1] Kruse, F.A., A.B. Lefkoff, J.W. Boardman, K.B. Heidebrecht, A.T. Shapiro, P.J. Barloon, and A.F.H. Goetz. "The Spectral Image Processing System (SIPS)—Interactive Visualization and Analysis of Imaging Spectrometer Data." Remote Sensing of Environment 44, no. 2–3 (May 1993): 145–63. https://doi.org/10.1016/0034-4257(93)90013-N.