Шаг 1: импортируйте цветные инфракрасные каналы из многоспектрального файла изображения

Этот пример находит, что растительность в Тематическом изображении Картопостроителя LANDSAT, покрывающем часть Парижа, Франция, сделала доступным любезность Space Imaging, LLC. Семь спектральных каналов (полосы) хранятся в одном файле в формате LAN Эрдаса. Файл LAN, paris.lan, содержит с 7 каналами 512 512 изображение Landsat. 128-байтовый заголовок сопровождается пиксельными значениями, которые являются полосой, чередованной с методической точностью (BIL) в порядке увеличения числа полосы. Пиксельные значения хранятся как 8-битные целые числа без знака в прямом порядке байтов.

Первый шаг должен считать полосы 4, 3, и 2 из файла LAN с помощью функции MATLAB® multibandread.

Каналы 4, 3, и 2 покрывают почти инфракрасное излучение (NIR), видимый красный, и видимые зеленые части электромагнитного спектра. Когда они сопоставлены с красными, зелеными, и синими плоскостями, соответственно, изображения RGB, результатом является инфракрасное излучение стандартного цвета (CIR) составной объект. Итоговый входной параметр к multibandread задает, который полосы читать, и в который порядок, так, чтобы можно было создать составной объект на одном шаге.

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

Переменный CIR 512 512 3 массивами класса uint8. Это - изображение RGB, но с фиктивными цветами. Когда изображение отображено, красные пиксельные значения показывают канал NIR, зеленые значения показывают видимый красный канал, и синие значения показывают видимый зеленый канал.

В изображении CIR водные сооружения являются очень темными (река Сена), и зеленая растительность кажется красной (парки и тенистые деревья). Большая часть внешнего вида изображений - то, вследствие того, что здоровая, богатая хлорофиллом растительность имеет высокий коэффициент отражения в почти инфракрасное излучение. Поскольку канал NIR сопоставлен с красным каналом в составном изображении, любая область с высокой плотностью растительности кажется красной в отображении. Значимым примером является область яркого красного цвета на левом крае, большой парк (Bois de Boulogne) расположенный запад центрального Парижа в повороте реки Сены.

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 в отдельные переменные. Также полезно преобразовать от класса uint8 классифицировать single так, чтобы те же переменные могли использоваться в расчете NDVI ниже, а также в графике рассеивания.

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

Просматривая два канала вместе как полутоновые изображения, вы видите, как отличающийся они смотрят.

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

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

С одним простым вызовом plot команда в MATLAB, можно создать график рассеивания, отображающий одну точку на пиксель (как синий крест, в этом случае), с его x-координатой, определенной его значением в красном канале и его y-координатой значением его значение в канале NIR.

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® Array Arithmetic

Заметьте из графика рассеивания, что взятие отношения уровня NIR к красному уровню было бы одним способом определить местоположение пикселей, содержащих плотную растительность. Однако результат был бы шумным для темных пикселей с маленькими значениями в обоих каналах. Также заметьте, что различие между NIR и красными каналами должно быть больше для большей плотности хлорофилла. Нормированный индекс растительности различия (NDVI) мотивирован этим вторым наблюдением. Это берет (NIR - красный) различие и нормирует его, чтобы помочь балансировать эффекты неровного освещения, такие как тени облаков или выступов. Другими словами, на попиксельной основе вычитают значение красного канала от значения NIR, образовывают канал и делятся на их сумму.

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

Заметьте, как арифметические операторы массивов в MATLAB позволяют вычислить целое изображение NDVI в одной простой команде. Вспомните это переменные R и NIR имейте класс single. Этот выбор использует меньше устройства хранения данных, чем класс double но различающийся целочисленный класс также позволяет получившемуся отношению принимать сглаженную градацию значений.

Переменный ndvi 2D массив класса single с теоретической максимальной областью значений [-1 1]. Можно задать эти теоретические пределы при отображении ndvi как полутоновое изображение.

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

Река Сена кажется очень темной в изображении NDVI. Большой легкой областью около левого края изображения является парк (Bois de Boulogne), отмеченный ранее.

Шаг 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-красном графике рассеивания, перерисовав график рассеивания с пикселями выше порога в контрастирующем (зеленом) цвете и затем вновь отобразив порог изображение 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')