Сегментация изображений с помощью Python

Оглавление

Введение в сегментацию изображений

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

Пример кода для этой статьи можно найти в репозитории Github. Мы привели советы по использованию кода во всех статьях.

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

Специалисты по обработке данных и медицинские исследователи могут использовать этот подход в качестве шаблона для любых сложных наборов данных, основанных на изображениях (например, астрономических данных), или даже для больших наборов данных, не связанных с изображениями. В конце концов, изображения - это матрицы значений, и нам повезло, что у нас есть отсортированный экспертами набор данных, который можно использовать в качестве "базовой истины". В этом процессе мы раскроем и опишем несколько инструментов, доступных через пакеты обработки изображений и научные пакеты Python (opencv, scikit-image и scikit-learn). Мы также будем активно использовать библиотеку numpy для обеспечения последовательного хранения значений в памяти.

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

Approach

Очистка

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

  • размытые или расфокусированные области
  • несбалансированность переднего и заднего планов (исправляется модификацией гистограммы)

Сегментация

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

Некоторые примеры:

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

Валидация

Мы начинаем с набора истинных данных, которые уже были сегментированы вручную. Для количественной оценки работы алгоритма сегментации мы сравниваем "истинные" данные с предсказанной бинарной сегментацией, показывая точность наряду с более эффективными метриками. Точность может быть аномально высокой, несмотря на низкое количество истинно положительных (TP) или ложноотрицательных (FN) результатов. В таких случаях для бинарной классификации лучше использовать метрики F1 Score и MCC, которые

подробно рассмотрят достоинства и недостатки этих метрик.

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

Загрузка и визуализация данных

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

Module Reason
numpy Histogram calculation, array math, and equality testing
matplotlib Graph plotting and Image visualization
scipy Image reading and median filter
cv2 (opencv-python) Alpha compositing to combine two images
skimage (scikit-image) Image thresholding
sklearn (scikit-learn) Binary classifier confusion matrix
nose Testing

Отображение графиков Sidebar: Если вы выполняете код примера в секциях из командной строки или испытываете проблемы с бэкендом matplotlib, отключите интерактивный режим, удалив вызов plt.ion(), и вместо этого вызывайте plt.show() в конце каждой секции, некомментируя предложенные вызовы в коде примера. В качестве бэкенда для вывода изображений будет использоваться либо 'Agg', либо 'TkAgg'. Графики будут отображаться в том виде, в котором они представлены в статье.

Импорт модулей

import cv2
import matplotlib.pyplot as plt
import numpy as np
import scipy.misc
import scipy.ndimage
import skimage.filters
import sklearn.metrics

# Turn on interactive mode. Turn off with plt.ioff()
plt.ion()

В этом разделе мы загружаем и визуализируем данные. Данные представляют собой изображение ткани мозга мыши, окрашенной индийскими чернилами, полученное методом сканирующей микроскопии Knife-Edge Scanning Microscopy (KESM). Это изображение размером 512 x 512 является подмножеством, называемым tile. Полный набор данных имеет размер 17480 x 8026 пикселей, глубину 799 срезов и размер 10 Гб. Таким образом, мы напишем алгоритмы для обработки тайла размером 512 x 512, что составляет всего 150 Кб.

Отдельные тайлы могут быть отображены для работы на многопроцессорной/многопоточной (т.е. распределенной) инфраструктуре, а затем сшиты вместе для получения полного сегментированного изображения. Конкретный метод сшивания здесь не демонстрируется. Вкратце, сшивание заключается в индексации полной матрицы и последующем соединении плиток в соответствии с этим индексом. Для объединения числовых значений можно использовать map-reduce. Map-Reduce позволяет получить такие метрики, как сумма всех оценок F1 по всем плиткам, которые затем можно усреднить. Просто добавьте результаты в список, а затем выполните собственную статистическую сводку.

Темные круглые/эллиптические диски слева - сосуды, остальное - ткань. Таким образом, в этом наборе данных мы имеем два класса:

  • передний план (сосуды) - обозначен как 255
  • задний план (ткани) - обозначен как 0

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

Загрузка и визуализация изображений на рисунке выше

grayscale = scipy.misc.imread('grayscale.png')
grayscale = 255 - grayscale
groundtruth = scipy.misc.imread('groundtruth.png')
plt.subplot(1, 3, 1)
plt.imshow(255 - grayscale, cmap='gray')
plt.title('grayscale')
plt.axis('off')
plt.subplot(1, 3, 2)
plt.imshow(grayscale, cmap='gray')
plt.title('inverted grayscale')
plt.axis('off')
plt.subplot(1, 3, 3)
plt.imshow(groundtruth, cmap='gray')
plt.title('groundtruth binary')
plt.axis('off')

