Свой маленький GIS: WPF-приложение для мультиспектральных и гиперспектральных снимков

Привет, Хабр. Меня зовут Алексей, я C#-разработчик. В этой статье хочу рассказать о своём дипломном проекте очень запавшем мне в душу, который я делал на тему обработки изображений, GIS и дистанционного зондирования Земли. Даже спустя годы мне интересна данная тема и она по-прежнему остаётся очень перспективной в различных отраслях.

Идея была в том, чтобы собрать небольшое настольное приложение, которое умеет работать с реальными спутниковыми данными: Landsat 8, Sentinel-2 и AVIRIS. То есть открывать не готовую RGB-картинку, а набор спектральных каналов, собирать из них естественные и псевдоцветные изображения, считать растровые индексы, выделять эталоны прямо на снимке, классифицировать пиксели, сегментировать изображение и пробовать более исследовательские вещи вроде EMD-разложения.

В итоге получилась учебно-исследовательская программа, но с полным рабочим циклом: от чтения спутникового архива до сохранения информативного результата обработки. Ниже расскажу, зачем вообще нужны такие снимки, какие особенности есть у разных спутниковых данных, что делает приложение и какие алгоритмы оказались самыми интересными.

Зачем это всё

Когда мы смотрим на обычную фотографию, у нас есть привычные три канала: красный, зеленый и синий. Этого достаточно, чтобы получить картинку, похожую на то, что увидел бы человек. Но спутниковая съемка почти никогда не ограничивается человеческим зрением.

Мультиспектральные и гиперспектральные снимки хранят информацию в десятках или даже сотнях каналов. Часть этих каналов находится в видимом диапазоне, часть — в ближнем инфракрасном, коротковолновом инфракрасном и других областях спектра. Для глаза это невидимо, но для анализа очень полезно. Например, растительность хорошо проявляется в ближнем инфракрасном диапазоне. Вода, песок, городская застройка, облака, сельскохозяйственные культуры и загрязнения тоже имеют разные спектральные признаки. Поэтому спутниковый снимок можно использовать не только как красивую картинку сверху, но и как набор измерений, из которых можно извлекать практическую информацию.

Такие данные применяются в сельском хозяйстве, мониторинге лесов и водоемов, оценке последствий пожаров и наводнений, поиске нефтяных загрязнений, геологии, строительстве, метеорологии и многих других задачах. Главный плюс дистанционного зондирования в том, что мы можем наблюдать большие территории регулярно, не вмешиваясь в происходящие процессы и не отправляя человека туда, где это дорого или опасно.

Но есть и обратная сторона. Сырые данные неудобны. Они большие, приходят в разных форматах, требуют понимания каналов и метаданных, а для анализа часто нужно писать свои инструменты. Поэтому в дипломе я хотел сделать программу, которая закрывает хотя бы базовый исследовательский сценарий: загрузил снимок, выбрал каналы, посмотрел композицию, применил алгоритм и получил результат.

Что получилось

Приложение написано на C#, для работы с изображениями используется OpenCvSharp, для ускорения обработки — TPL, для интерфейса — обычный десктопный WPF с меню, панелью инструментов, прогрессом операций и просмотром изображения с масштабированием.

Программа умеет:

собирать RGB-изображение из выбранных каналов;
работать с данными спутников Landsat 8, Sentinel-2 и AVIRIS (данные в интернете в открытом доступе);
строить естественные и псевдоцветные композиции;
обрезать область снимка;
считать растровые выражения вида (Ch1 - Ch2) / (Ch1 + Ch2);
строить гистограммы каналов;
добавлять пользовательские эталоны по выделенной области;
классифицировать изображение по эталонам;
считать площадь найденных классов в пикселях, метрах и процентах;
сегментировать изображение по спектральному углу;
выполнять экспериментальное EMD-разложение для AVIRIS;
сохранять результаты в TIFF;

С точки зрения пользователя сценарий довольно прямой. Сначала выбирается тип спутника и папка со снимком. Затем указываются номера каналов, из которых нужно собрать изображение. После этого снимок появляется в интерфейсе, и к нему можно применять операции.

