[Из песочницы] Детектирование поваленного леса на основе аэрофотосъёмки10.04.2019 19:03

Анализ / обзор

Первостепенная задача — выявление местоположения и характера завалов для поисково-спасательных отрядов. Независимо от постановки задач, необходимо планировать и отслеживать маршруты движения отряда в лесу. Неправильно направленный отряд может потеряться в буреломе от нескольких часов, до нескольких недель, вплоть до летального исхода. Завалы могут образоваться в любое время года, по самым разным причинам. Их нельзя отмечать на топографических картах, которые чертят раз в несколько десятков лет, как раз из-за их непостоянности. Завалы — временное, и неожиданно наступающее явление.
Другая задача, — исследование и инвентаризация лесов (лесотаксация), качественно решить которую возможно лишь автоматическими средствами. Это одновременно сканирование огромных площадей в поисках подходящих объектов на снимках аэрофотосъемки (АФС) и точный анализ выявленных сегментов снимков высокого разрешения.
На этой мозаике масштаб 1:15554 (или 15,5 км на сантиметр.):
Ниже — увеличенный участок этого снимка. Обратите внимание на выделенные фрагменты, которые «светятся» интенсивным голубым и синим цветом. Именно так выглядит сухостой и ветровал, то есть мертвые деревья.
Происходит это благодаря простому физическому эффекту: живые деревья всегда содержат хлорофилл, придающий листьям (или иголкам для сосен и елей) различные оттенки зеленого цвета. Причина — особенности отражения света хлорофиллом с поверхности растений в красном и ближнем инфракрасном диапазонах. Именно эти каналы используются при анализе растительного покрова и вычислении вегетационного индекса (NDVI).
Благодаря тому, что в мертвых деревьях этот пигмент отсутствует, они не отражают «красные» диапазоны. Именно такой особенностью и определяется интенсивный синий цвет бурелома на мультиспектральных снимках.
Поиск бурелома необходим для санитарных служб. Во многих странах проводится плановая уборка леса, в том числе и в России. Это необходимо для более «здорового» роста другой флоры.
Также, весь древесный уголь — продукт, полученный из сухостоя, то есть мертвых деревьев. Ведь в России, в зависимости от региона, налагается штраф за рубку зеленного леса. (Например, в Челябинской области на 2014 год за одно срубленное зеленное дерево полагался штраф в 8000 рублей).

Задача

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

Пример того, как может выглядеть решение, представлен на изображении ниже:

Исходные данные

Мозаика в формате GeoTIFF (1.5 Гб) находится по ссылке здесь (см первую каринку). Также снимок можно взять у эксперта.

Теоретическая часть:

Минимальная сегментация на основе связующего дерева
Сегментация изображения стремится разбить цифровое изображение на области пикселей со схожими свойствами, например однородностью. Представление области более высокого уровня упрощает задачи анализа изображений, такие как подсчет объектов или обнаружение изменений, поскольку атрибуты области (например, средняя интенсивность или форма) можно сравнить более легко, чем необработанные пиксели.
Мотивация для основанных на графике методов
Чтобы ускорить сегментацию больших изображений, работа может быть разделена между несколькими процессорами. Одним из способов достижения этого является разделение изображений на фрагменты, которые обрабатываются независимо. Однако области, которые охватывают границу фрагмента, могут быть разделены или потеряны, если фрагменты не соответствуют требованиям минимального размера алгоритма сегментации. Тривиальный обходной путь включает в себя перекрывающиеся фрагменты, то есть позволяет каждому процессору учитывать дополнительные пиксели вокруг границы своего фрагмента. К сожалению, это увеличивает вычислительную нагрузку, поскольку процессоры по обе стороны границы тайла выполняют избыточную работу. Кроме того, гарантированно сохраняются только объекты, меньшие, чем перекрытие мозаики, что означает, что длинные объекты, такие как реки на аэрофотоснимках, все еще могут быть разделены. В некоторых случаях результаты независимых плиток можно объединить, чтобы приблизить истинные результаты. [3] Альтернатива существует в форме основанных на графе методов сегментации. Информация о связях, присущая графам, позволяет выполнять независимую работу над частями исходного изображения и повторно соединять их, чтобы получить точный результат, как если бы обработка происходила глобально.
Возможность склеивания независимых подизображений мотивирует добавление информации о подключении к пикселям. Это можно рассматривать как график, узлами которого являются пиксели, а ребра представляют связи между пикселями. Простым и сравнительно экономичным вариантом этого является сеточный граф, в котором каждый пиксель связан со своими соседями в четырех основных направлениях. Поскольку отношение соседства пикселей является симметричным, результирующий график является ненаправленным, и необходимо сохранять только половину его краев (например, восточного и южного соседа каждого пикселя). Последний шаг требует кодирования информации о сходстве пикселей в весовых значениях ребер, чтобы исходное изображение больше не требовалось. В простейшем случае вес ребер вычисляется как разность интенсивностей пикселей.
Минимальное остовное дерево (MST) — это подмножество ребер графа графа с минимальным весом, так что все узлы соединены. В 2004 году Felzenszwalb представил метод сегментации, основанный на алгоритме MST Крускала. Края рассматриваются в порядке возрастания веса; их пиксели конечной точки объединяются в область, если это не вызывает цикл на графике, и если пиксели «похожи» на пиксели существующих областей. Обнаружение циклов возможно в почти постоянное время с помощью непересекающейся структуры данных. О сходстве пикселей судят по эвристике, которая сравнивает вес с порогом на сегмент. Алгоритм выводит несколько дизъюнктных MST, то есть лес; Каждое дерево соответствует сегменту. Сложность алгоритма квазилинейна, потому что сортировка ребер возможна за линейное время с помощью сортировки счетчиков
В 2009 году Вассенберг и соавт. разработал алгоритм, который вычисляет несколько независимых минимальных покрывающих лесов и затем соединяет их вместе. Это позволяет выполнять параллельную обработку без разделения объектов на границах листов. Вместо фиксированного порогового значения веса используется начальная маркировка подключенного компонента для оценки нижней границы порогового значения, что может снизить как избыточную, так и недостаточную сегментацию. Измерения показывают, что реализация на порядок превосходит последовательный алгоритм Felzenszwalb.

