Изобретаем успех: софт и стартапы
75,72
рейтинг
16 декабря 2015 в 18:11

Разработка → Как за 5233 человеко-часа создать софт для микротомографа



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

А еще сегодня 16 декабря, день рождения Иоганна Радона, австрийского математика, ректора Венского университета, который в 1917 году ввел интегральное преобразование функции многих переменных, родственное преобразованию Фурье, используемое сегодня во всех томографах.

Иоганн Радон был профессором 6 университетов (а в одном из них даже без кафедры), был президентом Австрийского математического общества. В Австрии в честь него назвали «Институт вычислительной и прикладной математики» и медаль.

О том, как проходила разработка софта для томографа и какие задачи решались в процессе — под катом.


Микротомограф





Ученые из Томского государственного университета создали микротомограф. Он позволяет с точностью до микрона узнать о внутренней структуре различных материалов. Например, алмазов.


Рис. Качество алмазов

Томограф может просветить материал с разрешением до микрона. Это в 100 раз тоньше человеческого волоса. После сканирования программа создает 3D-модель, где можно посмотреть не только на внешнюю сторону детали, но и узнать, что у нее внутри.

Устройство томографа
Список основных технических характеристик прибора.

• Количество разрешающих элементов детектора: 2048 x 2018 ячеек при размере одного элемента не более 13,3 x 13,3 мкм.
• Разрешение: 13 мкм.
• Габаритные характеристики: 504 x 992,5 x 1504 мм.
• Масса: 450 кг.
• Поле зрения: 1 мм.
• Рабочий диапазон длин волн: 0,3–2,3 А.
• Точность позиционирования электромеханических модулей движения в системе позиционирования РМТ: ±1 мкм.
• Соответствие требованиям безопасности: ГОСТ 12.1.030-81.
• Защита от рентгеновского излучения: 1–3 мкЗв/ч.


Рис. Внешний вид


Рис. Принцип работы электро-механической части


Рис. Источник рентгеновского излучения

• Напряжение: 20–160 кВ.
• Сила тока: 0–250 мкА.
• Мощность: 10 Вт.
• Диаметр фокального пятна: 1–5 мкм.


Рис. Позиционирование

• Ход ротора: СЭМД — 360°, ЛЭМД — 100 мм.
• Точность позиционирования, не менее: ±0,5 мкм.
• Скорость перемещения ротора: от 0,01 до 20 мм/с.
• Мощность: 0,7 кВт при напряжении 70 В.


Рис. Детектор рентгеновского излучения на базе ПЗС-матрицы

• Чувствительная область ПЗС: 2048 x 2048 пикселей.
• Геометрический размер пикселя: 13 x 13 мкм.
• Геометрический размер чувствительной области ПЗС: 26,6 x27,6 мм.
• Встроенный двухскоростной АЦП: 16 бит, 100 кГц и 16 бит 2МГц.



Рис. Иллюстрация оцифровки и реконструкции

Используемые математические алгоритмы.
  • Обратное преобразование Радона.
  • Марширующие/шагающие кубы.
  • Фильтр Гаусса.
  • Фильтрация/свертка, нормализация проекций.

Реализация и технологии.
  • C++/Qt.
  • Ubuntu, Windows.
  • CUDA.
  • Объемный-воксельный рендеринг модели.
  • Собственный формат хранения изображений с оттенками серого глубиной 12 бит.
  • Разделение реконструкции на клиент-сервер.
  • Сохранение объемных данных в виде «Октодеревьев».

Трудозатраты: 5233 человеко-часа.
Отладка производилась на сырых данных (проекционных снимках), полученных с томографа.

Алгоритмы


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

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

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


Для построения 3D-моделей поверхностей объектов в стандартном формате для просмотра в 3D-редакторах Компас, SolidWork, 3D Max Studio использовался алгоритм «Марширующие (шагающие) кубы». Суть алгоритма в том, что он пробегает скалярное поле, на каждой итерации просматривает 8 соседних позиций (вершины куба, параллельного осям координат) и определяет полигоны, необходимые для представления части изоповерхности, проходящей через данный куб. Далее, на экран выводятся полигоны, образующие заданную изоповерхность.

Под фильтром Гаусса в проекте понимается матричный фильтры обработки изображений с использованием матрицы свертки. Матрица свёртки – это матрица коэффициентов, которая «умножается» на значение пикселей изображения для получения требуемого результата. Фильтр используется для сглаживания воксельных данных и проекций срезов, что, в свою очередь, позволяет улучшить качество генерируемых 3D-моделей (Хабрапост 1, Хабрапост 2).

Реализация и технологии


В ходе работы также был создан ряд специализированных технических решений: библиотека для объемно-воксельного рендеринга модели; запись видео при выполнении операций с моделью; собственный формат хранения изображений с оттенками серого глубиной 12 бит; сохранение объемных данных в виде «Октодеревьев»; алгоритмы полигонизации 3D-модели.

