Геокоррекция

Цель решения задачи — находить в реальном времени местоположение и ориентацию текущего снимка на подложке.

1. Подложка может представлять собой изображение со спутника или беспилотного летательного аппарата.
2. Подложка может быть снята в другое время дня, время года, и на другой высоте (обычно высота подложки больше высоты текущей картинки).
3. Направление камеры при съемке вертикально вниз.
4. Стек технологий допустимы в решении задачи: a. ОС - Linux like b. Язык программирования - Python/C++ c. GUI - не требуется
5. Скорость обработки - в реальном времени
6. Требования к железу - без использования GPU(опционально)

Нелинейность оптики

Снимок с любой камеры имеет искажения, в центре объекты практически не растягиваются, а по краям растяжение минимально. Чтобы по искажению восстановить "равномерный" снимок, чаще всего используют полиномиальную коррекцию. Пример её использования описан в этой статье. Само восстановление называется ортотрансформированием или ортокоррекцией.

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

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

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

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

Опорные объекты

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

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

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

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

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

Если преобразование по контуру найдено, то можно делать дополнительную проверку, удовлетворяют ли преобразованию все точки объекта.

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

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

Минимальное отсечение

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

Для решения этой задачи можно реализовать подобие алгоритма отсечения. Каждая точка на изображении камеры должна соответствовать ровно одной точке на подложке.

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

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

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

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

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

Поиск афинного преобразования

Любое афинное преобразование в 2D описывается матрицей 3x3:

$$A = \begin{bmatrix} R_{11} & R_{12} & t_x\\ R_{21} & R_{22} & t_y\\ 0&0&1 \end{bmatrix}$$

$R_{ij}$ - компоненты совмещённой матрицы поворота и масштаба, а $t_x$ и $t_y$ - смещения вдоль соответствующих осей

Обозначим координаты $i$-ой опорной точки на подложке как $p_i$, а соответствующую ей точку на изображении с камеры $v_i$, тогда они должны удовлетворять уравнению

$$A_tp_i=v_i$$

где $A_t$-матрица преобразования из системы координат подложки в систему координат камеры.

Т.к. в матрице $A_t$ неизвестны только шесть элементов и каждое из равенств является линейной комбинацией, то для однозначного вычисления матрицы преобразования необходимо составить шесть уравнений. Для этого нам понадобится три опорных точки на подложки и три соответствующие им на экране камеры. Каждая пара даст два уравнения: по x и по y координатам.

Запишем эти уравнения в матричном виде:

$$A_t\begin{bmatrix} p_{1_x}&p_{2_x}&p_{3_x}\\ p_{1_y}&p_{2_y}&p_{3_y} \\ 1&1&1 \end{bmatrix}=\begin{bmatrix} v_{1_x}&v_{2_x}&v_{3_x}\\ v_{1_y}&v_{2_y}&v_{3_y}\\ 1&1&1 \end{bmatrix}$$

Обозначим соответствующие матрицы значений векторов как $P$ и $V$ и перепишем предыдущее уравнение:

$$A_tP=V$$

Если определитель $\det P$ далёк от нуля (матрица далека от сингулярной), то можно просто умножить обе части на $P^{-1}$:

$$A_t=VP^{-1}$$

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

$$A_tPP^T=VP^T$$

$$A_t(PP^T)(PP^T)^{-1}=VP^T(PP^T)^{-1}$$

$$A_t=VP^T(PP^T)^{-1}$$

Проверка нарушения преобразования

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

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

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

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

Восстановление положения

Чтобы восстановить положение, угол поворота и масштаб, можно использовать формулы отсюда:

$$\left[\begin{array}{ccc} \mathrm{a} & \mathrm{b} & \mathrm{tx}\\ \mathrm{c} & \mathrm{d} & \mathrm{ty}\end{array}\right]=\left[\begin{array}{ccc} s_{x}\cos\psi & -s_{x}\sin\psi & x_{c}\\ s_{y}\sin\psi & s_{y}\cos\psi & y_{c}\end{array}\right]$$

где $s_x$ , $s_y$ коэффициенты масштабирования вдоль соответствующих осей, $x_c$ , $y_c$ - смещения, $\psi$ - угол поворота

Тогда:

$$x_{c}=\mathrm{tx}$$

$$y_{c}=\mathrm{ty}$$

$$s_{x}=\mathrm{sign(a)\,}\sqrt{\mathrm{a}^{2}+\mathrm{b}^{2}}$$

$$s_{y}=\mathrm{sign(d)\,}\sqrt{\mathrm{c}^{2}+\mathrm{d}^{2}}$$

$$\tan\psi=-\frac{\mathrm{b}}{\mathrm{a}}=\frac{\mathrm{c}}{\mathrm{d}}$$

Тогда угол поворота

$$\psi = {\rm atan2}(-b,a) = {\rm atan2}(c,d)$$