Пример снимка Sentinel после применения классификации

Почему данные разных спутников так отличаются

Одна из первых вещей, которая быстро становится понятной при работе со спутниковыми снимками: «открыть файл» недостаточно. У разных источников разные форматы поставки, структура папок, метаданные, разрешение каналов и даже способ хранения значений.

Landsat 8

Landsat 8 поставляет каналы отдельными TIFF-файлами и метаданные в XML. Это довольно удобный вариант для старта: каждый канал можно открыть отдельно, а затем собрать нужную комбинацию.

В программе для Landsat используются номера каналов. Например, классическая композиция в естественных цветах собирается из каналов 4, 3 и 2: красный, зеленый и синий. Для анализа растительности уже интереснее брать ближний инфракрасный и коротковолновой инфракрасный диапазоны (какой спектр под каким номером легко найти в интернете для любого спутника). В коде для Landsat есть отдельный класс, который умеет находить нужный файл канала в директории и читать пространственное разрешение из метаданных.

Sentinel-2

Sentinel-2 обычно приходит в структуре SAFE. Там есть директория GRANULE, внутри нее данные изображения, а сами каналы лежат в JP2. Дополнительная особенность Sentinel-2 в том, что каналы могут иметь разное пространственное разрешение: 10, 20 или 60 метров. Для простого просмотра это уже создает проблему. Если пользователь выбрал три канала, а один из них имеет другое разрешение, их нельзя просто сложить в RGB-массив пиксель к пикселю. Поэтому в программе каналы Sentinel-2 приводятся к одному размеру.

Работает это так: за целевой размер берется первый выбранный канал, а остальные каналы интерполируются до его количества строк и столбцов. Если первый канал имеет размер W x H, то для каждого пикселя результата с координатами (x, y) библиотека берет значение из исходного канала в дробной координате:

$x_{src} = x \cdot \frac{W_{src}}{W_{dst}}$ $y_{src} = y \cdot \frac{H_{src}}{H_{dst}}$

При линейной интерполяции значение нового пикселя получается как взвешенная сумма ближайших исходных пикселей. В коде это выглядит так:

int targetRows = bands[0].Rows;int targetCols = bands[0].Cols;for (int i = 1; i < bands.Length; i++){    if (bands[i].Rows == targetRows && bands[i].Cols == targetCols)        continue;    Mat resized = new Mat();    Cv2.Resize(bands[i], resized, new Size(targetCols, targetRows),        0, 0, InterpolationFlags.Linear);    bands[i].Dispose();    bands[i] = resized;}

Корректно ли так делать? Для визуализации и быстрых экспериментов — я считаю это нормальный компромисс, иначе мы вообще не сможем собрать три канала в одну картинку. Но если задача научная и по результату потом принимаются решения, лучше относиться к этому как к упрощению, при понижении разрешения часть деталей сглаживается. Поэтому для строгой обработки лучше явно выбирать единую сетку, учитывать геопривязку, метод ресемплинга и не смешивать каналы разных разрешений без понимания последствий.

AVIRIS

AVIRIS — это уже настоящий гиперспектральный снимок. В отличие от мультиспектральных Landsat и Sentinel-2, здесь каналов намного больше. В используемых данных AVIRIS — 224 канала, а спектральный диапазон покрывает примерно 366-2500 нм. Такие данные можно представить как куб: две координаты отвечают за положение пикселя на изображении, а третья — за спектр в этой точке. Для каждого пикселя у нас получается не просто три числа RGB, а целая сигнатура.

Дополнительная сложность в том, что гиперспектральные данные могут храниться по-разному. Часто встречаются форматы BIL, BIP и BSQ:

— BIL хранит строки каналов последовательно;

— BIP хранит сначала все каналы одного пикселя, затем следующего;

— BSQ хранит сначала весь первый канал, потом весь второй и так далее.

BIP – сначала записываются значения всех каналов 1-го пикселя, затем 2-го и т. д.