Предварительная обработка

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

Медианный фильтр с размером ядра 3

median_filtered = scipy.ndimage.median_filter(grayscale, size=3)
plt.imshow(median_filtered, cmap='gray')
plt.axis('off')
plt.title('median filtered image')

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

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

Визуализация гистограммы интенсивностей пикселей

counts, vals = np.histogram(grayscale, bins=range(2 ** 8))
plt.plot(range(0, (2 ** 8) - 1), counts)
plt.title('Grayscale image histogram')
plt.xlabel('Pixel intensity')
plt.ylabel('Count')

Сегментация

После удаления шумов можно применить модуль фильтров skimage и попробовать все пороговые значения, чтобы выяснить, какие методы порогового преобразования работают лучше. Иногда гистограмма интенсивности пикселей на изображении не является бимодальной. В этом случае может быть использован другой метод пороговой обработки, например, адаптивный метод пороговой обработки, основанный на локальных интенсивностях пикселей в пределах формы ядра. Полезно посмотреть, какие результаты дают различные методы порогового вычисления, и для этого удобно использовать skimage.filters.thresholding.try_all_threshold().

Попробовать все пороговые методы

result = skimage.filters.thresholding.try_all_threshold(median_filtered)

Простейший подход к пороговой обработке использует ручную установку порога для изображения. С другой стороны, при использовании автоматизированного порогового метода на изображении его численное значение вычисляется лучше, чем человеческий глаз, и может быть легко воспроизведено. Для нашего изображения в данном примере, похоже, хорошо работают Оцу, Йен и метод "Треугольник". Остальные результаты для данного случая заметно хуже.

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

Пороговая обработка и визуализация Оцу

threshold = skimage.filters.threshold_otsu(median_filtered)
print('Threshold value is {}'.format(threshold))
predicted = np.uint8(median_filtered > threshold) * 255
plt.imshow(predicted, cmap='gray')
plt.axis('off')
plt.title('otsu predicted binary image')

Если приведенные выше простые методы не подходят для бинарной сегментации изображения, то для сегментации изображений можно использовать UNet, ResNet с FCN или различные другие контролируемые методы глубокого обучения. Для удаления мелких объектов из-за шума на переднем плане сегментированного изображения можно также попробовать skimage.morphology.remove_objects().

Валидация

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

Матрица смешения

Для получения элементов матрицы путаницы используем sklearn.metrics.confusion_matrix(), как показано ниже. Функция матрицы путаницы Scikit-learn возвращает 4 элемента матрицы путаницы, учитывая, что на вход подается список элементов с бинарными элементами. Для крайних случаев, когда все является одним двоичным значением(0) или другим(1), sklearn возвращает только один элемент. Мы оборачиваем функцию матрицы путаницы sklearn и пишем свою собственную с учетом этих крайних случаев следующим образом:

get_confusion_matrix_elements()

def get_confusion_matrix_elements(groundtruth_list, predicted_list):
    """returns confusion matrix elements i.e TN, FP, FN, TP as floats
	See example code for helper function definitions
    """
    _assert_valid_lists(groundtruth_list, predicted_list)

    if _all_class_1_predicted_as_class_1(groundtruth_list, predicted_list) is True:
        tn, fp, fn, tp = 0, 0, 0, np.float64(len(groundtruth_list))

    elif _all_class_0_predicted_as_class_0(groundtruth_list, predicted_list) is True:
        tn, fp, fn, tp = np.float64(len(groundtruth_list)), 0, 0, 0

    else:
        tn, fp, fn, tp = sklearn.metrics.confusion_matrix(groundtruth_list, predicted_list).ravel()
        tn, fp, fn, tp = np.float64(tn), np.float64(fp), np.float64(fn), np.float64(tp)

    return tn, fp, fn, tp

Точность

Точность является общей метрикой валидации в случае бинарной классификации. Она рассчитывается как

accuracy

где TP = True Positive, TN = True Negative, FP = False Positive, FN = False Negative

get_accuracy()

def get_accuracy(groundtruth_list, predicted_list):

    tn, fp, fn, tp = get_confusion_matrix_elements(groundtruth_list, predicted_list)
    
    total = tp + fp + fn + tn
    accuracy = (tp + tn) / total
    
    return accuracy

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

