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

TL; DR;

U-NET - это искусственная нейронная сеть на основе ConvNets, способная генерировать визуальную информацию. Благодаря вычислительному обучению и четко определенной формуле оптимизации удалось получить приемлемые результаты (~ 0,9 по Dice Score) для сегментации костей и печени на КТ-сканировании.

Вступление

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

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

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

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

Фон Теория

Искусственные нейронные сети появились в 1940-х годах, когда исследования Маккалока и Питтса привели к модели нейрона, разработанной на основе перцептрона Розенблатта.

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

Также в 1989 году ЛеКун опубликовал свою работу о сверточных нейронных сетях. Проще говоря, CNN использует свертку для вычисления весов нейронных сетей, таким образом уменьшая размерность по слоям, сохраняя при этом важные пространственные отношения.

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

Проще говоря, свертка - это математическая операция, которая принимает две функции f и g и создает третью функцию h, которая описывает, как фигура одна функция влияет на другую (ref: Wikipedia).

Дискретная свертка достигается по следующей формуле.

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

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

U-Nets появился в статье 2015 года от Ronneberger et al. и в 2016 в рамках работы Кристи и др. для автоматической сегментации печени на изображениях компьютерной томографии. Отличная идея U-Net заключается в том, что он может получать изображение на входе и создавать другое изображение на выходе, что довольно удобно для создания изображений сегментации.

Идея U-Net заключается в том, что посредством обучения сеть сможет синтезировать релевантную информацию в первой половине сети, минимизируя функцию затрат, связанную с желаемой операцией, а во второй половине она сможет построить изображение.

Для тех, кто знаком с ИНС применительно к машинному обучению, эта концепция кажется очень знакомой с концепцией автоэнкодеров, как объясняется в официальном документе 2012 года от Пьера Бальди. Действительно, идея синтеза релевантной информации и преобразования размерности N➡N очень похожа, но можно утверждать, что обе идеи преследуют разные цели. Тем не менее, если вы знакомы с автоэнкодерами, вам будет легче следовать документам Роннебергера и Христа.

Подробнее о ConvNets

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

Для борьбы с этими сценариями ConvNets могут применять дополнительные уровни, такие как объединение и выпадение.

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

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

Более подробные сведения об обоих уровнях можно найти в объяснении CS231n о CNN.

Задача

Целью этой публикации является применение U-Net от Роннебергера с использованием Tensorflow с Keras на компьютерной томографии для сегментации печени и костей.

Материалы

Для этого выбран набор данных компьютерной томографии 3D IRCAD. Этот набор данных содержит анонимную компьютерную томографию 20 пациентов (10 мужчин и 10 женщин). Изображения предоставлены авторами в формате DICOM и VTK в формате 512x512 пикселей. Одним из недостатков этого набора данных для нашей цели является то, что 75% пациентов имеют опухоль печени, что не является основной целью сегментации для этого сообщения и может представлять собой шум в данных, которые мы пытаемся выделить.

Для кодирования мы в основном будем использовать Python, TensorFlow, Keras и Jupyter для применения U-Net.

Поскольку изображения из выбранного набора данных имеют формат DICOM, мы будем использовать библиотеку python для извлечения этих данных, и для удобства мы будем использовать SciPy ndimage для обработки, изменения размера и PIL для сохранения изображения.

Тензорфлоу и Керас

Tensorflow - фреймворк машинного обучения (см. Tensorflow.io). TF имеет несколько структур для разработки сложных нейронных сетей на основе различных моделей, таких как CNN и RNN. Он обеспечивает реализацию алгоритмов оптимизации, определенных типов NN (например, LSTM), перекрестную проверку, пакетное обучение. Он также предоставляет Tensorboard, который представляет собой инструмент для визуализации сетевого графика и отслеживания процессов обучения и тестирования.

Tensorflow использует подход, основанный на графах, когда пользователь строит граф для модели, в которой каждый узел представляет операцию, ввод или вывод. Затем этот график составляется для повышения производительности. Граф присоединяется к сеансу, который может выполнять свою операцию на CPU, GPU или других узлах сетевой обработки. Как выбор аппаратного устройства, так и сетевая кластеризация легко выполняются, следуя некоторым ручным инструкциям.

Keras - это высокоуровневый API нейронной сети, поддерживающий Tensorflow и Theano (см. Keras.io). Keras отвечает за построение общих слоев на основе нейронных сетей, то есть: плотных, объединяемых, сглаженных и т. Д. Такая функция значительно ускоряет разработку модели.

