В белом-белом городе на белой-белой улице стояли белые-белые дома… А как быстро вы можете найти все крыши домов на этой фотографии?

Все чаще можно слышать про планы правительства провести полную инвентаризацию объектов недвижимости с целью уточнения кадастровых данных. Для первичного решения этой задачи можно применить простой способ, основанный на расчете площади крыш капитальных строений по аэрофотоснимкам и дальнейшее сопоставление с кадастровыми данными. К сожалению, ручной поиск и расчет занимает много времени, а поскольку новые дома сносятся и строятся непрерывно, то расчет требуется повторять снова и снова. Сразу возникает гипотеза, что этот процесс можно автоматизировать с помощью  алгоритмов машинного обучения, в частности, Computer Vision. В этой статье я расскажу о том, как мы в «НОРБИТ» ^[1] решали эту задачу и с какими сложностями столкнулись.

Спойлер — у нас получилось ^[2]. В основе разработанного ML сервиса находится модель глубокого машинного обучения, основанная на свёрточных (convolution) нейронных сетях. Сервис принимает на вход снимки с беспилотных летательных аппаратов, на выходе генерирует GeoJSON файл с разметкой найденных объектов капитального строительства с привязкой к географическим координатам.

В результате это выглядит так:

Проблемы

Начнем с технических проблем, с которыми мы столкнулись:

• есть существенное отличие между зимними и летними аэрофотоснимками (модель, обученная только на летних фотографиях совершенно не способна находить крыши зимой);
• для решения задачи определения площади капитального строения очень важно точное определение границ крыш, а не просто их обнаружение и локализация;
• часто возникала ситуация, когда несколько частей одной крыши оказывались на соседних тайлах (квадратных кропах облета), а поскольку процесс инференса (запуск обученной модели для новых изображений) выполняется для каждого тайла в отдельности, то на выходе вместо цельной крыши мы получали несколько ее кусков;
• модель демонстрирует хорошие результаты, начиная с определенного разрешения снимков, когда нейронная сеть может четко выделять контуры объектов (это важно учитывать при планировании аэрофотосъемки). Для инференса космоснимков нужно немного иначе обучать модель;
• некоторые объекты (например, бассейны и грузовики) модель иногда ошибочно считала крышами.

А еще беспилотники иногда приносят вот такие фотографии:

Также хотелось бы отметить и проблемы, которые могли бы быть, но нас не коснулись:

• у нас не было задачи выполнять инференс в ограниченное время (например, непосредственно в момент полета), что сразу решило все возможные проблемы с производительностью;
• на входе для обработки мы сразу получали качественные изображения в высоком разрешении (используются объективы с фокусным расстоянием 21мм на высоте 250 м, что составляет 5см/px) от нашего заказчика, компании «Шахты», могли использовать их экспертизу в геопозиционировании объектов на картах, а также имели возможность устанавливать определенный набор требований к будущим полетам беспилотников, что в итоге сильно снизило вероятность появления совсем уж уникальных тайлов, которых не было в обучающей выборке;

Первый вариант решения задачи, обводка с помощью Boundary box 

Несколько слов о том, какие инструменты были использованы нами  для создания решения.

• Anaconda — удобная система управления пакетами для языков Python и R.
• Tensorflow — открытая программная библиотека для машинного обучения, разработанная компанией Google.
• Keras — надстройка над фреймворками Deeplearning4j, TensorFlow и Theano.
• OpenCV — библиотека алгоритмов компьютерного зрения, обработки изображений и численных алгоритмов общего назначения с открытым кодом.
• Flask — фреймворк для создания веб-приложений на языке программирования Python.

В качестве ОС использовали Ubuntu 18.04. С драйверами на GPU (NVIDIA) в Ubuntu все в порядке, поэтому задача обычно решается одной командой:

> sudo apt install nvidia-cuda-toolkit

Подготовка тайлов

Первой задачей, с которой мы столкнулись, было разбиение изображений облета на тайлы (2048х2048 px). Можно было бы написать свой скрипт, но тогда бы пришлось думать о сохранении географической привязки каждого тайла. Проще было использовать готовое решение, например, GeoServer — это ПО с открытым исходным кодом, позволяющее публиковать на сервере геоданные. Помимо этого GeoServer решил для нас еще одну задачу – удобное отображение результата автоматической разметки на карте. Это можно сделать локально, например, и в qGIS, но для распределенной команды и демонстрации удобнее web-ресурс.

