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

Этот пример находит растительность на изображении тематического картографа LANDSAT, охватывающем часть Парижа, Франция, которое предоставлено компанией Space Imaging, LLC. Семь спектральных каналов (диапазонов) хранятся в одном файле в формате LAN Erdas. файл LAN, paris.lan, содержит 7-канальное изображение Ландсата 512 на 512. За 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-канал отображается на красный канал в составном изображении, любая область с высокой плотностью растительности отображается на дисплее красным цветом. Заметный пример - область ярко-красного цвета на левом краю, большой парк (Буа-де-Булонь), расположенный к западу от центрального Парижа в пределах изгиба реки Сена.

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

Анализируя различия между НДК и красными каналами, можно количественно определить этот контраст в спектральном содержании между растительными участками и другими поверхностями, такими как дорожное покрытие, голый грунт, здания или вода.

Шаг 2: Построение графика NIR-красного спектрального рассеяния

График рассеяния является естественным местом для начала при сравнении 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 ®

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

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

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

Переменная ndvi является 2-D массивом класса single с теоретическим максимальным диапазоном [-1 1]. Эти теоретические пределы можно задать при отображении ndvi в градациях серого.

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

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

Шаг 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')