Формат DICOM

Формат DICOM (Digital Imaging and Communications in Medicine) обычно используется для передачи данных медицинских изображений. Помимо самого изображения, формат содержит идентификационные и информационные теги. Пиксельные данные внутри файла DICOM могут быть сжаты с использованием форматов JPEG, JPEG2000, JPEG без потерь, RLE и LZW (см. Формат Wikipedia DICOM).

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

Выбранная библиотека также называется dicom, и ее установка проста:

pip install dicom

Чтение из файлов DICOM достигается за счет:

import dicom
def getDicomFile(filename):
    return dicom.read_file(filename)

Но, как упоминалось ранее, формат DICOM содержит не только данные изображения, но и другую информацию.

Чтобы получить данные изображения из файла DICOM, можно получить доступ к элементу pixel_data

import dicom
def getPixelDataFromDicom(filename):
    return dicom.read_file(filename).pixel_array

Набор данных 3D IRCAD также содержит ручную сегментацию печени, костей, опухолей и т. Д., Выполненную медицинскими специалистами для всех изображений 20 пациентов.

Данные DICOM из 3D IRCAD анонимны, но для справки сохраняется мало информации, например пол пациента и примерный возраст, как показано на официальной веб-странице и на рисунке ниже. Для этой конкретной работы мы не будем использовать эту информацию, но она может быть актуальной для будущих работ в том же наборе данных.

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

U-Net, представленная в статье Роннебергера, рассматривает изображение 572x572 пикселей как входное и выходное. На слое с самым низким разрешением он имеет 32x32 пикселя до сведения.

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

Пустые пары изображений (CT + Mask) отбрасывались, так как не помогали в тренировочном процессе.

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

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

В 3D IRCAD данные пациента распределены по двум папкам. В папке PATIENT_DICOM находятся реальные данные компьютерной томографии, а в MASKS_DICOM - соответствующие файлы сегментации.

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

Структура папок зависит от пациента, но обычно аналогична структуре папок пациента №1:

.
├── 3Dircadb1.1
│   ├── LABELLED_DICOM
│   ├── MASKS_DICOM
│   │   ├── artery
│   │   ├── bone
│   │   ├── leftkidney
│   │   ├── leftlung
│   │   ├── liver
│   │   ├── liverkyst
│   │   ├── livertumor01
│   │   ├── livertumor02
│   │   ├── livertumor03
│   │   ├── livertumor04
│   │   ├── livertumor05
│   │   ├── livertumor06
│   │   ├── livertumor07
│   │   ├── portalvein
│   │   ├── rightkidney
│   │   ├── rightlung
│   │   ├── skin
│   │   ├── spleen
│   │   └── venoussystem
│   ├── MESHES_VTK
│   ├── PATIENT_DICOM

Перед обучением нам сначала нужно (1) преобразовать изображения в более простой для работы формат и (2) создать входные и целевые наборы.

Преобразование из формата DICOM выполняется путем получения pixel_array и его сохранения в более удобном формате. В данном случае я выбрал формат PNG.

Предварительная обработка изображения выглядит следующим образом:

import dicom
from scipy.ndimage.interpolation import zoom
from scipy.misc import imsave
def add_gaussian_noise(inp, expected_noise_ratio=0.05):
        image = inp.copy()
        if len(image.shape) == 2:
            row,col= image.shape
            ch = 1
        else:
            row,col,ch= image.shape
        mean = 0.
        var = 0.1
        sigma = var**0.5
        gauss = np.random.normal(mean,sigma,(row,col)) * expected_noise_ratio
        gauss = gauss.reshape(row,col)
        noisy = image + gauss
        return noisy
def normalize(img):
        arr = img.copy().astype(np.float)
        M = np.float(np.max(img))
        if M != 0:
            arr *= 1./M
        return arr
def preprocess(filename, resize_ratio=0.25):
    img = dicom.read_file(filename).pixel_array
    img = normalize(zoom(img, resize_ratio))
    img = add_gaussian_noise(img)
    return img
### filelist contains all *.dcm files from 3D IRCAD folder
for dicom_file in filelist:
    pp_image = preprocess(dicom_file)    
    imsave(dicom_file.replace("dcm","png"), pp_image, "png")

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