image
Объемный-воксельный рендеринг модели в проекте использовался для просмотра модели с возможностью вращения и масштабирования в реальном времени. Воксель — это трехмерный пиксель. Так же с помощью воксельного рендеринга оператору предоставляются удобные инструменты для определения области просмотра с автоматическим увеличением уровня детализации и возможность расположения секущей плоскости под любым углом в два клика. На основе секущей плоскости в дальнейшем можно получить изображение среза с максимальным разрешением.
imageОктодерево (дерево октантов, восьмеричное дерево, англ. octree) — тип древовидной структуры данных, в которой у каждого внутреннего узла ровно восемь «потомков». Восьмеричные деревья чаще всего используются для разделения трёхмерного пространства, рекурсивно разделяя его на восемь ячеек. В проекте октодерево позволяет выполнять отображение данных в режиме предварительного просмотра, когда нет необходимости в данных, которые не видны пользователю. Так, например, объемный рендеринг получает набор данных для отображения с детализацией в зависимости от выбранной области с помощью октодерева, что обеспечивает высокий FPS, когда видна вся модель, и в то же время увеличение детализации, если выбрана какая-то меньшая часть модели.

Отладка


Далее следовал этап оптимизации. Для одного объекта томограф выдаёт 360 снимков, каждый разрешением до 8000*8000. Поскольку объём обрабатываемых данных велик, решение задачи «в лоб» было бы совершенно неудовлетворительно. Это учитывалось ещё на этапе проектирования, однако после получения первой версии алгоритмы пришлось несколько раз оптимизировать и адаптировать. Задание требовало, чтобы время создания трехмерной модели микроструктуры не превышало 2-х часов, поэтому этап оптимизации был заложен изначально. В ходе тестирования первой версии мы столкнулись с тем, что использование стандартного формата хранения изображений проекций не подходит для проекта. Входные данные содержат изображения формата TIFF с 16-битным кодированием оттенков серого. Такая глубина цвета для расчётов избыточна, а дискового пространства, сетевого канала, оперативной памяти и процессорного времени обработка таких изображений требует много. С другой стороны, стандартной 8-битной глубины цвета нам оказалось недостаточно для сохранения точности реконструкции. Поэтому был разработан формат хранения изображений с 12-битной глубиной цвета.

В технический проект было заложено горизонтальное масштабирование вычислений. Реконструкция 3D-модели, то есть основная вычислительная задача, разделялась на небольшие пакеты-задания, которые центральный модуль ПО распределял по сети серверов в кластере. На серверах применялась технология CUDA, позволяющая задействовать вычислительные мощности графических процессоров для расчётов. Время, требуемое для расчёта одной модели, сокращается пропорционально количеству серверов в кластере, так как вычислительные задачи идеально распараллеливаются, и все серверы загружены на 100%.

image

Архитектура CUDA применима не только для высокопроизводительных графических вычислений, но и для различных научных вычислений с использованием видеокарт nVidia. Ученые и исследователи широко используют CUDA в различных областях, включая астрофизику, вычислительную биологию и химию, моделирование динамики жидкостей, электромагнитных взаимодействий, компьютерную томографию, сейсмический анализ и многое другое. В CUDA имеется возможность подключения к приложениям, использующим OpenGL и Direct3D. CUDA — кроссплатформенное программное обеспечение для таких операционных систем как Linux, Mac OS X и Windows.

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


Карты, поддерживаемые нашим ПО:


Задание предусматривало, что ПО должно работать в среде Microsoft Windows XP/Vista/7, Linux. В связи с этим кроссплатформенность решения была заложена изначально. В качестве языка разработки был выбран C++/Qt, что позволило иметь единый исходный код и собирать ПО под разные ОС.





Видео























