С 12-го июля  по 21-ого августа проходило мероприятие Большая математическая мастерская. Я был руководителем проекта целью которого стала разработка метода оценки пигментного состава зерна ячменя на основе анализа цифровых изображений с помощью алгоритмов компьютерного зрения и машинного обучения. Цвет оболочки зерен злаков – важный признак, характеризующий содержащиеся в ней пигменты и метаболиты. Наличие пигментов в оболочке влияет на различные технологические свойства зерна. Показано, что растения с темной окраской зерна являются более холодо- и засухоустойчивыми, а также обладают повышенной устойчивостью к действию патогенов. Такие свойства окрашенных растений связаны с высоким содержанием антиоксидантов, а также с дополнительной механической прочностью оболочек зерна. Темная окраска колоса ячменя может быть обусловлена синтезом и накоплением двух групп пигментов. Голубая и фиолетовая окраска зерна связана с синтезом антоцианов. Серую и черную окраску ячменя обуславливают пигменты меланины. Данные пигменты могут накапливаться в оболочках зерна самостоятельно, либо совместно, поэтому визуально определить, накопление каких именно пигментов обуславливает темный цвет зерна, затруднительно. Для точного определения наличия/отсутствия пигментов используются химические и генетические методы, которые, однако, дороги и трудозатраты. Поэтому сознание нового метода для быстрой оценки наличия пигментов в зернах является актуальной задачей, решение которой поможет при исследовании механизмов генетического контроля окраски зерен.

Данные

Материал

Для обучения моделей были выбраны семена 41 образца ячменя (Hordeum vulgare) с темной окраской оболочек зерна, а также 38 образцов с неокрашенным зерном. Материал был получен из коллекции ячменя Всероссийского института генетических ресурсов растений имени Н.И. Вавилова (ВИР), коллекции ячменя Института Цитологии и Генетики СО РАН (ИЦиГ), а также был использован материал популяции Oregon Wolfe Barleys (OWB). Для отложенной выборки был отдельно выбран 21 образец ячменя из коллекции ВИР с различными комбинациями пигментов в зерне.

Протокол получения изображения зерен