Для выполнения разбиения на тайлы нужно в настройках указать требуемый масштаб и размер.

Для переводов между системами координат мы использовали библиотеку pyproj:

from pyproj import Proj, transform

class Converter:
    P3857 = Proj(init='epsg:3857')
    P4326 = Proj(init='epsg:4326')
...
    def from_3857_to_GPS(self, point):
        x, y = point
        return transform(self.P3857, self.P4326, x, y)
    def from_GPS_to_3857(self, point):
        x, y = point
        return transform(self.P4326, self.P3857, x, y)
...

В результате удалось легко сформировать из всех полигонов один большой слой и наложить его сверху на подложку. 

GeoServer – не единственное приложение, позволяющее решить нашу задачу. В качестве альтернативы, например, можно рассмотреть ArcGIS for Server, однако этот продукт является проприетарным, потому использовать его мы не стали.

Далее предстояло на каждом тайле находить все видимые крыши. Первым подходом к решению задачи было использование object_detection ^[4] из набора models/research Tensorflow. Этим способом можно находить и локализовывать прямоугольным выделением (boundary box) классы на изображениях. 

Разметка обучающих данных 

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

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

Сейчас доступно множество сервисов, позволяющих произвести разметку изображений – как с открытым исходным кодом для установки на свой компьютер/сервер, так и корпоративные решения, включающие работу внешних асессоров, например Яндекс.Толока. В этом проекте мы использовали самый простой VGG Image Annotator ^[5]. В качестве альтернативы ему можно попробовать coco-annotator ^[6] или label-studio ^[7]. Последний мы обычно используем для разметки текста и аудиофайлов.

Для обучения на разметке различных аннотаторов обычно нужно выполнить небольшое перекладывание полей, пример для VGG ^[8].

Чтобы правильно вычислить площадь крыши, которая попала в область прямоугольного выделения, необходимо, соблюдения нескольких условий:

• все границы крыш должны быть параллельны/перпендикулярны на каждом тайле границам обводки. В противном случае будет проблема с превышением величины вычисленной площади:

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

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

...

lines = cv2.HoughLinesP(edges, 1, np.pi/180, 50, minLineLength=minLineLength, maxLineGap=5)
if lines is not None:
    length = image.shape[0]
    angles = []
    for x1, y1, x2, y2 in lines[0]:
        angle = math.degrees(math.atan2(y2 — y1, x2 - x1))
        angles.append(angle)
    parts_angles.append(angles)
    median_angle = np.median(angles)
...

# разбиение тайла на кропы 