В проекте для AVIRIS реализовано чтение данных из бинарного файла и заголовка. Пользователь может указать номера каналов, а программа соберет из них изображение. Можно выбрать три канала и получить RGB-композицию, а можно выбрать один канал и посмотреть его в оттенках серого.

Сборка изображения из каналов

Самая базовая операция в приложении — собрать видимую картинку из выбранных каналов. На первый взгляд звучит просто: взяли три массива и положили один в R, второй в G, третий в B, все осмысленные комбинации легко загуглить, но на практике появляются детали. Формально, исходный снимок можно представить как набор каналов: I(x, y, k), где x и y — координаты пикселя, k — номер спектрального канала.

Чтобы получить обычное RGB-изображение, мы выбираем три канала r, g, b и строим новый пиксель: RGB(x, y) = (I(x, y, r), I(x, y, g), I(x, y, b)). Например, для естественных цветов у Landsat-8 удобно брать каналы 4, 3 и 2. Для псевдоцветной композиции можно поставить ближний инфракрасный канал в красную компоненту. Тогда растительность становится заметнее, потому что хорошо отражает в NIR-диапазоне.

Снимок Landsat, собранный из 5, 6 и 4 каналов соответственно

На картинке комбинация ближнего, среднего ИК-каналов и красного видимого канала позволяет четко различить границу между водой и сушей и подчеркнуть скрытые детали плохо видимые при использовании только каналов видимого диапазона. С большой точностью будут детектироваться водные объекты внутри суши. Эта комбинация отображает растительность в различных оттенках и тонах коричневого, зеленого и оранжевого. Эта комбинация дает возможность анализа влажности и полезны при изучении почв и растительного покрова. В целом, чем выше влажность почв, тем темнее она будет выглядеть, что обусловлено поглощением водой излучения ИК диапазона.

Но уже на этапе сборки без сложностей не обойтись:

Во-первых, данные могут быть 16-битными, а для отображения в интерфейсе часто нужен 8-битный диапазон. Поэтому значения масштабируются и ограничиваются диапазоном 0-255. В простом виде это можно записать так:

$v_8 = \operatorname{clamp}(v_{16} \cdot coef,\ 0,\ 255)$

Где coef — коэффициент яркости, подобранный для конкретного источника данных.

Во-вторых, разные спутники требуют разного чтения. Landsat открывается через TIFF-декодер, Sentinel-2 читается как JP2, AVIRIS разбирается из бинарного файла, поэтому сборка изображения разведена по типам спутников.

Для Landsat данные приходят как 16-битный TIFF. В моей реализации я беру старший байт каждого значения, умножаю его на коэффициент яркости и раскладываю выбранные каналы по BGR-порядку OpenCV:

vecs[i] = new Vec3b(  ClipToByte(bytes[2][i * 2 + 1]*satellite.brightCoef),  ClipToByte(bytes[1][i * 2 + 1]*satellite.brightCoef),  ClipToByte(bytes[0][i * 2 + 1]*satellite.brightCoef));

Здесь может смутить порядок 2, 1, 0, но это обычная особенность OpenCV: цветное изображение хранится как BGR, а не RGB.

AVIRIS читается иначе. Там данные лежат в бинарном файле, и для каждого пикселя нужно достать значения нужных каналов. В моем случае значение канала занимает два байта:

int c1 = blue + j / 3 * bands * 2;int c2 = green + j / 3 * bands * 2;int c3 = red + j / 3 * bands * 2;res[j] = ClipToByte(BitConverter.ToInt16(new byte[] { bytes[c1 + 1], bytes[c1] }, 0) / 256.0 * sat.brightCoef);

Это уже ближе к работе не с картинкой, а с большим бинарным массивом, где важно правильно посчитать смещения.

В-третьих, обычная RGB-композиция не всегда самая полезная. Для Sentinel-2 можно взять каналы 4, 3 и 2 и получить изображение, близкое к естественному цвету. Но можно собрать псевдоцветную композицию, где растительность, вода или влажные участки станут заметнее. В этом и состоит одна из приятных особенностей мультиспектральных данных: мы не обязаны смотреть только «как глазом».

