Краткое руководство по профилированию линии горизонта городской панорамы с помощью Python в несколько строк кода

Позвольте мне сказать очевидное: линии горизонта — это красиво.

С самого детства я мечтал жить в США, и одной из причин для этого была идея находиться среди прекрасных панорам Нью-Йорка, Чикаго или Сан-Франциско.

Когда я вырос (и стал жить в США), то понял, что Штаты могут предложить гораздо больше, чем архитектурные очертания на фоне неба, но я по-прежнему их люблю.

В этой статье мы узнаем, как извлечь профиль линии горизонта из фотографии. Примерно так:

Рисунок (1): Изображение автора, данные отсюда (Лицензия: CC0: Public Domain)

Давайте начнем.

0. Идея

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

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

Тогда почему определить небо проще, чем небоскребы?

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

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

Я хочу доказать свои слова на примере. Возьмем эту фотографию.

Рисунок (2): Изображение автора, данные отсюда (Лицензия: CC0: Public Domain)

Отлично. Теперь вырежем фрагмент из этой области:

Рисунок (3): Изображение автора, данные отсюда (Лицензия: CC0: Public Domain)

Затем возьмем отрезок от 0 до 50 и выведем стандартное отклонение:

Рисунок (4): Изображение автора, данные отсюда (Лицензия: CC0: Public Domain)

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

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

1. Алгоритм

Алгоритм для выполнения того, что вы видите на Рисунке (1), следующий:

Рисунок (5): Изображение автора

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

1.1 Преобразование изображения в черно-белое

Я понимаю, что вы и так все это знаете. Но важно сказать, зачем мы это делаем. Как вы понимаете, все этапы по блюрингу и фильтрации имеют смысл, когда вы применяете их к матрице. Цветное изображение технически является тензором, потому что оно имеет число строк X, число столбцов X 3, значения каналов (красный зеленый и синий). Черно-белое изображение — это матрица, состоящая из числа строк и числа столбцов.

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

1.2 Этапы блюринга

Оба этих этапа — применение медианного и нормализованного бокс-фильтра — используются для удаления шумов из сигналов с сохранением границ.

1.3 Фильтр Лапласа

Главным героем этого проекта, конечно же, является фильтр Лапласа. Этот фильтр считается второй производной дискретного пространства по времени.

Зачем нам вообще нужна вторая производная по времени?

Мы уже говорили, что между небом и зданием существует разница в стандартном отклонении. Это изменение стандартного отклонения возникает в определенной точке, которая является границей изображения (и небоскреба).

Следовательно, мы хотим увидеть быстрое изменение изображения. В частности, нужно, чтобы такое изменение было максимальным. Это означает, что нам нужна точка (или соседние точки), где вторая производная по времени равна нулю.

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

Вторая производная по времени выглядит следующим образом:

Это ядро (фильтр), которое мы собираемся применить к изображению, и оно даст нам картинку для второй производной.

1.4 Применение порога 1/0

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

Вот почему мы применяем порог 1/0.

1.5 Фильтр Erode

Фильтр Erode (англ. — подвергать эрозии, размывать) — это то, что мы используем для сглаживания имеющегося изображения. Идея, которая стоит за этим, заключается в том, что мы хотим сделать изображение более четким. Говоря технически, есть ядро, которое “проходит” по изображению и заменяет значения на минимальные. Опять же, поскольку сейчас у нас изображение 1/0, это просто делает его более контрастным.

1.6 Установите маску на 0, пока не будет найден последний индекс

Этот этап достаточно сложен для объяснения, но его очень легко понять. Когда вы все это проделаете, возможно, в колонке вашего изображения будет последовательность из 0 и 1. Это не имеет особого смысла, поскольку не может быть что-то вроде "небоскреб - небо - небоскреб снова". По этой причине мы находим наш самый большой индекс со значением 0 в нашей колонке и устанавливаем все равным 0, пока не найдем это значение. Тогда все остальное равно 0.

2. Практическая реализация

Объяснять это было весьма долго, но реализовать легко.

Давайте сделаем это пошагово:

2.1 Импорт библиотек:

import numpy as np 
import pandas as pd 
from os import listdir
from PIL import Image
import matplotlib.pyplot as plt
from os.path import isfile, join
import cv2
from scipy.signal import medfilt
from scipy import ndimage
import numpy as np
from matplotlib import pyplot as plt
import os
from sklearn.metrics import plot_confusion_matrix
from sklearn.metrics import confusion_matrix

2.1 Импорт данных:

Я получил данные с сайта Kaggle. Набор данных является открытым и не имеет авторских прав (CC0: Public Domain). В частности, мною была загружена только часть датасета с изображениями 12 городов по 10 небоскребов в каждом:

mypath = 'data/images'
subfolders = [ f.path for f in os.scandir(mypath) if f.is_dir() ]
images = []
images_baw = []
labels = []
label=0
size = 256
string_labels = []
for folder in subfolders: 
    onlyfiles = [f for f in listdir(folder) if isfile(join(folder, f))]
    for file in onlyfiles:
        image_file = Image.open(folder+'/'+file).resize((size,size))
        images.append(np.array(image_file))
        images_baw.append(np.array(image_file.convert('1')))
        labels.append(label)
        string_labels.append(folder.split('/')[-1])
    label=label+1
labels = np.array(labels)
images = np.array(images)
image_baw = np.array(images_baw)

2.2 Визуализация данных

plt.figure(figsize=(20,20))
for i in range(1,13):
    plt.subplot(4,3,i)
    identify_label = np.where(labels==i-1)[0]
    identify_label = np.random.choice(identify_label)
    plt.title('City label = %s'%(string_labels[identify_label]))
    plt.imshow(images[identify_label])

2.3 Определение функций:

Вся теория, изложенная выше, реализуется в виде следующих строк:

def cal_skyline(mask):
    h, w = mask.shape
    for i in range(w):
        raw = mask[:, i]
        after_median = medfilt(raw, 19)
        try:
            first_zero_index = np.where(after_median == 0)[0][0]
            first_one_index = np.where(after_median == 1)[0][0]
            if first_zero_index > 20:
                mask[first_one_index:first_zero_index, i] = 1
                mask[first_zero_index:, i] = 0
                mask[:first_one_index, i] = 0
        except:
            continue
    return mask


def get_sky_region_gradient(img):

    h, w, _ = img.shape

    img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

    img_gray = cv2.blur(img_gray, (9, 3))
    cv2.medianBlur(img_gray, 5)
    lap = cv2.Laplacian(img_gray, cv2.CV_8U)
    gradient_mask = (lap < 6).astype(np.uint8)

    kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (9, 3))

    mask = cv2.morphologyEx(gradient_mask, cv2.MORPH_ERODE, kernel)
    # plt.imshow(mask)
    # plt.show()
    mask = cal_skyline(mask)
    after_img = cv2.bitwise_and(img, img, mask=mask)

    return after_img

2.4 Применение алгоритма

Это применение алгоритма для всего датасета:

zero_images = []
zero_images_plot = []
for K in range(len(images)):
    zero_image_plot = np.zeros((size,size))
    zero_image = zero_image_plot.copy()
    for i in range(size):
        zero_image_plot[i,tot_profiles[K][i]-2:tot_profiles[K][i]+2] = 1
        zero_image[i,tot_profiles[K][i]] = 1 
    zero_images_plot.append(zero_image_plot.T)
    zero_images.append(zero_image.T)

Вот пример:

K=0
plt.figure(figsize=(10,10))
plt.suptitle('City = %s'%(string_labels[K]),fontsize=20,y=0.72)
plt.subplot(1,3,1)
plt.title('Original Image')
plt.imshow(images[K])
plt.subplot(1,3,3)
plt.title('Masked Image')
plt.imshow(zero_images_plot[0])
plt.subplot(1,3,2)
plt.title('Profiled Skyline')
plt.imshow(get_sky_region_gradient(images[K]))

<matplotlib.image.AxesImage at 0x7f987f0ba520>

Довольно круто, правда?

2.5 Преобразование в сигнал

Для того чтобы получить правильный облик Рисунка (1), мы применим следующую функцию:

def signal_from_profile(K):
    x = np.arange(size)
    y = np.array(size-np.array(tot_profiles[K]))
    return x,y

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

plt.figure(figsize=(13,10))
#plt.suptitle('City = %s'%(string_labels[K]),fontsize=20)
i=0
for q in range(3):
    K = np.random.choice(len(images))
    #print(K)
    skyline_signal_x,skyline_signal_y =signal_from_profile(K)
    plt.subplot(4,4,i+1)
    plt.title('Original Image')
    plt.imshow(images[K])
    plt.subplot(4,4,i+3)
    plt.title('Masked Image')
    plt.imshow(zero_images_plot[K])
    plt.subplot(4,4,i+2)
    plt.title('Profiled Skyline')
    plt.imshow(get_sky_region_gradient(images[K]))
    plt.subplot(4,4,i+4)
    plt.title('Extracted Signal')
    plt.plot(skyline_signal_x,skyline_signal_y)
    plt.tight_layout()
    i=i+4

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

• Оно объясняет, как использовать фильтр Лапласа для определения границ в режиме неглубокого обучения.

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

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

Как видно, в городах A и B разные профили, и, благодаря извлеченному сигналу, мы можем углубить это исследование:

• Извлечение среднего значения, медианы и стандартного отклонения линий горизонта

• Использовать глубокое обучение для классификации линии горизонта города

• Провести статистическое исследование зависимости горизонта от времени (как меняется горизонт со временем?)

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

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

Недавно в OTUS прошел открытый урок, на котором мы поговорили о том, что такое Машинное Обучение, какие задачи решают Data Scientists и чем ML отличается от классического программирования; почему специалисты по ML сегодня так востребованы и где применяют современные методы ML. Если вам интересна эта тема, предлагаем посмотреть запись этого урока.

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