Для фенотипирования образцов по интенсивности окраски зерна был разработан протокол получения цифровых двумерных изображений. На белый фон помещается чашка Петри с зернами и цветовая палитра ColorChecker Mini Classic target (https://xritephoto.com/camera), что позволяет произвести цветокоррекцию изображения. Такая коррекция позволяет устранить искажения цвета на изображении, возникающие из-за разных условий освещения. Другим преимуществом использования цветовой шкалы является ее стандартный размер 63 x 108 мм, что позволяет оценить масштаб изображения. Цветные изображения зерен получаются с помощью цифровой фотокамеры Canon EOS 600D, объектив Canon EF 100mm f/2.8 Macro USM в формате jpg, разрешение 18Mpx. На рисунке 1 пример изображения, полученного в результате выполнения протокола. В среднем при съемке в чаше Петри находилось 30 зерен принадлежащих одному образцу. Для каждого образца было получено 9 изображений: три случайных выборки (повторности) зерен, для каждой из которых зерна перемешивались по 3 раза.

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

Разметка данных

Рисунок 2. Пример образца с нанесенной разметкой зерен и границ чашки Петри.
Рисунок 2. Пример образца с нанесенной разметкой зерен и границ чашки Петри.

Для разработки моделей предсказания пигментного состава зерна по цифровому изображению пигментный состав оболочек зерна был исследован при помощи качественных химических реакций на антоцианы и меланины в лаборатории. Все образцы были генотипированы с использованием ПЦР-маркеров. Для 212 изображений из 59 образцов была выполнена ручная разметка зерен и границ чаши Петри. Для этого мы использовали программу LabelMe. Пример фрагмента размеченного изображения показан на рисунке 2.

Схема валидации и стратификация данных

Рисунок 3. Распределение количества образцов с пигментами и без пигментов
Рисунок 3. Распределение количества образцов с пигментами и без пигментов

Мы разбили наши изображения на 3 подвыборки: train (60% данных: 423 изображения, 47 образцов), для обучения модели; valid (20% данных: 144 изображения, 16 образцов), для выбора лучшей модели в процессе обучения; test (20% данных: 144 изображений, 16 образцов)  для оценки точности выбранной модели. Для сбалансированности распределения multilable классов в каждой подвыборке мы использовали алгоритм итеративной стратификации (метод iterative_train_test_split() модуль model_selection библиотеки skmultilearn). Для финальной оценки точности мы использовали отложенную (hold out) выборку из 29 образцов включающею 261 изображений. Эти образцы были получены уже после процесса обучения наших моделей, и баланс меток классов для них не соответствовал тренировочной выборке. Для каждого образца из нашего набора данных имеется по 9 изображений. Была использована стратификация по образцам: все изображения одного образца всегда находятся в одной подвыборке.

Оценка качества моделей классификации

Для каждого изображения модель предсказывает 2 бинарных числа, каждое из которых характеризует наличие или отсутствие той или иной метки (антоциан и меланин). Чтобы оценить точность метода на выборке из тестовых изображений мы сравнивали для каждого изображения предсказанный набор таких чисел и истинный набор, так что, если тестовая выборка содержала M изображений, мы проводили 2M таких бинарных сравнений и на их основе подсчитывали величины true positive (TP), true negative (TN), общее количество positive (1, P) и negative (0, N) значений. На основании этих величин вычислялось значение accuracy (ACC) согласно формуле:

ACC = (TP + TN) / (P + N)

Модель предсказания масок зерен

Для получения масок зерен мы обучили сегментацию. Для сегментации взяли U-net с энкодером ResNet-18. Обучались на 261 изображениях, для которых маски зерен были размечены вручную. Разбили наши изображения на 3 подвыборки: train, valid и test случайным образом в соотношении 60%, 20% и 20%, соответственно. При обучении использовали loss-функцию Dice, модель обучали 10 эпох, размер батча 4, оптимизатор Аdam. Точность на test - 0.985 в метрике IoU. Используя данную модель мы получили маски зерен для всех наших изображений.

Извлечение цветовых дескрипторов из масок зерен

Средние значения и дисперсия цвета

Для описания цветовых характеристик зерен мы использовали представление цвета в виде четырех моделей: RGB, HSV, Lab, YCrCb. Каждая из них представляет цвет в виде трех компонент; значения компонент одного пространства могут быть получены путем преобразований значений компонент другого. Диапазон значений цветовых каналов от 0 до 255, кроме канала Hue в цветовой модели HSV, диапазон значений которого от 0 до 180. Цветовые дескрипторы вычислялись независимо для каждой из этих компонент.

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

Глобальная цветовая гистограмма (GCH)

Следующий тип цветовых характеристик – глобальная цветовая гистограмма. Она оценивается путем разбиения трёхмерного цветового пространства на трехмерные бины, размер которых одинаков и определяется числом бинов на компоненту, d. В итоге получается общее количество d3 бинов. Каждый бин определяется индексами (i,j,k =1,…d), соответствующих трем компонентам и интервалами интенсивности по каждой компоненте в пределах ((i-1)∙xdi, ixdi; (j-1)∙xdj, jxdj; (k-1)∙xdk, kxdk), где xdc это размер бина по компоненте c. Мы использовали два типа гистограмм с d=4 (всего 64 дескриптора) и d=8 (всего 256 дескрипторов).

Дескриптор доминантных цветов

Дескриптор доминантных цветов обеспечивает компактное описание репрезентативных цветов на изображении или области изображения. Дескриптор определяется как: F={{ci,pi,vi},s}, (i=1,…,N), где N это число доминантных цветов. Каждый доминантный цвет ci, это вектор соответствующих компонентов цветового пространства (например, 3-D вектор пространства RGB). Процент pi, нормализованный на значения между 0 и 1, это процент пикселей на изображении или области изображения соответствующей цвету ci, и . Опциональная дисперсия цвета vi, описывает вариабельность цвета пикселей в кластере вокруг соответствующего репрезентативного (доминантного) цвета. Пространственная когерентность s – это число, представляющее общую пространственную однородность доминирующих цветов на изображении.

Фильтрация признаков

Из каждого изображения было извлечено 2380 признаков (цветовых дескрипторов). Для  фильтрации малозначимых признаков были удалены признаки:

  1. значения которых одинаковы для всех изображений

  2. значения которых на менее 20% изображений не превышают значения 0,01

  3. признаки с корреляцией Спирмена более 0,97

После фильтрации осталось 345 дескриптора вместо 2380.

Анализ главных компонент

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

Рисунок 4. Диаграмма PCA. Первая компонента объясняет  35,2% дисперсии; вторая 9,9% дисперсии.
Рисунок 4. Диаграмма PCA. Первая компонента объясняет 35,2% дисперсии; вторая 9,9% дисперсии.

На представленной на рисунке 4 диаграмме PCA желтые точки соответствуют изображениям зерен без пигментов. Они хорошо отделимы от изображений в котором присутствуют пигменты. Изображения на которых представлены образцы зерен содержащих один из двух пигментов распадаются на два кластера: образцы содержащие антоцианы (синие точки) образуют кластер в левом нижнем углу диаграммы; образцы содержащие меланины (черные точки) падают в кластер растянутый вдоль второй главной компоненты, расположенный преимущество в центре диаграммы. В этот же кластер попадают образцы содержащие оба пигмента (красные точки).

Интерпретация метода главных компонент позволила сделать вывод о том, что в первую компоненту основной вклад вносит канал “Cr” среднего цвета и дескриптора доминантного цвета в цветовой модели YCrCb. Канал “Cr” это - красная цветоразностная компонента, варьирующаяся от сине-зеленой области до красно-пурпурной.

Методы классификации пигментного состава

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

Random forest на фильтрованных признаках

Рисунок 5. Random forest на фильтрованных признаках
Рисунок 5. Random forest на фильтрованных признаках

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

Классификация пигментного состава на основе CNN

Рисунок 6. Классификация пигментного состава на основе CNN
Рисунок 6. Классификация пигментного состава на основе CNN

На исходных изображениях при помощи модели сегментации мы вырезали фрагмент изображения, который соответствует чашке Петри и подали на вход сети ResNet18. В процессе обучения были применены различные техники аугментации: повороты, смещения, зеркальные отражения, изменения яркости и чёткости изображения и т.д. Модель имеет слои batch-нормализации и dropout=0.5 в полносвязном слое. В качестве функции потерь использовалась бинарная кросс-энтропия комбинированная с функцией активации Sigmoid(). В качестве оптимизатора использовался Adam со начальной скоростью обучения lr=5e-5.

Сегментации с головой для классификации

Рисунок 7. Модель на основе сегментации с головой для классификации
Рисунок 7. Модель на основе сегментации с головой для классификации

Идея данного подхода (рисунок 7) добавить модель на основе сегментации с архитектурой U-Net дополнительный выходной слой (иначе говоря голову) для классификации. В качестве энкодера мы использовали архитектуру EfficientNet-B0. В качестве loss-функций мы использовали комбинацию бинарной кросс-энтропии BCEWithLogitsLoss() с уменьшением по среднему значению и DiceLoss() с функцией активации Sigmoid(). Для оптимизации параметров сети в процессе обучения алгоритм Adam c начальной скоростью обучения lr=1e-4. Длительность обучения 10 эпох, размер батча 2. Сеть обучалась на RGB изображениях, трансформированных до размера 640x640 px с сохранением пропорций.

2-канальная сегментация

Рисунок 8. Модель 2-канальной сегментации
Рисунок 8. Модель 2-канальной сегментации

Еще один способ для создания модели классификации на основе модели для сегментации заключался в следующем: пусть наша сеть будет выдавать не стандартную одноканальную черно-белую маску, а двух-канальную, где первый канал сегментирует изображение, если образец содержит меланин, и не сегментирует если не содержит меланина. Второй канал аналогичен первому но уже нужен для второго пигмента - антоциана (рисунок 8). За основу импользовалась модель для сегментации на основе архитектуры U-net с энкодером ResNet34. В качестве loss-функций мы использовали DiceLoss() с функцией активации Sigmoid(). Для выбора лучшей модели использовалась модифицированная IoU метрика. В отличие от стандартной IoU метрики наша метрика при значении IoU меньшем определенного порога возвращает 0. Порог подбирался на валидационной выборке отдельно для каждого пигмента. Для оптимизации параметров сети использовался алгоритм Adam. Количество эпох 10, размер батча 4.

Результаты

Результаты для каждого из рассмотренных подходов представлены в таблице:

Valid

Test

Hold out

Random forest на фильтрованных признаках

0.896

0.903

0.652

Классификация пигментного состава на основе CNN

0.938

0.934

0.817

Сегментации с головой для классификации

0.906

0.962

0.852

2-канальная сегментация

0.917

0.903

0.845

Лучший результат на отложенной выборке демонстрирует модель сегментации с головой для классификации. Можно сделать вывод о лучшей обобщающей способности моделей которые одновременно решают несколько различных задач (multi task learning). Худший результат у случайного леса на цветовых дескрипторах. Результаты на hold out выборке заметно хуже чем на test. Есть несколько предположений почему так:

  1. Баланс меток различных классов в train, valid и test был одинаковый и не совпадал с соотношением, который пришел к нам в hold out. В частности, изображений с зернами без пигментов в hold out было в 1.5 раза меньше чем в обучающей выборке, а для классификации такие зерна оказываются наиболее простым случаем.

  2. На основании цветовых дескрипторов, которые мы извлекали был обучен бинарный классификатор который отличал зерна из hold out от других зерен с accuracy = 1. Это означает, что между данными сериями изображений есть существенные различия которые могут объясниться тем, что в hold out были выбраны зерна из других коллекций или протокол съемки этих изображений был немного другим. Этим можно объяснить наибольшую просадку в качестве классификации у модели случайного леса. Мы намеренно приняли решения не нормализировать изображение на основе цветовой палитры ColorChecker перед извлечением дескрипторов. Это, вероятно, было ошибкой. Нейронные сети лучше справились с задачей за счет использования аугментаций в процессе обучения.

Благодарности

Хотел бы выразить благодарности Математическому центру в Академгородке за организацию мероприятия. Глаголевой Анастасии Юрьевне за предоставленный биологический материал. Куратором проекта Евгению Комышеву, Константину Генаеву и участникам проекта Никите Артеменко, Игорю Бусову, Евгению Заварзину, Анне Ивлевой и Михаилу Кожекину за продуктивную работу над задачами проекта.

Комментарии (6)


  1. Tomasina
    02.09.2021 14:42

    А возможен анализ проросших семян по цвету листьев - на признаки их болезни или неправильного питания?


    1. genaev Автор
      02.09.2021 15:17
      +1

      @Tomasina при наличии хорошего дата сета, думаю можно. Мы делали что-то подобное для грибковых заболеваниях у пшеницы. Там заболевания определялись в том числе и у молодых растений на стадии проростков. Вот статья про это https://www.mdpi.com/2223-7747/10/8/1500


  1. pdkdrp
    02.09.2021 22:10

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

    https://www.foxtechfpv.com/ms600-series-multi-spectral-camera.html

    Или такой:

    https://agrowing.com/products/alpha-7riv-sextuple/

    И все обучение сведется к порогову детектору.


    1. Evgenius2020
      03.09.2021 07:28

      Интересная идея. В нашей работе мы оцениваем наличие 0, 1 или 2 пигментов - изображения при этом будут иметь различные спектры. Выходит, нам нужно перенастраивать спектрометр или использовать сразу несколько? Наша задача не в сегментации (хотя это важный этап). Нам нужно также определить наличие пигментов.
      Вопрос в том, можно ли подобрать подходящие настройки для камер - это - полноценная исследовательская задача.. В целом, если верить результатам, задачу можно решить - тот же RandomForest классификатор по цветовым дескрипторам может решить задачу (если опустить ситуацию с hold out).


    1. genaev Автор
      03.09.2021 07:39

      @pdkdrp да, конечно. Следует ожидать, что использование мутиспектральных данных поможет получить более точную модель по сравнению с моделями на обычных RGB изображениях. Смущает, однако, стоимость такого оборудования. Не каждая лаборатория может позволить себе такую камеру. Есть данные сколько стоит такой объективчик на Sony? У нас был опыт работы с гиперспектральной камерой Specim IQ для анализа растений. Вот статья на эту тему https://doi.org/10.18699/VJ19.587


      1. pdkdrp
        03.09.2021 09:21

        >>Есть данные сколько стоит такой объективчик на Sony?

        Мультиспектральные камеры для растениеводства 2-6K$, причем, если глянуть мою первую ссылку, то похоже заметная часть цены - это софт, система калибровки и т. д.

        Узкополосный фильтр стоит ~100$, т.е. найти гиперспектрометр через знакомых (хабр?), затем, по снятым спектрам подобрать фильтры (n x 100$), прицепить их хоть на пачку RPI HQ или через распечатанный переходник завести на матрицу всевдозеркалки... вроде адекватно по деньгам получается.

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