Калькулятор растров

Следующий полезный инструмент — вычисление выражений по каналам. Самый известный пример — вегетационный индекс NDVI: (NIR — Red) / (NIR + Red). Он работает за счёт того, что здоровые растения активно поглощают красный свет (для фотосинтеза) и сильно отражают инфракрасный свет. В программе это сделано обобщенно: пользователь может указать выражение через Ch1, Ch2, Ch3, а приложение заменяет их на значения и считает формулу для каждого пикселя. Внутри используется разбор выражения в обратную польскую запись. Например: (Ch1 — Ch2) / (Ch1 + Ch2).

Дальше результат нужно как-то отобразить. В коде положительные и отрицательные значения раскладываются по цветам: один канал усиливается при росте значения индекса, другой — при уменьшении. Получается не просто числовой растр, а сразу наглядная картинка, где области с разными значениями хорошо различимы.

Практически это удобно для быстрых экспериментов. Можно не зашивать отдельной кнопкой каждый индекс, а дать пользователю формулу. Сначала он считает NDVI, затем пробует EVI индекс, потом CVI.

Пример NDVI: чем зеленее, тем гуще растительность

Классификация по эталонам

В приложении реализована классификация с обучением, то есть пользователь сам выбирает эталонные области на снимке. Например, выделяет прямоугольник с водой, затем участок леса, затем песок или городскую застройку. Для каждого эталона можно выбрать цвет, которым будут отмечаться найденные пиксели.

После выделения области программа считает ее спектральное представление. В текущей реализации эталон хранится как небольшая матрица значений, выбранных по каналам. Затем каждый пиксель изображения сравнивается с каждым эталоном. Если он проходит порог по выбранному критерию, пиксель перекрашивается в цвет класса.

Дополнительно считается статистика:

— сколько пикселей попало в каждый класс;

— какая это площадь с учетом разрешения спутника;

— какой процент изображения занимает класс.

Пример классификации на снимке Sentinel (с порогом не игрался, получилось как получилось)

Сам обход изображения устроен довольно просто, но важная часть — параллельность. Спутниковые снимки большие: например, Sentinel-2 может быть 10980х10980 пикселей. Последовательный обход таких массивов это безумие, поэтому классификация должна быть хоть как-то распараллелена, а прогресс отправляется в интерфейс через BackgroundWorker.

Самое интересное здесь не в цикле, а в критериях сходства.

Спектральный угол

Первый критерий — спектральный угол. Мы берем пиксель как вектор и эталон как другой вектор. Если у нас три выбранных канала, то это векторы в трехмерном пространстве. Дальше считаем угол между ними:

$\cos{\theta} = \frac{\mathbf{a} \cdot \mathbf{b}}{\|\mathbf{a}\| \cdot \|\mathbf{b}\|}$ $\theta = \arccos{\left(\frac{\mathbf{a} \cdot \mathbf{b}}{\|\mathbf{a}\| \cdot \|\mathbf{b}\|}\right)}$

Если расписать для трех каналов, получится:

$\mathbf{a} \cdot \mathbf{b} = a_R b_R + a_G b_G + a_B b_B$ $\|\mathbf{a}\| = \sqrt{a_R^2 + a_G^2 + a_B^2}$ $\|\mathbf{b}\| = \sqrt{b_R^2 + b_G^2 + b_B^2}$

Если угол маленький, значит форма спектральной сигнатуры похожа. Это важное отличие от простого сравнения яркости. Два пикселя могут быть один темнее другого, но относительное соотношение каналов у них останется близким. Для спутниковых данных это часто полезно, потому что освещенность, тени и атмосферные эффекты могут менять абсолютные значения.

В коде это выглядит как обычное скалярное произведение, нормы векторов и Math.Acos. Перед вычислением косинус зажимается в диапазон от -1 до 1, чтобы избежать проблем с погрешностями.