F1 score

Оценка F1 варьируется от 0 до 1 и рассчитывается следующим образом:

F1 Score

при этом 0 - худший, а 1 - лучший прогноз. Теперь займемся вычислением оценки F1 с учетом крайних случаев.

get_f1_score()

def get_f1_score(groundtruth_list, predicted_list):
    """Return f1 score covering edge cases"""

    tn, fp, fn, tp = get_confusion_matrix_elements(groundtruth_list, predicted_list)
    
    if _all_class_0_predicted_as_class_0(groundtruth_list, predicted_list) is True:
        f1_score = 1
    elif _all_class_1_predicted_as_class_1(groundtruth_list, predicted_list) is True:
        f1_score = 1
    else:
        f1_score = (2 * tp) / ((2 * tp) + fp + fn)

    return f1_score

Хорошим показателем F1 считается значение, превышающее 0,8, что говорит о том, что прогнозирование идет успешно.

MCC

MCC обозначает коэффициент корреляции Мэтьюса и рассчитывается следующим образом:

mcc equation

Он лежит между -1 и +1. -1 - абсолютно противоположная корреляция между истиной и предсказанием, 0 - случайный результат, когда некоторые предсказания совпадают, а +1 - когда абсолютно все совпадает между истиной и предсказанием, что приводит к положительной корреляции. Следовательно, нам нужны более эффективные метрики валидации, такие как MCC.

При вычислении MCC числитель состоит только из четырех внутренних ячеек (кросс-произведение элементов), а знаменатель - из четырех внешних ячеек (точечное произведение элементов) матрицы путаницы. В случае, когда знаменатель равен 0, MCC сможет заметить, что ваш классификатор идет в неправильном направлении, и сообщит вам об этом, установив неопределенное значение (т.е. numpy.nan). Но для того чтобы получить корректные значения и при необходимости усреднить MCC по разным изображениям, мы установили MCC равным -1 - наихудшему из возможных значений в этом диапазоне. Другие крайние случаи включают в себя все элементы, правильно определенные как передний и задний план, для которых MCC и F1 score установлены в 1. В противном случае MCC устанавливается равным -1, а F1 - 0.

Для получения более подробной информации о MCC и крайних случаях, это хорошая статья. Для более детального понимания того, почему MCC лучше, чем точность или F1 score, хорошо подходит Википедия здесь.

get_mcc()

def get_mcc(groundtruth_list, predicted_list):
    """Return mcc covering edge cases"""   

    tn, fp, fn, tp = get_confusion_matrix_elements(groundtruth_list, predicted_list)
    
    if _all_class_0_predicted_as_class_0(groundtruth_list, predicted_list) is True:
        mcc = 1
    elif _all_class_1_predicted_as_class_1(groundtruth_list, predicted_list) is True:
        mcc = 1
    elif _all_class_1_predicted_as_class_0(groundtruth_list, predicted_list) is True:
        mcc = -1
    elif _all_class_0_predicted_as_class_1(groundtruth_list, predicted_list) is True :
        mcc = -1

    elif _mcc_denominator_zero(tn, fp, fn, tp) is True:
        mcc = -1

    # Finally calculate MCC
    else:
        mcc = ((tp * tn) - (fp * fn)) / (
            np.sqrt((tp + fp) * (tp + fn) * (tn + fp) * (tn + fn)))
    
    return mcc

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

>>> validation_metrics = get_validation_metrics(groundtruth, predicted)

{'mcc': 0.8533910225863214, 'f1_score': 0.8493358633776091, 'tp': 5595.0, 'fn': 1863.0, 'fp': 122.0, 'accuracy': 0.9924278259277344, 'tn': 254564.0}

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

F1 score составляет 0,84. Таким образом, в данном случае нам, скорее всего, не нужен более сложный пороговый алгоритм для бинарной сегментации. Если бы все изображения в стеке имели схожее распределение гистограмм и шумов, то можно было бы использовать Otsu и получить удовлетворительные результаты предсказания.

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

Теперь визуализируем и посмотрим, как распределены по изображению элементы матрицы путаницы TP, FP, FN, TN. Это показывает, где порог улавливает передний план (сосуды), когда их нет (FP), а где истинные сосуды не обнаружены (FN), и наоборот.

Визуализация валидации

Для визуализации элементов матрицы путаницы выясним, куда именно на изображении попадают элементы матрицы путаницы. Например, мы находим массив TP (т.е. пикселей, правильно распознанных как передний план) путем нахождения логического "и" между истинным и предсказанным массивами. Аналогично, для нахождения массивов FP, FN, TN мы используем логические булевы операции, которые обычно называются Bit blit