for x in range(0, image.shape[0]-1, image.shape[0] // count_crops):
    for y in range(0, image.shape[1]-1, image.shape[1] // count_crops):
        get_line(image[x:x+image.shape[0]//count_crops, y:y+image.shape[1]//count_crops, :])
...

# вычисление медианы угла на всех кропах

np.median([a if a>0 else 90+a for a in np.array(parts_angles).flatten()])

После нахождения угла изображение поворачиваем с помощью аффинного преобразования:

h, w = image.shape[:2]
image_center = (w/2, h/2)

if size is None:
    radians = math.radians(angle)
    sin = math.sin(radians)
    cos = math.cos(radians)
    size = (int((h * abs(sin)) + (w * abs(cos))), int((h * abs(cos)) + (w * abs(sin))))
    rotation_matrix = cv2.getRotationMatrix2D(image_center, angle, 1)
    rotation_matrix[0, 2] += ((size[0] / 2) — image_center[0])
    rotation_matrix[1, 2] += ((size[1] / 2) — image_center[1])
else:
    rotation_matrix = cv2.getRotationMatrix2D(image_center, angle, 1)

cv2.warpAffine(image, rotation_matrix, size)

Полный код примера тут ^[9]. Вот как это выглядит:

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

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

Итоговое решение на базе Mask-RCNN

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

Материалов о принципах работы, видах и задачах свёрточных нейронных сетей, в том числе и на русском языке, уже предостаточно, поэтому мы не станем углубляться в них в этой статье. Остановимся только на одной конкретной реализации, Mask-RCNN – архитектуры для локализации и выделения контуров объектов на изображениях. Есть и другие отличные решения со своими достоинствами и недостатками, например, UNet, но лучшего качества удалось добиться именно на Mask-RCNN.

В процессе своего развития она прошла несколько стадий. Первая версия R-CNN была разработана в 2014 году. Принцип ее работы заключается в выделении на изображении небольших областей, для каждой из которых производится оценка вероятности наличия в этой области целевого объекта. R-CNN отлично справлялась с поставленной задачей, однако скорость ее работы оставляла желать лучшего. Закономерным развитием стали сети Fast R-CNN и Faster R-CNN, получившие улучшения в алгоритме обхода изображений, что позволило значительно повысить быстродействие. На выходе в Faster R-CNN появляется разметка прямоугольным выделением, обозначающим границы объекта, что не всегда достаточно для решения задачи. 

В Mask R-CNN добавляется еще и попиксельное наложение маски, позволяющее получить точное очертание объекта.

Наглядно boundary box и маски можно увидеть на результате работы модели (включен фильтр по минимальной площади здания):

Условно можно выделить 4 этапа в работе этой сети:

• стандартное для всех сверточных нейронных сетей выделение признаков на изображении, таких как линии, изгибы, контрастные границы и другие;
• Region Proposal Network (RPN) сканирует небольшие фрагменты изображения, именуемые anchors (якоря) и определяет, содержит ли этот якорь признаки, характерные для целевого класса (в нашем случае крыши);
• Region of Interest Classification and Bounding Box. На этом этапе сеть, основываясь на результатах предыдущего этапа пытается выделить на фотографии большие прямоугольные области, предположительно содержащие целевой объект;
• Segmentation Masks. На этом этапе из прямоугольной области, полученной при наложении boundary box, получается маска искомого объекта.

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

Использование исключительно только RGB изображений не позволило добиться необходимой точности распознавания (модель пропускала целые здания, была средняя ошибка в 15% в вычислении площади крыш), поэтому мы подавали на вход модели дополнительные полезные данные, например карты высот, полученные методом фотограмметрии. ^[10] 

Метрики, применяемые для оценки качества модели

При определении качества моделей мы чаще всего пользовались метрикой Intersection over Union (IoU)

Пример кода для вычисления IoU с использованием библиотеки geometry.shapely:

from shapely.geometry import Polygon

true_polygon = Polygon([(2, 2), (2, 6), (5, 6), (5, 2)])
predicted_polygon = Polygon([(3, 3), (3, 7), (6, 7), (6, 3)])
print(true_polygon.intersection(predicted_polygon).area / true_polygon.union(predicted_polygon).area)

>>> 0.3333333333333333

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

Tensorboard предоставляет данные по множеству различных метрик. Для нас наиболее интересными являются:

• val_mrcnn_bbox_loss — показывает, насколько хорошо модель локализует объекты (т. е. накладывает boundary box);
• val_mrcnn_mask_loss — показывает, насколько хорошо модель сегментирует объекты (т. е. накладывает маску).

Обучение модели и валидация

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

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

Создание фонового сервиса для распознавания

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

Ubuntu использует systemd, и пример будет приведен именно для этой системы. Код самого сервиса можно посмотреть здесь ^[11]. Пользовательские юниты находятся в директории /etc/systemd/system, там мы и создадим файл нашего сервиса.

cd /etc/systemd/system

sudo touch my_srv.service

Отредактируем файл:

sudo vim my_srv.service

Юнит systemd состоит из трех секций:

• [Unit] — описывает порядок и условие запуска (например, можно указать процессу дожидаться запуска определенной службы и только после этого стартовать самому);
• [Service] — описывает параметры запуска;
• [Install] — описывает поведение сервиса при добавлении его в автозагрузку.

В результате наш файл будет выглядеть таким образом:

[Unit]
Description=my_test_unit

[Service]
WorkingDirectory=/home/user/test_project
User=root
ExecStart=/home/user/test_project/venv/bin/python3 /home/user/test_project/script.py

[Install]
WantedBy=multi-user.target

Теперь перезагрузим конфигурацию systemd и запустим наш сервис:

sudo systemctl daemon-reload
sudo systemctl start my_srv.service

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

Выводы

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

В результате оценки точности модели на тестовых данных получены следующие значения: количество найденных крыш – 91%, точность обводки крыш полигонами – 94%.

Удалось добиться приемлемого качества работы моделей на летних и зимних облётах, но качество распознавания может снижаться на снимках сразу после снегопада.

Теперь даже Сиднейский оперный театр не ускользнет от глаз нашей модели. 

Мы планируем выложить этот сервис с обученной моделью на нашем демостенде. Если интересно попробовать работу сервиса на своих собственных фото, присылайте заявки на ai@norbit.ru.

Автор: Дмитрий Тимаков

Источник ^[12]