double scalar = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];double k1 = Math.Sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);double k2 = Math.Sqrt(b[0] * b[0] + b[1] * b[1] + b[2] * b[2]);double cos = scalar / (k1 * k2);cos = Math.Max(-1, Math.Min(1, cos));double angle = Math.Acos(cos) * 180 / Math.PI;

На практике спектральный угол хорошо подходит для выделения объектов с похожей спектральной формой. В экспериментах он давал более четкие границы, чем часть простых метрик, хотя работал медленнее евклидова расстояния.

Евклидово расстояние

Второй критерий — евклидово расстояние. Здесь все интуитивно: считаем расстояние между точкой-пикселем и точкой-эталоном в пространстве признаков.

$d = \sqrt{ (r_1 - r_2)^2 + (g_1 - g_2)^2 + (b_1 - b_2)^2 }$

В реализации разница по каналам дополнительно нормируется относительно значений эталона, чтобы один канал не перетягивал результат слишком сильно.

Плюс евклидова расстояния — скорость и простота. Его легко объяснить, легко отладить, легко подобрать порог. Минус — чувствительность к яркости. Если объект того же типа оказался в тени или снят с другой интенсивностью, расстояние может вырасти, хотя спектрально он все еще похож на эталон.

В моих замерах на снимке Landsat 8 размером 8081х8161 пикселей евклидово расстояние оказалось самым быстрым из трех критериев:

Критерий	1 эталон	2 эталона
Евклидово расстояние	6.9с	12с
Спектральный угол	9.45с	16.55с
Барицентрические координаты	51.5с	82.65с

Эти цифры не стоит воспринимать как универсальный бенчмарк: они зависят от компьютера, размера снимка и реализации. Но соотношение хорошо показывает характер методов.

Барицентрические координаты

Самым необычным алгоритмом в проекте для меня стала классификация через барицентрические координаты. Я не встречал такой подход в прикладных материалах по обработке спутниковых снимков, поэтому было особенно интересно попробовать. Обычно барицентрические координаты объясняют через треугольник: есть вершины A, B и C, а положение точки X можно выразить через веса относительно этих вершин. Если точка внутри треугольника, веса ведут себя «спокойно», а если выходит далеко наружу, координаты начинают показывать отклонение.

В приложении идея переносится в пространство каналов. Эталон задает три точки, которые можно воспринимать как вершины в RGB-пространстве выбранной композиции. Текущий пиксель — это точка, которую мы хотим проверить. Если его барицентрические координаты относительно эталона не выходят за заданный порог, считаем пиксель похожим на эталон.

Знаки каждой координаты при разном расположении точки, на гранях АВС соответствующие координаты равны 0

В матричном виде это можно записать так:

$\begin{bmatrix} \lambda_1 \\ \lambda_2 \\ \lambda_3 \end{bmatrix} = E^{-1}X$ $E = \begin{bmatrix} e_{1R} & e_{2R} & e_{3R} \\ e_{1G} & e_{2G} & e_{3G} \\ e_{1B} & e_{2B} & e_{3B} \end{bmatrix}, \quad X = \begin{bmatrix} x_R \\ x_G \\ x_B \end{bmatrix}$

Где E — матрица эталона, X — вектор текущего пикселя, а lambda — барицентрические координаты. В классическом геометрическом случае для точки внутри треугольника выполняется:

$\lambda_1 + \lambda_2 + \lambda_3 = 1$

А сами коэффициенты показывают вклад каждой вершины. В моей задаче это скорее не строгая геометрия треугольника на плоскости, а практический критерий: насколько текущий пиксель выражается через выбранные эталонные компоненты без сильного выброса по коэффициентам.

В коде для каждого эталона заранее строится обратная матрица:

Matrix m = new Matrix(3, 3);for (int i = 0; i < 3; i++)    for (int j = 0; j < 3; j++)        m[i, j] = model.data[j, i];Matrix im = m.CreateInvertibleMatrix();

А затем для каждого пикселя выполняется умножение:

Matrix res = inverted_matrix[k] * r;if (Math.Abs(res[0, 0]) < threshold &&    Math.Abs(res[1, 0]) < threshold &&    Math.Abs(res[2, 0]) < threshold){    // пиксель относится к эталону}

Метод получился самым тяжелым по времени, потому что матричные операции выполняются для большого количества пикселей. Зато он интересен тем, что оценивает пиксель не просто как точку рядом с центром эталона, а относительно формы, заданной самим эталоном.

На снимке Sentinel-2 с выделенными эталонами песка, воды и леса метод дал наглядный результат: вода выделялась с четкими контурами, песчаные участки тоже определялись достаточно уверенно. Порог около 8 в моих экспериментах выглядел рабочим стартовым значением, хотя, конечно, универсального числа здесь нет.

Результат классификации на барицентрических координатах, 3 эталона

Сегментация как граф пикселей

Классификация отвечает на вопрос: «похож ли пиксель на один из известных эталонов?». Сегментация решает другую задачу: разбить изображение на области, внутри которых пиксели похожи друг на друга. В проекте реализована графо-ориентированная сегментация по идее Felzenszwalb и Huttenlocher. Изображение рассматривается как неориентированный взвешенный граф:

вершины — пиксели;
ребра — связи между соседними пикселями;
вес ребра — отличие между пикселями;

В качестве отличия я использовал спектральный угол. Это логично для мульти- и гиперспектральной задачи: нам важно сравнивать не только яркость, но и направление спектрального вектора.

Алгоритм в общих чертах такой:

Сначала изображение сглаживается гауссовым фильтром, чтобы уменьшить шум.
Для соседних пикселей строятся ребра с весом, равным спектральному углу.
Ребра сортируются по весу.
Изначально каждый пиксель считается отдельным компонентом.
Компоненты объединяются, если вес ребра между ними достаточно мал относительно внутренних перепадов.
Маленькие компоненты дополнительно склеиваются на постобработке.

Для хранения компонент используется структура непересекающихся множеств, то есть union-find. Это хорошо ложится на задачу: нужно быстро находить, в какой компоненте находится пиксель, и объединять две компоненты. Математически здесь есть две важные величины. Первая — внутреннее различие компоненты, то есть максимальный вес ребра внутри уже собранного сегмента, вторая — порог объединения, который зависит от размера компоненты. Две компоненты объединяются, если ребро между ними легче этого порога:

$w(C_1, C_2) \le \operatorname{MInt}(C_1, C_2)$

Интуитивно это означает: маленькие компоненты можно объединять смелее, а крупные должны иметь более убедительное сходство на границе.

Пользователь может управлять тремя параметрами:

sigma — степень гауссова размытия;
threshold — насколько охотно сегменты объединяются;
min_size — минимальный размер сегмента;

Если увеличить детализацию, сегментов станет больше, границы будут мельче. Если поднять пороги и минимальный размер, результат станет более обобщенным. В дипломном эксперименте на AVIRIS-снимке размером 753х1924 пикселя сегментация занимала около 13 секунд.

Сегментация при sigma = 0.5, threshold = 100, min_size = 50

Сравнение результатов при разных параметрах

EMD: разложение спектра пикселя на моды

Еще одна экспериментальная часть проекта — EMD (Empirical Mode Decomposition). Это метод разложения нелинейных и нестационарных сигналов на эмпирические моды. В контексте гиперспектрального изображения сигналом можно считать спектральную сигнатуру одного пикселя: значения интенсивности по 224 каналам AVIRIS. То есть для каждого пикселя у нас есть не просто цвет, а последовательность чисел по спектру. EMD позволяет разложить эту последовательность на составляющие: первые моды отвечают за более высокочастотные изменения, поздние — за низкочастотные.

Если проводить аналогию со звуком, то спектральная сигнатура пикселя похожа на короткую мелодию, сыгранную не во времени, а по длинам волн. В ней может быть общий плавный тренд, резкие локальные всплески, шум и небольшие колебания. EMD пытается разложить эту «мелодию» на несколько слоев: отдельно быстрые колебания, отдельно более медленные, отдельно остаточный тренд.

В классическом виде разложение можно записать так:

$f(t) = \sum_{i=1}^{N} \operatorname{IMF}_i(t) + r(t)$

Где f(t) — исходный сигнал, IMF_i(t) — эмпирические моды, а r(t) — остаток. Для гиперспектрального снимка вместо времени t можно мысленно подставить номер спектрального канала.

То есть каждый пиксель превращается в одномерный сигнал по спектральной оси.

Упрощенно алгоритм в программе выглядит так:

Берем спектральную сигнатуру пикселя.
Считаем сглаженный остаток через бегущее окно.
Вычитаем остаток из текущего сигнала и получаем моду.
Повторяем процесс до выбранного номера моды.
Отображаем выбранную моду как изображение.

В реализации EMD работает с AVIRIS, потому что там достаточно каналов, чтобы сама идея имела смысл. В моей реализации сглаженный остаток считается через бегущее окно шириной w = 3:

$R_i = \frac{1}{w} \sum_{j=i-\lfloor w/2 \rfloor}^{i+\lfloor w/2 \rfloor} f_j$

Если окно выходит за границы спектра, берется крайнее доступное значение. После этого мода считается как разность:

В коде это выглядит примерно так:

for (int b = 0; b < bands; b++){    short Ri = 0;    for (int k1 = -w / 2; k1 <= w / 2; k1++)    {        int x = b + k1;        if (x < 0) x = 0;        if (x >= bands) x = bands - 1;        Ri = (short)(F[x][i] + Ri);    }    Ri = (short)(Ri / w);    mode[b][i] = (short)(F[b][i] - Ri);    F[b][i] = Ri;}

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

Слева исходное изображение, далее по порядку 1-3 моды 31 канала

Пиксель, выбранный для эксперимента с EMD

а - исходная гистограмма пикселя морской воды, б-д это 1-4 моды, е - одна из последних, 15-я мода — а — исходная гистограмма пикселя морской воды, б-д это 1-4 моды, е — одна из последних, 15-я мода

Для первых трех мод в дипломных экспериментах получилось такое время:

Номер моды	Время
1	27.2с
2	49.7с
3	70.75с

Метод тяжелый, но перспективный. Его можно использовать для поиска шумовых составляющих, выделения тонких спектральных признаков, последующего сжатия или подготовки признаков для классификации.

Я бы не стал называть текущую реализацию EMD промышленной. Скорее это исследовательский модуль внутри программы: он показывает, что гиперспектральный снимок можно анализировать не только как набор картинок по каналам, но и как множество сигналов.

Немного про производительность

Снимки в дистанционном зондировании быстро учат уважать размер данных. Даже простая операция «пройтись по всем пикселям» становится заметной, когда изображение имеет десятки или сотни миллионов пикселей. Поэтому в проекте многие операции сделаны через линейные массивы и параллельные циклы. Вместо обращения к пикселям через высокоуровневые методы в каждом шаге изображение превращается в массив Vec3b[] или byte[], а затем обрабатывается напрямую. Это не делает программу магически быстрой, но убирает лишние накладные расходы.

Например, классификация работает не через GetPixel/SetPixel, а через массив пикселей:

Vec3b[] array = imageInfo.GetBytes();int progressStep = Math.Max(1, array.Length / 100);Parallel.For(0, array.Length, (i) =>{    if (!IsZero(array[i]))    {        for (int k = 0; k < models.Count; k++)        {            double angle = GetAngle(modelVector, array[i]);            if (Math.Abs(angle) < threshold)            {                Interlocked.Increment(ref coverCount[k]);                array[i] = new Vec3b(                    models[k].userColor.B,                    models[k].userColor.G,                    models[k].userColor.R);                break;            }        }    }    int currentProgress = Interlocked.Increment(ref progress);    if (currentProgress % progressStep == 0)        backgroundWorker.ReportProgress(            Math.Min(99, currentProgress / progressStep));});

Parallel.For распределяет пиксели между потоками, а сама запись результата происходит в тот же массив, без создания отдельной структуры на каждый пиксель. Для оценки памяти можно прикинуть порядок величин. Цветной снимок 10980x10980 в CV_8UC3 занимает: 10980*10980*3 = 362мб. И это только одно изображение без временных массивов, исходных каналов и результата обработки. Если держать одновременно несколько Mat, три канала и копии под классификацию, память начинает расходоваться очень быстро.

Временные матрицы, которые больше не нужны, освобождаются явно:

Mat resized = new Mat();Cv2.Resize(bands[i], resized, new Size(targetCols, targetRows), 0, 0, InterpolationFlags.Linear);bands[i].Dispose();bands[i] = resized;

Для обычных маленьких картинок это выглядело бы занудством, но для спутниковых данных мы схлопочем OutOfMemory.

Где это может пригодиться

Такой проект можно воспринимать в двух смыслах.

Первый — образовательный. Если человек только входит в тему GIS и ДЗЗ, ему полезно руками собрать RGB-композицию из каналов, попробовать Sentinel-2, увидеть разницу между обычным и псевдоцветным отображением, выделить эталон и посмотреть, как меняется результат при разных метриках. После этого слова «спектральная сигнатура» и «канал NIR» перестают быть абстракцией.

Второй — исследовательский. Программа не заменяет большие GIS-пакеты, но позволяет быстро проверять идеи алгоритмов. Например:

какой критерий лучше выделяет воду на конкретном снимке;
как влияет порог спектрального угла на классификацию леса;
какие каналы лучше подходят для псевдоцветной композиции;
что даст сегментация при разных параметрах;
есть ли полезные признаки в первых EMD-модах.

Потенциальные прикладные сценарии тоже вполне понятны: мониторинг сельскохозяйственных культур, оценка площадей водных объектов, поиск выгоревших территорий, выделение песка или открытой почвы, предварительная разметка областей для дальнейшего анализа.

Самое приятное в таких задачах — результат видно сразу. Меняешь каналы, порог или эталон, запускаешь обработку и получаешь новую карту областей. Это сильно помогает понимать данные.

Что можно улучшить

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

Во-первых, можно добавить нормальную радиометрическую и атмосферную коррекцию. Сейчас программа в основном работает с теми значениями, которые прочитала из файлов, и масштабирует их для отображения.

Во-вторых, классификацию можно расширить до полноценной работы со спектральными сигнатурами по большему числу каналов. Сейчас основной сценарий завязан на выбранную трехканальную композицию, что удобно для визуализации, но не раскрывает весь потенциал гиперспектральных данных.

В-третьих, можно добавить экспорт не только картинки, но и геопривязанного результата. Для GIS-задач важно сохранять пространственную привязку, чтобы результат можно было открыть в QGIS или другом пакете.

В-четвертых, EMD явно просится на оптимизацию. Там много работы с каналами и пикселями, и текущая реализация скорее демонстрирует идею, чем выжимает максимум производительности.

Заключение

Для меня этот проект оказался хорошим способом потрогать дистанционное зондирование на уровне реальных данных и алгоритмов. У каждого снимка есть структура, каналы, метаданные, разрешение, формат хранения. У каждого алгоритма есть свои плюсы, ограничения и цена по времени.

В итоге получилось WPF-приложение, которое умеет читать данные Landsat 8, Sentinel-2 и AVIRIS, собирать изображения из каналов, считать растровые выражения, классифицировать снимки по эталонам, сегментировать их и экспериментировать с EMD.

Больше всего мне понравились две вещи. Первая — насколько выразительными становятся спутниковые данные, когда начинаешь менять каналы и смотреть за пределы обычного RGB. Вторая — барицентрическая классификация: необычный метод, который интересно ложится на задачу сравнения пикселей с эталоном и дает наглядный результат.

Если коротко, дистанционное зондирование — это область, где даже небольшое приложение быстро приводит к серьезным вопросам: как хранить большие данные, как сравнивать спектры, как отличать объект от фона, как не потерять смысл при визуализации. И именно этим она мне нравится.

Исходный код проекта, в описании есть ссылка на снимки для проб и ошибок: https://github.com/Alexey42/HSI

ссылка на оригинал статьи https://habr.com/ru/articles/1031848/