Скрипт (реализация):

import skimage.io
import io, json, codecs
from IPython.display import Image

# Read image into scimage package
img = skimage.io.imread(sample_data_file_name)
skimage.io.imsave('original.png', img)

# Display original image
display(Image(filename='original.png',width=800, height=800))

from osgeo import gdal, osr
import numpy
from skimage.segmentation import felzenszwalb
from skimage.segmentation import mark_boundaries
from skimage.measure import regionprops
import csv

# Prepare result structure
result = {
    "trees": []
}

# Open image with gdal
ds = gdal.Open(sample_data_file_name)
xoff, a, b, yoff, d, e = ds.GetGeoTransform()

# Get projection information from source image
ds_proj = ds.GetProjectionRef()
ds_srs = osr.SpatialReference(ds_proj)

# Get the source image's geographic coordinate system (the 'GEOGCS' node of ds_srs)
geogcs = ds_srs.CloneGeogCS()

# Set up a transformation between projected coordinates (x, y) & geographic coordinates (lat, lon)
transform = osr.CoordinateTransformation(ds_srs, geogcs)

# Convert multi-channel image it into red, green and blueb[, alpha] channels 
red, green, blue = numpy.rollaxis(numpy.array(img), axis=-1)
alpha = 0

# Mask: threshold + stops canny detecting image boundary edges
mask = blue > 240
mask1 = red < 90
mask2 = green < 200

# Create mask for edge detection
skimage.io.imsave('mask.png', mask * 255)
skimage.io.imsave('mask1.png', mask1 * 255)
skimage.io.imsave('mask2.png', mask2 * 255)


# Use Felzenszwalb algo to find segements
segments_fz = felzenszwalb(numpy.dstack((mask, mask1, mask2)),
                               scale=3700,
                               sigma=2,
                               min_size=4) 




# Build labeled mask to show where trees were dectected
segmented_img = mark_boundaries(numpy.dstack((mask, mask1, mask2)), segments_fz)
skimage.io.imsave('mask_labeled.png', segmented_img)
segmented_img1 = mark_boundaries(mask1, segments_fz)
skimage.io.imsave('mask_labeled1.png', segmented_img1)
segmented_img2 = mark_boundaries(mask2, segments_fz)
skimage.io.imsave('mask_labeled2.png', segmented_img2)

# Count trees and save image of each tree clipped from masked image
k=0

for idx, tree in enumerate(regionprops(segments_fz)):
    
    # If area matches that of a stanard tree, count it
    if (tree.area >= 2 and tree.area <= 100):
        
        # Incrment count
        k+=1
        
        # Create tree thumbnail
        x, y = (int(numpy.average([tree.bbox[0],
                                tree.bbox[2]])),
                                int(numpy.average([tree.bbox[1],
                                tree.bbox[3]])))
        sx, ex = max(x - 35, 0), min(x + 35, img.shape[0] - 1)
        sy, ey = max(y - 35, 0), min(y + 35, img.shape[1] - 1)
        img_tree = img[sx:ex, sy:ey]
        skimage.io.imsave('tree-id_' + str(idx + 1) + '-count_' + str(k) +'.png', img_tree)

        # Get global coordinates from pixel x, y coords
        projected_x = a * y + b * x + xoff
        projected_y = d * y + e * x + yoff
        
        # Transform from projected x, y to geographic lat, lng
        (lat, lng, elev) = transform.TransformPoint(projected_x, projected_y)
        
             
        
        # Add tree to results cluster
        result["trees"].append({
            "id": idx + 1,
            "tree_count": k,
            "lat": lat,
            "lng": lng
        })
    
Author: Nikita Permyakov

Profile in LinkedIn: /nikitapermikov

Выводы

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

Источники и ссылки:

Современные проблемы дистанционного зондирования Земли из космоса. 2016. Т. 13. № 1. С. 105–116. Некоторые приложения сегментации снимков ДЗЗ.Е. С. Иванов;
R. Haralick, L. Shapiro: Image Segmentation Techniques. CVGIP 29 (January 1985)
J. Iivarinen, M. Peura, J. Sarela, A. Visa: Comparison of combined shape descriptors for irregular objects. In: BMVC (1997) pp. 430–439
M.-H. Chen, T. Pavlidis: Image seaming for segmentation on parallel architecture. PAMI Vol. 12(6), June 1990, pp. 588–594
P. Felzenszwalb, D. Huttenlocher: Efficient Graph-Based Image Segmentation. IJCV 59(2) (September 2004)
G. Harfst, E. Reingold: A Potential-Based Amortized Analysis of the Union-Find Data Structure. SIGACT 31 (September 2000) pp. 86–95
J. Wassenberg, W. Middelmann, P. Sanders: An Efficient Parallel Algorithm for Graph-Based Image Segmentation. In: Computer Analysis of Images and Patterns, LNCS Vol. 5702 (September 2009) pp. 1003–1010