U-Net с TF + Keras

Окончательная архитектура этого эксперимента состояла из 10 сверточных слоев и была получена от 0,87 до 0,93 по Dice Score (в зависимости от некоторых гиперпараметров, которые мы рассмотрим позже). Christ et al. В бумаге 2016 года использовалось 19 слоев, и результат составил от 0,93 до 0,94 по Dice Score.

Окончательная версия, как показано на изображении выше, реализованная для этого поста, не является ванильной версией из статьи Роннебергера и др. 2015 года, ни Христом и др. 2016 года. Основные отличия заключаются во включении отсеиваемых слоев после каждой фазы объединения и количестве слоев. Стратегия отсева была выбрана из-за ограниченного количества информации из 3D IRCAD. Без отсева слоев тренировка потребовала большего количества взаимодействий, чтобы сойтись (если конвергенция) из-за переобучения. Количество слоев было уменьшено для ускорения процесса разработки, а не по техническим причинам. Но, по крайней мере, для указанного выше набора данных и объектов результаты были вполне удовлетворительными.

Как упоминалось выше, U-Net для этого проекта состояла из 10 сверточных слоев. Иными словами, из оригинальных произведений был добавлен дропаут.

Создание такой структуры с помощью Keras выглядит так:

def get_model(optimizer, loss_metric, metrics, lr=1e-3):
    inputs = Input((sample_width, sample_height, 1))
    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(inputs)
    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv1)
    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)
    drop1 = Dropout(0.5)(pool1)

    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(drop1)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)
    drop2 = Dropout(0.5)(pool2)

    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(drop2)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)
    drop3 = Dropout(0.3)(pool3)

    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(drop3)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    pool4 = MaxPooling2D(pool_size=(2, 2))(conv4)
    drop4 = Dropout(0.3)(pool4)

    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(drop4)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv5)

    up6 = concatenate([Conv2DTranspose(256, (2, 2), strides=(2, 2), padding='same')(conv5), conv4], axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)

    up7 = concatenate([Conv2DTranspose(128, (2, 2), strides=(2, 2), padding='same')(conv6), conv3], axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)

    up8 = concatenate([Conv2DTranspose(64, (2, 2), strides=(2, 2), padding='same')(conv7), conv2], axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)

    up9 = concatenate([Conv2DTranspose(32, (2, 2), strides=(2, 2), padding='same')(conv8), conv1], axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)

    conv10 = Conv2D(1, (1, 1), activation='sigmoid')(conv9)

    model = Model(inputs=[inputs], outputs=[conv10])

    model.compile(optimizer=optimizer(lr=lr), loss=loss_metric, metrics=metrics)
return model

Наилучшие результаты были достигнуты с оптимизатором Adam. Adam Optimizer - это расширение для стохастического градиентного спуска, в котором применяются идеи RMSProp и AdaGrad для адаптации скорости обучения во время обучения. Ссылка на статью для Адама здесь.

Коэффициент Соренсена – Дайса используется для сравнения двух разных статистических выборок, но он также используется в областях обработки изображений для сопоставления сходства между двоичными изображениями (см. Wikipedia).

Функция потерь задается отрицательным значением коэффициента Дайса.

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

#smooth = 1.
# Dice Coefficient to work with Tensorflow
def dice_coef(y_true, y_pred):
    y_true_f = K.flatten(y_true)
    y_pred_f = K.flatten(y_pred)
    intersection = K.sum(y_true_f * y_pred_f)
    return (2. * intersection + smooth) / (K.sum(y_true_f) + K.sum(y_pred_f) + smooth)

def dice_coef_loss(y_true, y_pred):
    return -dice_coef(y_true, y_pred)

# Dice Coefficient to work outside Tensorflow
def dice_coef_2(y_true, y_pred):
    side = len(y_true[0])
    y_true_f = y_true.reshape(side*side)
    y_pred_f = y_pred.reshape(side*side)
    intersection = sum(y_true_f * y_pred_f)
    return (2. * intersection + smooth) / (sum(y_true_f) + sum(y_pred_f) + smooth)

Учитывая модель, представленную выше, и функцию потерь, у нас есть модель для обучения:

model = get_model(optimizer=Adam, loss_metric=dice_coef_loss, metrics=[dice_coef], lr=1e-3)

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

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

Другие результаты

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

Дальнейшая работа

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

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

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

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

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

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

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