Больше проектов:
Как за 5233 человеко-часа создать софт для микротомографа
SDK для внедрения поддержки электронных книг в формате FB2
Управление доступом к электронным документам. От DefView до Vivaldi
Подробнее о разработке софта рентгеновского томографа
«Сфера»: как мониторить миллиарды киловатт-часов
Разработка простого плагина для JIRA для работы с базой данных
В помощь DevOps: сборщик прошивок для сетевых устройств на Debian за 1008 часов
Автообновление службы Windows через AWS для бедных
Автор: @chookcha
Edison
рейтинг 75,72
Изобретаем успех: софт и стартапы

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

  • +4
    Впечатляющая работа!!!
  • 0
    а погонять эту программу можно?
  • +6
    С удовольствием бы почитал про обращение преобразования Радона с точечным источником и о реализации этого всего с использованием CUDA
  • +1
    Извините, не понял. Это ваш продукт или вы сделали эту прогу для кого то? Сколько стоит сам томограф и прога к нему?
  • +5
    Вот это вещь! Ролики впечатлили. Интересно, а разглядеть что изображено под слоем «сотри и выиграй» лотерейного билета на нем можно?
  • +1
    Мде интересно, а лотерейные билеты (которые ногтем тереть) эта машина просветит? Там же должен быть рельеф?
    • +3
      Забавный hivemind у вас и товарища выше.
    • +2
      Отличие в плотности цифр под слоем минимальная, рентген не возьмёт. тем более их как раз и разрабатывали в том числе и для того чтобы нельзя было просветить рентгеном. Структура билетиков многослойная и в каждом слое свой мусор, который не отличается материалом от самих цифр — на фоне этого мусора на рентгеновском снимке надпись просто не будет видно.
      И ещё, я бы предположил что краска которую надо стирать меняет цвет после рентгена. Это было бы логичным шагом.
      • 0
        Спасибо :(
        Что ж… значит придется по старинке: запираться в сортире с светящим мобильником и пультом от теле. Долго-долго смотреть… обычно после женщины в красном что-то, да проглядывает через черную краску.
      • 0
        Томограмма могла бы побороться с многослойностью, будь разрешение достаточным
        • 0
          У меня есть сомнения на этот счёт. Эти билетики делают специально чтобы было трудно просвечивать, куда не ткни, а луч вынужден проходить множество неоднородностей из такого же материала как и сама надпись. С какой стороны леса не посмотри, а на опушку в глубине леса никак не посмотришь.
          • 0
            Эй, это же томограмма! Это и есть способ смотреть «со всех сторон леса на любую опушку в глубине леса»
            • 0
              А если опушка заросла такими же деревьями, другими по форме но точно такими же по проникающей способности для рентгена?
  • +6
    Решили подобную задачу (для большого томографа) в 2000 году. Отличие — пациент облучался не 360 раз, а 4. По 4 снимкам восстанавливали трехмерную картинку исходя из принципа максимальной энтропии. Решали методом Метрополиса. Изоповерхность строили так же, но вокселов не знали, картинки были статичные, строились медленно.
    • 0
      Подробности есть? В продакшн пошло? За снижение лучевой нагрузки все томографостроители борятся же!
      • 0
        Заказчик — медицинская контора из Бостона. Проект в рамках МНТЦ. На нас было лишь R&D. Про конечный железный продукт ничего не могу сказать.
        • 0
          Алгоритмы насколько выкуплены и защищены? Есть право / возможность выложить в open source? Хотя бы в виде статьи с формулами?
          • +1
            Source code принадлежит заказчику. Было два доклада и 1 публикация в ВАНТ (Вопросы Атомной Науки и Техники), пытался найти ссылку — не нашел. Нашел только выпуск про построение изоповерхности восстановленных данных от томографа, что тривиально в наше время. Спрошу у соавтора, если он еще жив.
      • 0
        Ну не все же борятся. Есть ведь медицина, а есть промышленность.
        Ближайший аналог того томографа, что в статье, это, к примеру Phoenix Nanotom, он в серийном продакшене:
        image
        У него паспортное разрешение до 0,2 микрона, но, правда и стоимость соответствующая. Вот видео

        На самом деле алгоритмы реконструирования сами по себе хорошо известны и относительно несложные, однако дьявол кроется в деталях — если применять их «в лоб», то мы получим артефакты преобразования (в виде концентрических окружностей, полос, и т.п.) и довольно большая часть не кода не столько сама реконструкция, сколько борьба с артефактами. Ну и стабильность железа важна — зачастую всё монтируют на монолитной гранитной плите, ставят подшипники на воздушной подушке, термостабилизированный детектор, встроенный кондиционер и всё такое.
        • 0
          Шикарно. Даже не думал, почему-то, что это ещё и 3D сканеры ведь, идеально подходящие для дальнейшей 3D печати!
          • 0
            А, кстати, мысль. Ведь воспроизвести можно не только внешний контур, но и всю внутреннюю структуру. Пожалуй 3D принтер и станет моей следующей игрушкой, надо только жабу задушить.
  • 0
    Фантастика.
  • +2
    Здорово, просто молодцы. Хотя у меня по ходу чтения появилось немножко вопросов.

    Камера, которую вы использовали (PIXIS-XF надо полагать) отдаёт картинки 2048х2048, а в статье вы пишете «до 8000х8000». Это вы как получили? Вы перемещаете также образец или камеру по вертикали и горизонтали и склеиваете потом картинки?

    Демо изображения, которые в видео в конце статьи — они все получены из 360 проекций? Если так, то это хорошо, ведь 360 проекций с шагом в градус — довольно мало, обычно идут с шагом треть/четверть градуса, иначе будут артефакты реконструкции. Вроде есть формула оптимального количества проекций для заданного разрешения, но вот запамятовал.

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

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

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

    Спасибо.
  • 0
    Смотрелку сохранённого воксельного массива не хотите выложить в исходниках и сделать расширяемой энтузиастами? А то все они убогие и закрытые…
  • 0
    Коллеги, ответы на вопросы можно почитать в новой статье: Подробнее о разработке софта рентгеновского томографа.

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

Самое читаемое Разработка