get_confusion_matrix_intersection_mats()

def get_confusion_matrix_intersection_mats(groundtruth, predicted):
    """ Returns dict of 4 boolean numpy arrays with True at TP, FP, FN, TN
    """

    confusion_matrix_arrs = {}

    groundtruth_inverse = np.logical_not(groundtruth)
    predicted_inverse = np.logical_not(predicted)

    confusion_matrix_arrs['tp'] = np.logical_and(groundtruth, predicted)
    confusion_matrix_arrs['tn'] = np.logical_and(groundtruth_inverse, predicted_inverse)
    confusion_matrix_arrs['fp'] = np.logical_and(groundtruth_inverse, predicted)
    confusion_matrix_arrs['fn'] = np.logical_and(groundtruth, predicted_inverse)

    return confusion_matrix_arrs

Затем мы можем сопоставить пиксели в каждом из этих массивов разным цветам. Для приведенного ниже рисунка мы отобразили TP, FP, FN, TN в пространство CMYK (Cyan, Magenta, Yellow, Black). Аналогичным образом их можно сопоставить с цветами (Green, Red, Red, Green). Тогда мы получим изображение, в котором все, что выделено красным цветом, означает неверные прогнозы. Пространство CMYK позволяет различать TP, TN.

get_confusion_matrix_overlaid_mask()

def get_confusion_matrix_overlaid_mask(image, groundtruth, predicted, alpha, colors):
    """
    Returns overlay the 'image' with a color mask where TP, FP, FN, TN are
    each a color given by the 'colors' dictionary
    """
    image = cv2.cvtColor(image, cv2.COLOR_GRAY2RGB)
    masks = get_confusion_matrix_intersection_mats(groundtruth, predicted)
    color_mask = np.zeros_like(image)
    for label, mask in masks.items():
        color = colors[label]
        mask_rgb = np.zeros_like(image)
        mask_rgb[mask != 0] = color
        color_mask += mask_rgb
    return cv2.addWeighted(image, alpha, color_mask, 1 - alpha, 0)

alpha = 0.5
confusion_matrix_colors = {
   'tp': (0, 255, 255),  #cyan
   'fp': (255, 0, 255),  #magenta
   'fn': (255, 255, 0),  #yellow
   'tn': (0, 0, 0)     #black
   }
validation_mask = get_confusion_matrix_overlaid_mask(255 - grayscale, groundtruth, predicted, alpha, confusion_matrix_colors)
print('Cyan - TP')
print('Magenta - FP')
print('Yellow - FN')
print('Black - TN')
plt.imshow(validation_mask)
plt.axis('off')
plt.title('confusion matrix overlay mask')

Здесь мы используем opencv для наложения этой цветовой маски на исходное (неинвертированное) полутоновое изображение в виде прозрачного слоя. Это называется Альфа-композицией:

Заключительные замечания

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

Travis CI очень полезен для проверки того, работает ли ваш код на версиях модулей, описанных в требованиях, и проходят ли все тесты при слиянии новых изменений в master. Чистота кода, хорошая документированность и наличие юнит-тестов для всех утверждений - это лучшая практика. Такие привычки ограничивают необходимость поиска ошибок, когда сложный алгоритм строится поверх простых функциональных частей, которые можно было бы протестировать. В целом, документирование и модульное тестирование помогают другим людям быть в курсе ваших намерений относительно функции. Линтинг помогает улучшить читаемость кода, и flake8 - хороший пакет Python для этого.

Вот важные выводы из этой статьи:

  1. Подход сшивки и черепицы для данных, не помещающихся в памяти
  2. Опробование различных методов пороговой обработки
  3. Тонкости метрик валидации
  4. Визуализация валидации
  5. Лучшие практики

Существует множество направлений, по которым вы можете двигаться в своей работе или проектах. Применение одной и той же стратегии к различным наборам данных или автоматизация подхода к выбору валидации - это отличные варианты для начала. Далее, представьте, что вам нужно проанализировать базу данных, содержащую множество таких 10-гигабайтных файлов. Как можно автоматизировать этот процесс? Как проверить и обосновать результаты для человека? Как более качественный анализ может улучшить результаты в реальных сценариях (например, при разработке хирургических процедур и лекарств)? Задавая подобные вопросы, можно продолжать совершенствоваться в области статистики, науки о данных и машинного обучения.

Вернуться на верх