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

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

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

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

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

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

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

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

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

Река Сена выглядит очень темной на изображении НДВИ. Большой световой площадью около левого края изображения является отмеченный ранее парк (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')