Что такое цифровая обработка изображений? Зачем нам вообще знать про алгоритмы обработки, когда есть фотошоп и фильтры в телефоне? Или всё можно отдать нейросети и получить крутой результат? И при чём тут Julia, наконец? Будем разбираться!
Мы запускаем серию статей про обработку изображений с использованием языка Julia и вычислительной среды Engee. Задача – ответить на часто встречающиеся вопросы вроде
актуальности этого направления компьютерной науки;
задач, решаемых методами обработки изображений;
применения и реализации стандартных и «умных» алгоритмов.
В первой части ознакомимся с основами на примере сегментации спутникового снимка. Постараюсь выдержать баланс теории и практической реализации в коде 😉
Кстати, тематику передачи изображений а также тематики ЦОС, связи и радиолокации мы будем обсуждать на семинаре День Engee, который пройдёт уже через неделю.
Введение (он же ликбез)
Цифровая обработка изображений (ЦОИ) в понимании рядового юзера – это работа с картинками на компьютере с помощью математических алгоритмов. Если представить изображение как матрицу значений интенсивности его пикселей в каналах, то ЦОИ позволяет менять цвет, яркость, контраст, удалять шумы, находить объекты, повышать качество картинки.
Впрочем, для работы с личными фото погружаться в премудрости не требуется – огромное количество доступных редакторов в пару кликов сделают всё за нас. Но применение алгоритмов ЦОИ простирается гораздо дальше фильтров для любимых censored-грамов.
Количество обрабатываемой фото- и видеоинформации ещё никогда не было таким большим. Как вы, наверное, догадались, обработка основной части этой информации происходит автоматически, без участия оператора. Вот краткий список приложений ЦОИ:
Медицина и биология: улучшение рентгеновских снимков, МРТ, КТ, УЗИ, обнаружение опухолей, микроскопия клеток, отслеживание движения крови в сосудах.
Компьютерное зрение и робототехника: распознавание лиц, жестов, автомобилей, навигация роботов, дополненная реальность.
Дистанционное зондирование и геоинформатика: спутниковые и аэрофотоснимки – анализ состояния лесов, сельхозугодий, городской застройки, картография, мониторинг стихийных бедствий.
Промышленность и автоматизация: обнаружение дефектов на производстве (микрочипы, автомобили, пищевая промышленность), оптическая сортировка, роботизированная сборка.
Безопасность и видеонаблюдение, развлечения и мультимедиа, искусственный интеллект и машинное обучение, оцифровка произведений искусства и многое другое.
А если говорить о платформах, на которых алгоритмы разворачиваются, то тут весь спектр, начиная от десктопных приложений, продолжая суперкомпьютерами и кластерами GPU, заканчивая ПЛИС и компактными ASIC для обработки видеопотока на борту БПЛА.
Мы же в свою очередь для разработки используем уже отлично зарекомендовавшую себя среду Engee. Она постоянно развивается, не требуется установки, так как работает в браузере, поддерживает языки Julia и Python. А еще в ней можно применять GPU для расчетов.
Постановка задачи
Охватить всё перечисленное выше не представляется возможным, поэтому будем начинать знакомиться с алгоритмами на конкретном примере – сегментации изображения спутникового снимка.
Для чего обрабатывать спутниковые снимки, полагаю, понятно: картография, автоматизированный контроль для сельского хозяйства, экологии, безопасности, планирования застройки, обнаружение интересующих объектов, наконец. Но что такое сегментация?
Сегментация изображений – это разделение картинки на осмысленные части или объекты. Например, на изображении человека можно выделить лицо, волосы, одежду, фон и т. д.
В медицине: отделить опухоль от здоровых тканей на МРТ.
В беспилотных авто: найти дорогу, пешеходов, другие машины.
В фотографии: заменить фон (как в телеконференциях или редакторе).
В сельском хозяйстве: анализировать состояние растений по снимкам с дронов.
В результате применения алгоритма сегментации каждому пикселю цифрового изображения присваивается «лейбл» – знак принадлежности к тому или иному классу. Если в качестве «лейблов» использовать целые числа, то мы получим матрицу размером с исходное изображение, содержащую значения этих чисел на месте пикселей объектов или областей. Это так называемая семантическая сегментация.
Мы в примере будем рассматривать что-то подобное, но для частного случая, когда классов всего два. Эту задачу можно назвать выделением объекта на фоне, а, значит, у нас будет только классы «объект» и «фон». Результирующая матрица будет состоять из чисел 1 и 0 соответственно.
В первой части нашего увлекательного путешествия в мир ЦОИ для сегментации изображения мы рассмотрим относительно простой алгоритм, состоящий из этапов предобработки, бинаризации, морфологии и наложения результирующей матрицы (маски).
В этих этапах подробнее разберёмся по ходу решения. В качестве тестового изображения я взял спутниковый снимок Сокольнического парка и прилегающих жилых районов с целью отделить зелёный массив от построек. Что нужно сделать, понятно, теперь узнаем – как.
Разработка на языке Julia в Engee
Во-первых, почему Julia: я долгое время работал в MATLAB и обработкой изображений тоже занимался в MATLAB. Иногда доходило до того, что для какой-то стандартной задачи, вроде изменения контрастности или удаления шума на фотке, я открывал не фотошоп, а садился писать скрипт в знакомой среде. Благо, синтаксис позволял этапы обработки описывать одной строчкой кода. Ну а Julia и её расширения для ЦОИ максимально похожи на MATLAB, поэтому для меня порог вхождения был совсем низким.
Об основных особенностях языка писали здесь, а вот с точки зрения нашей основной темы стоит всё же пояснить, почему не обратиться к стандарту - Python с OpenCV и PyTorch.
Помимо матлабоподобного синтаксиса, Julia быстрее вышеупомянутого стандарта в 10-100 раз на нативных операциях с матрицами, имеет встроенную поддержку многопоточности без GIL, а JIT-компиляция и специализация типов делают её идеальной для задач, где важна скорость: обработка видео, 3D-реконструкция, нейросетевые модели реального времени. Для задач ЦОИ и компьютерного зрения используем:
Пакеты для базовой обработки
Images.jl – аналог OpenCV (фильтры, морфология, цветовые пространства),
ImageSegmentation.jl – классические методы (k-means, watershed),
ImageFeatures.jl – детекция ключевых точек (SIFT, ORB).
Глубокое обучение
Flux.jl – нейросетевой фреймворк с автоматическим дифференцированием,
Metalhead.jl – предобученные модели (ResNet, VGG),
Поддержка ONNX и PyTorch через ONNX.jl и Torch.jl.
Уникальные фичи
Бенчмаркинг в одну строку (@btime из BenchmarkTools),
Макросы для генерации SIMD-оптимизированного кода.
Ну а пользоваться скоростью и удобством Julia я буду в среде разработки Engee, а именно, в её интерактивных скриптах для технических расчётов. Я думаю, вы убедитесь в том, что писать код и выводить промежуточные результаты обработки тут так же удобно, как и в MATLAB.
Загрузка и визуализация
Наша программа (скрипт) будет начинаться с указания используемых библиотек. Используем библиотеки для общей ЦОИ, визуализации, изменения контраста, бинаризации и морфологии:
using Images, ImageShow, ImageContrastAdjustment, ImageBinarization, ImageMorphologyТеперь вызовем функцию load для загрузки исходного цветного изображения. Точку с запятой в конце строчки не ставим и тут же наблюдаем вывод. Напомню, мы обрабатываем область снимка района Москвы с застройкой и лесополосой:
I = load("$(@__DIR__)/map_small.jpg") # загрузка исходного цветного изображения
И вот перед нами встаёт вопрос: по какому принципу отнести пиксели к классу «фона» или «объекта»? Они должны как-то явно отличаться, но у пикселей не так много параметров. В нашем случае это только интенсивность в цветовых каналах. Давайте же рассмотрим сами каналы:
CV = channelview(I); # разложение изображения на каналы
[ RGB.(CV[1,:,:], 0.0, 0.0) RGB.(0.0, CV[2,:,:], 0.0) RGB.(0.0, 0.0, CV[3,:,:]) ]
В каждом из каналов находится матрица, соответствующая изображению в градациях серого (цвет на визуализации добавили просто для наглядности). Можно заметить, что зелёные насаждения выглядят темнее именно в синем канале (а не в зелёном 😉). Убедимся в этом, сравнив матрицу в синем канале с монохромным изображением (объединением интенсивностей трёх каналов), которое я получаю функцией Gray.(I):
Интенсивность «пикселей леса» в синем канале явно ниже, чем «пикселей города», и мы можем попробовать разделить их по порогу, то есть бинаризовать. Но не всё так просто – изображение содержит весьма мелкие детали, и при бинаризации «в лоб» результирующая маска будет содержать много шума. Поэтому нам нужна предобработка. И да, мы не будем осуществлять пока что сегментацию по цвету – продолжим работать с синим каналом, то есть с изображением в градациях серого.
Предобработка
В задачи предобработки входят изменение контраста и яркости, удаление шума, сглаживание текстуры, и проч. В нашем случае полезно будет по изображению пройтись сглаживающим фильтром: на тёмные области леса эта операция особо не повлияет, а вот город она уравняет по яркости.
Мы будем применять двумерную фильтрацию, то есть свёртку исходного изображения с матрицей меньшего размера, так называемым кернелом. Кстати, свёртка изображений – линейная операция, используемая в нейросетевых алгоритмах обработки изображений. Но пока говорим только о линейной фильтрации – умножить, сложить, передвинуть.
Сравним два кернела – с двумерным распределением Гаусса и равномерный, выполняющий, по сути, операцию нахождения среднего арифметического. Гауссов кернел у нас вышел размером 9х9 пикселей, а осредняющий – 7х7.
Сумма весов в кернеле равна 1, поэтому подобные сглаживающие фильтры не изменяют общую интенсивность изображения. Выполним двумерную фильтрацию функцией imfilter:
BLUE = CV[3,:,:]; # интенсивность синего канала RGB
gaussIMG = imfilter(BLUE, Kernel.gaussian(2)); # линейный фильтр Гаусса
kernel_average = ones(7,7) .* 1/(7^2); # кернел среднего арифметического
avgIMG = imfilter(BLUE, kernel_average); # линейный фильтр среднего
Мне больше нравится вариант справа, поэтому дальнейшая обработка будет проводится над результатом фильтрации по методу осреднения – матрицей avgIMG. Теперь наша задача отнести более светлые пиксели к классу «город», а более тёмные – к классу «лес».
Бинаризация
Бинаризация – это процесс формирования матрицы, содержащий бинарные значения, то есть 1 или 0. На изображении это соответствует белым и чёрным пикселям. Обычно чёрные пиксели относят «фону», а белые – «объектам» (но не обязательно так).
Самый простой способ бинаризовать матрицу интенсивности – сравнить с пороговым значением. Всё, что не превышает порог, станет нулями, а что превышает – единицами.
В библиотеке ImageBinarization.jl собраны более сложные методы, учитывающие неравномерный фон, изменения контраста изображения и прочее. Мы сравним метод Отсу, доступный посредством функции binarize, и простое сопоставление интенсивности пороговому значению.
binOTSU = binarize(avgIMG, Otsu()); # бинаризация методом Отсу
binTHRESH = avgIMG .> 0.25; # бинаризация по порогу 0.25
Результаты почти идентичны, и мы продолжим с бинарным изображением после простого сравнения с порогом. Вроде бы мы уже получили бинарную матрицу (маску), но на ней много шума: области имеют разрывы, «прорези», присутствуют мелкие объекты. Чтобы сделать маску более «гладкой» и свободной от шума, мы воспользуемся алгоритмами бинарной морфологии.
Морфологические операции
Это простые, но мощные методы обработки бинарных или полутоновых изображений, основанные на геометрии объектов. Они помогают убирать шум (мелкие артефакты), склеивать или разделять объекты, улучшать форму контуров.
Типичные операции: дилатация, раскрытие, закрытие и эрозия. Эти операции работают со структурным элементом (ядром) – небольшой матрицей (например, 3×3), которая «сканирует» изображение и изменяет пиксели по заданным правилам.
Подробнее про морфологию можно почитать тут или тут.
Нам доступны параметризованные функции для создания типичных структурных элементов типа «ромб» или «квадрат», впрочем, нам интересно создать круглый структурный элемент. Поэтому мы бинаризуем двумерное распределение Гаусса по порогу в 0.0025.
kernel_d = strel_diamond((13, 13), r=6); # структурный элемент "ромб"
kernel_b = strel_box((13, 13), r=4); # структурный элемент "квадрат"
se = Kernel.gaussian(3) .> 0.0025; # структурный элемент "пользовательский"
Теперь осуществим операцию морфологического закрытия функцией closing с пользовательским структурным элементом. Это позволит «закрыть» небольшие отверстия и щели между белыми пикселями, а также сгладить форму объектов. Объектами (или блобами) мы называем «объединённые» области белых пикселей, окружённых черными пикселями, или границей изображения. После операции закрытия мы воспользуемся функцией area_open и удалим объекты, размер которых не превышает 500 пикселей.
closeBW = closing(BW_AVG,se); # морфологическое закрытие
openBW = area_opening(closeBW; min_area = 500) .> 0; # удаление "малых" объектов
Теперь инвертируем нашу бинарную маску и удалим объекты размером менее 5000 пикселей. Инвертированную маску также сгладим операцией закрытия:
fillBW = area_opening(.!openBW; min_area = 5000); # удаление "средних" объектов
MASK_AVG = closing(fillBW,se); # морфологическое закрытие (сглаживание)
Результат сегментации
Я хочу визуально оценить, насколько хорошо полученная бинарная маска сегментирует исходное изображение. Подсветить объекты на фоне можно, добавив 30% интенсивности бинарной маски в красный и зелёный каналы исходного цветного изображения. Таким образом, получим полупрозрачную маску приятного жёлтого цвета:
sv1_AVG = StackedView(CV[1,:,:] + (MASK_AVG./3), CV[2,:,:] +
(MASK_AVG./3), CV[3,:,:]); # наложение "жёлтой маски"
forest_AVG = colorview(RGB, sv1_AVG);
sv2_AVG = StackedView(CV[1,:,:] + (.!MASK_AVG./3), CV[2,:,:] +
(.!MASK_AVG./3), CV[3,:,:]); # наложение "жёлтой маски"
city_AVG = colorview(RGB, sv2_AVG);
Моя оценка – неплохо. Мы наблюдаем выделение лесного массива. Впрочем, скрипт у нас получился совсем не универсальным, так как многие параметры подбирались не адаптивно, а вручную по конкретно этому тестовому изображению. Но мы и не задавались целью создать робастный «боевой» алгоритм для разнообразных случаев. Какие ещё подходы можно использовать?
Нелинейная фильтрация
Предобработка, которую мы осуществляли ранее, использовала линейные фильтры, но мы вправе внутри выбранной области «под кернелом» делать всё, что не запрещено математикой. Например, перемешивать элементы массива, находить максимумы, минимумы, медианные значения, градиенты и осуществлять многое-многое другое, что можно описать ёмким термином «нелинейная фильтрация».
Я вижу, что области леса на изображении более однородные, без резких изменений яркости, в то время, как город гораздо более «пёстрый». Нелинейным фильтром можно, например, найти среднеквадратическое отклонение внутри небольшой области изображения, энтропию или абсолютное значение разницы по столбцам матрицы. Для областей, где интенсивность меняется слабо, численное значение выхода нелинейного фильтра будет меньше. Подобные фильтры ещё называют текстурными, потому что они позволяют разделять области с одинаковой средней интенсивностью, но отличающиеся по текстуре.
Один из самых простых подобных фильтров – фильтр диапазона. Он возвращает разницу между максимальным и минимальным значением интенсивности внутри кернела. Его я и буду использовать, но сначала надо сделать изображение более контрастным. Я возьму монохромное изображение и применю к нему функцию растяжения гистограммы:
IMG = adjust_histogram(Gray.(I), LinearStretching(dst_minval = -1, dst_maxval = 1.8));
Пришло время воспользоваться фильтром диапазона. В MATLAB я бы просто вызвал функцию rangefilt, но в доступных пакетах Julia я такой не нашёл…
С одной стороны, придётся покопаться в документации, с другой – чуть-чуть покодить тоже полезно. В пакетe ImageFiltering.jl есть контейнер для пользовательских функций для операций над пикселями внутри окна (кернела) – mapwindow. Напишу для неё функцию расчёта значения диапазона:
function range(array) # функция нелинейной фильтрации
(minval, maxval) = extrema(array); # мин. и макс. значения массива
return maxval - minval; # диапазон (макс. - мин.)
endТеперь применю mapwindow c моей функцией range к контрастному изображению. А выход нелинейного фильтра уже бинаризуем методом Отсу:
rangeIMG = mapwindow(range, IMG, (7,7)); # функция `range` в окне 7х7
BW_RNG = binarize(rangeIMG, Otsu()); # бинаризация методом Отсу
Соберём последующие за бинаризацией операции в одну функцию myMorphology и применим её к бинарному изображению:
function myMorphology(BW)
se = Kernel.gaussian(3) .> 0.0025;
closeBW = closing(BW,se);
openBW = area_opening(closeBW; min_area = 500) .> 0;
fillBW = area_opening(.!openBW; min_area = 5000)
MASK = closing(fillBW,se);
return MASK
end
Собираем пользовательскую функцию сегментации
Изучив различные методы и их параметры, можно последовательность применяемых команд зафиксировать и собрать в одну функцию, которая будет получать на вход цветное изображение, а на выходе – давать бинарную маску. За основу возьмём предобработку с растяжением гистограммы и нелинейной фильтрацией:
function mySegRange(I)
CV = channelview(I);
IMG = adjust_histogram(Gray.(I), LinearStretching(dst_minval = -1, dst_maxval = 1.8));
rangeIMG = mapwindow(range, IMG, (7,7));
BW_RNG = binarize(rangeIMG, Otsu());
MASK_RNG = myMorphology(BW_RNG);
return MASK_RNG
endПрименим функцию к более крупному снимку местности:
I2 = load("$(@__DIR__)/map_big.jpg");
CV2 = channelview(I2);
myMASK = mySegRange(I2);
sview = StackedView(CV2[1,:,:] + (.!myMASK./3), CV2[2,:,:] + (.!myMASK./3), CV2[3,:,:]);
colorview(RGB, sview)
Практический смысл подобных функций – возможность осуществлять анализ и вычисления над выделенными объектами в автоматическом режиме. Зная линейные размеры пикселя на карте в метрах, можно, к примеру, определить линейные размеры и площадь объектов на карте. В нашем случае мы можем подсчитать площадь лесополосы, застройки, узнать их отношение и т.д.
Заключение
Итак, мы рассмотрели базовые операции цифровой обработки изображений на примере процесса сегментации спутникового снимка. При этом затронули как теоретические аспекты – последовательность типовых операций, так и аспекты их практической реализации на языке Julia.
Отдельно отмечу приятный опыт работы в среде Engee: реализация скрипта заняла не так много времени, результаты сохранил в своём облаке, а также быстро закинул в Сообщество. Если у вас есть желание самостоятельно пробежаться по скрипту и попробовать его модифицировать для обработки ваших изображений, то весь код доступен по ссылке.
В следующей части мы рассмотрим работу с цветом и «умные» алгоритмы сегментации. До новых встреч!
