手順 1: 複数スペクトル イメージ ファイルからのカラー赤外チャネルのインポート

この例では、Space Imaging, LLC の提供により、フランス、パリの一部を撮影した LANDSAT Thematic Mapper イメージ内の植生を検出します。7 つのスペクトル チャネル (帯域) は Erdas LAN 形式の 1 つのファイルに格納されます。LAN ファイル paris.lan には、7 チャネルの 512×512 の Landsat イメージが含まれています。128 バイトのヘッダーに続くピクセル値は、帯域番号の昇順によるライン挟み込みバンド並び (BIL) 形式です。ピクセル値は、リトルエンディアン バイトの順番で符号なし 8 ビット整数として保存されます。

最初の手順では、MATLAB® 関数 multibandread を使用して LAN ファイルから帯域 4、3 および 2 を読み取ります。

チャネル 4、3、および 2 は、電磁スペクトルの近赤外線 (NIR)、可視赤色光、および可視緑色光の部分を網羅します。それぞれが RGB イメージの赤、緑、および青の平面にマッピングされる場合、その結果は標準のカラー赤外 (CIR) 合成イメージになります。multibandread への最後の入力引数では、読み取る帯域と、その順序を指定します。そのため、1 ステップで合成イメージを作成できます。

CIR = multibandread('paris.lan',[512, 512, 7],'uint8=>uint8',...
                    128,'bil','ieee-le',{'Band','Direct',[4 3 2]});

変数 CIR は uint8 クラスの 512×512×3 の配列です。これは RGB イメージですが、フォールス カラーをもっています。イメージを表示する場合、赤のピクセル値は NIR チャネル、緑の値は可視赤色チャネル、青の値は可視緑色チャネルをそれぞれ示します。

CIR イメージでは、水は非常に暗い特徴をもち (セーヌ川)、緑の植生は赤く見えます (公園および街路樹)。イメージの外観の多くは、正常でクロロフィルが豊富な植生のため、近赤外線領域での反射率が高くなっています。NIR チャネルは合成イメージの赤のチャネルにマッピングされるので、植生密度が高い領域はディスプレイに赤く表示されます。注目すべき例は、セーヌ川の湾曲部内のパリ中心部の西側に位置する大きな公園 (ブーローニュの森) である左端の明るい赤の領域です。

imshow(CIR)
title('CIR Composite')
text(size(CIR,2),size(CIR,1) + 15,...
  'Image courtesy of Space Imaging, LLC',...
  'FontSize',7,'HorizontalAlignment','right')

NIR チャネルと赤のチャネル間の差異を解析することにより、このコントラストを、植生地域とその他の外見 (路面、露出土壌、建物、水など) の間のスペクトル量で定量化できます。

手順 2: NIR-Red スペクトル散布図の構築

NIR チャネル (赤のピクセル値で表示) と可視赤色チャネル (緑のピクセル値で表示) を比較する場合、散布図から始めます。これは、これらのチャネルを元の CIR 合成イメージから抽出し、独立した変数にするのに便利です。これはまた、散布図だけでなく、以下の NDVI 計算でも同じ変数を使用できるので、uint8 クラスを single クラスに変換するときにも役立ちます。

NIR = im2single(CIR(:,:,1));
R = im2single(CIR(:,:,2));

差異を見比べるため、2 つのチャネルを一緒にグレースケール イメージとして表示します。

imshow(R)
title('Visible Red Band')

imshow(NIR)
title('Near Infrared Band')

MATLAB での plot コマンドを呼び出すだけで、赤のチャネルの値により決定されるその x 座標と、NIR チャネルの値による y 座標と共に、1 ピクセルに 1 つの点 (この場合は青い十字) を表示する散布図を作成できます。

plot(R,NIR,'+b')
ax = gca;
ax.XLim  = [0 1];
ax.XTick = 0:0.2:1;
ax.YLim  =  [0 1];
ax.YTick = 0:0.2:1;
axis square
xlabel('red level')
ylabel('NIR level')
title('NIR vs. Red Scatter Plot')

パリの景色の散布図の外観は、青葉がある温暖な市街地を特徴としています。NIR と赤の値がほぼ等しい対角線の近くに一連のピクセルがあります。この "グレーの境界" には、路面や多数の屋根などの特徴が含まれています。上側と左側は、多くの場合に NIR 値が赤の値をかなり上回る別の一連のピクセルです。このゾーンは、基本的にすべての緑の植生を取り囲みます。

手順 3: MATLAB® 配列演算による植生インデックスの計算

散布図から、NIR レベル対赤のレベルの比率を読み取ることが、密集した植生を含むピクセルを位置付ける方法の 1 つであることに気付きます。ただし、その結果は、両方のチャネルの小さな値を持つ暗いピクセルについてはノイズを含みます。また、NIR チャネルと赤のチャネルの間の差異は、クロロフィル密度が高いほど大きくなることに注意してください。正規化植生指標 (NDVI) は、この 2 番目の観察によって変動します。これにより (NIR - Red) 差異を読み取って正規化し、雲や丘の影のような不均一な照度の影響を相殺するのに役立てます。つまり、ピクセルごとに、赤のチャネルの値を NIR チャネルの値から減算し、その合計で除算します。

ndvi = (NIR - R) ./ (NIR + R);

MATLAB の配列算術演算子により、1 つの簡単なコマンドだけで、NDVI イメージ全体を計算することが可能になることに注目してください。変数 R および NIR には single クラスがあることを思い出してください。この選択により使用されるメモリ量は double クラスより少ないですが、整数クラスとは異なり、結果の比率で値の滑らかなグラデーションを推測することも可能です。

変数 ndvi は、理論上の最大範囲が [-1 1] である single クラスの 2 次元配列です。ndvi をグレースケール イメージとして表示すると、これらの理論上の範囲を指定できます。

figure
imshow(ndvi,'DisplayRange',[-1 1])
title('Normalized Difference Vegetation Index')

セーヌ川は、NDVI イメージでは非常に暗く表示されます。イメージの左端付近の大きな明るい領域は、前述の公園 (ブーローニュの森) です。

手順 4: 植生の配置 -- NDVI イメージのしきい値処理

多くの植生を含んでいる可能性が最も高いピクセルを特定するには、単純なしきい値を NDVI イメージに適用します。

threshold = 0.4;
q = (ndvi > threshold);

したがって、選択したピクセルの割合は以下のようになります。

100 * numel(NIR(q(:))) / numel(NIR)

ans = 5.2204

つまり、約 5% になります。

公園や植生の他の小さな領域は、論理 (バイナリ) イメージ q を表示したときに、既定で白く表示されます。

imshow(q)
title('NDVI with Threshold Applied')

手順 5: スペクトル量と空間量のリンク

スペクトル量と空間情報をリンクするには、しきい値より上のピクセルを NIR-Red の散布図に配置します。これにより、しきい値より上のピクセルをもつ散布図は高コントラストの色 (緑) で再描画され、しきい値 NDVI イメージは同じ青緑色のスキームで再表示されます。予想どおり、しきい値より上の NDVI 値をもつピクセルは、その他のピクセルの左上に現れ、CIR 合成イメージ表示のより赤いピクセルと一致します。

散布図を作成し、しきい値処理された NDVI を表示します。

figure
subplot(1,2,1)
plot(R,NIR,'+b')
hold on
plot(R(q(:)),NIR(q(:)),'g+')
axis square
xlabel('red level')
ylabel('NIR level')
title('NIR vs. Red Scatter Plot')

subplot(1,2,2)
imshow(q)
colormap([0 0 1; 0 1 0]);
title('NDVI with Threshold Applied')