Хочу поподробнее рассказать об интересном проекте компании Edison. Перед разработчиками поставили задачу написать софт для микротомографа, они с этим отлично справились, а потом запихивали в этот томограф семечки, болты, конденсаторы и моль. А серьезным дядям этот томограф нужен, чтобы проверять алмазы и не покупать дырявые.
А еще сегодня 16 декабря, день рождения Иоганна Радона, австрийского математика, ректора Венского университета, который в 1917 году ввел интегральное преобразование функции многих переменных, родственное преобразованию Фурье, используемое сегодня во всех томографах.
Иоганн Радон был профессором 6 университетов (а в одном из них даже без кафедры), был президентом Австрийского математического общества. В Австрии в честь него назвали «Институт вычислительной и прикладной математики» и медаль.
О том, как проходила разработка софта для томографа и какие задачи решались в процессе — под катом.
Микротомограф
Ученые из томского государственного университета создали микротомограф. Он позволяет с точностью до микрона узнать о внутренней структуре различных материалов. Например, алмазов.
Томограф может просветить материал с разрешением до микрона. Это в 100 раз тоньше человеческого волоса. После сканирования программа создает 3D-модель, где можно посмотреть не только на внешнюю сторону детали, но и узнать, что у нее внутри.
Используемые математические алгоритмы.
- Обратное преобразование Радона.
- Марширующие/шагающие кубы.
- Фильтр Гаусса.
- Фильтрация/свертка, нормализация проекций.
Реализация и технологии.
- C++/Qt.
- Ubuntu, Windows.
- CUDA.
- Объемный-воксельный рендеринг модели.
- Собственный формат хранения изображений с оттенками серого глубиной 12 бит.
- Разделение реконструкции на клиент-сервер.
- Сохранение объемных данных в виде «Октодеревьев».
Трудозатраты: 5233 человеко-часа.
Отладка производилась на сырых данных (проекционных снимках), полученных с томографа.
Алгоритмы
В первую очередь, необходимо было подобрать математический аппарат для решения задачи. Основной алгоритм — обратное преобразование Радона — был заложен в постановку задачи,
однако его потребовалось адаптировать под особенности нашей работы и использовать несколько дополнительных, вспомогательных алгоритмов. Например, ввиду того, что объект просвечивался одной «лампочкой», пришлось адаптировать формулы обратного преобразования Радона под конические, а не прямые проекции. Стандартный алгоритм подразумевает, что объект просвечивается пучком параллельных рентгеновских лучей, исходящих от бесконечно удалённого источника. В реальности источник лучей точечный, поэтому и пучок лучей имеет форму конуса. В связи с этим в алгоритм обратного преобразования Радона потребовалось ввести преобразование координат из конической системы в прямоугольную.
На первой стадии вычислений производится предварительная фильтрация/свертка, нормализация проекций. Это необходимо для того, чтобы на проекциях приглушить шумы, а плотности более явно выделить.
Для построения 3D-моделей поверхностей объектов в стандартном формате для просмотра в 3D-редакторах Компас, SolidWork, 3D Max Studio использовался алгоритм «Марширующие (шагающие) кубы». Суть алгоритма в том, что он пробегает скалярное поле, на каждой итерации просматривает 8 соседних позиций (вершины куба, параллельного осям координат) и определяет полигоны, необходимые для представления части изоповерхности, проходящей через данный куб. Далее, на экран выводятся полигоны, образующие заданную изоповерхность.
Под фильтром Гаусса в проекте понимается матричный фильтры обработки изображений с использованием матрицы свертки. Матрица свёртки – это матрица коэффициентов, которая «умножается» на значение пикселей изображения для получения требуемого результата. Фильтр используется для сглаживания воксельных данных и проекций срезов, что, в свою очередь, позволяет улучшить качество генерируемых 3D-моделей (Хабрапост 1, Хабрапост 2).
Реализация и технологии
В ходе работы также был создан ряд специализированных технических решений: библиотека для объемно-воксельного рендеринга модели; запись видео при выполнении операций с моделью; собственный формат хранения изображений с оттенками серого глубиной 12 бит; сохранение объемных данных в виде «Октодеревьев»; алгоритмы полигонизации 3D-модели.
Объемный-воксельный рендеринг модели в проекте использовался для просмотра модели с возможностью вращения и масштабирования в реальном времени. Воксель — это трехмерный пиксель. Так же с помощью воксельного рендеринга оператору предоставляются удобные инструменты для определения области просмотра с автоматическим увеличением уровня детализации и возможность расположения секущей плоскости под любым углом в два клика. На основе секущей плоскости в дальнейшем можно получить изображение среза с максимальным разрешением.
Октодерево (дерево октантов, восьмеричное дерево, англ. octree) — тип древовидной структуры данных, в которой у каждого внутреннего узла ровно восемь «потомков». Восьмеричные деревья чаще всего используются для разделения трёхмерного пространства, рекурсивно разделяя его на восемь ячеек. В проекте октодерево позволяет выполнять отображение данных в режиме предварительного просмотра, когда нет необходимости в данных, которые не видны пользователю. Так, например, объемный рендеринг получает набор данных для отображения с детализацией в зависимости от выбранной области с помощью октодерева, что обеспечивает высокий FPS, когда видна вся модель, и в то же время увеличение детализации, если выбрана какая-то меньшая часть модели.
Отладка
Далее следовал этап оптимизации. Для одного объекта томограф выдаёт 360 снимков, каждый разрешением до 8000*8000. Поскольку объём обрабатываемых данных велик, решение задачи «в лоб» было бы совершенно неудовлетворительно. Это учитывалось ещё на этапе проектирования, однако после получения первой версии алгоритмы пришлось несколько раз оптимизировать и адаптировать. Задание требовало, чтобы время создания трехмерной модели микроструктуры не превышало более 2-х часов, поэтому этап оптимизации был заложен изначально. В ходе тестирования первой версии мы столкнулись с тем, что использование стандартного формата хранения изображений проекций не подходит для проекта. Входные данные содержат изображения формата TIFF с 16-битным кодированием оттенков серого. Такая глубина цвета для расчётов избыточна, а дискового пространства, сетевого канала, оперативной памяти и процессорного времени обработка таких изображений требует много. С другой стороны, стандартной 8-битной глубины цвета нам оказалось недостаточно для сохранения точности реконструкции. Поэтому был разработан формат хранения изображений с 12-битной глубиной цвета.
В технический проект было заложено горизонтальное масштабирование вычислений. Реконструкция 3D-модели, то есть основная вычислительная задача, разделялась на небольшие пакеты-задания, которые центральный модуль ПО распределял по сети серверов в кластере. На серверах применялась технология CUDA, позволяющая задействовать вычислительные мощности графических процессоров для расчётов. Время, требуемое для расчёта одной модели, сокращается пропорционально количеству серверов в кластере, так как вычислительные задачи идеально распараллеливаются, и все серверы загружены на 100%.
Архитектура CUDA применима не только для высокопроизводительных графических вычислений, но и для различных научных вычислений с использованием видеокарт nVidia. Ученые и исследователи широко используют CUDA в различных областях, включая астрофизику, вычислительную биологию и химию, моделирование динамики жидкостей, электромагнитных взаимодействий, компьютерную томографию, сейсмический анализ и многое другое. В CUDA имеется возможность подключения к приложениям, использующим OpenGL и Direct3D. CUDA — кроссплатформенное программное обеспечение для таких операционных систем как Linux, Mac OS X и Windows.
В нашем проекте CUDA применяется для основного процесса — реконструкции объёмных данных из проекций. Так как графические процессоры имеют специализированный набор команд, то вычисления реконструкции хорошо ложились именно на графические карты через технологию CUDA. На CPU данная задача решается дольше как на этапе разработки, так и на этапе выполнения.
Карты, поддерживаемые нашим ПО:
- Nvidia Tesla K80 24GB (научная);
- EVGA GeForce GTX TITAN X 12GB (игровая).
Задание предусматривало, что ПО должно работать в среде Microsoft Windows XP/Vista/7, Linux. В связи с этим кроссплатформенность решения была заложена изначально. В качестве языка разработки был выбран C++/Qt, что позволило иметь единый исходный код и собирать ПО под разные ОС.
Видео
Комментарии (25)
Uranix
16.12.2015 20:53+6С удовольствием бы почитал про обращение преобразования Радона с точечным источником и о реализации этого всего с использованием CUDA
schroeder
16.12.2015 22:16+1Извините, не понял. Это ваш продукт или вы сделали эту прогу для кого то? Сколько стоит сам томограф и прога к нему?
inwardik
16.12.2015 22:27+5Вот это вещь! Ролики впечатлили. Интересно, а разглядеть что изображено под слоем «сотри и выиграй» лотерейного билета на нем можно?
nikitasius
16.12.2015 23:11+1Мде интересно, а лотерейные билеты (которые ногтем тереть) эта машина просветит? Там же должен быть рельеф?
Alexeyslav
17.12.2015 10:29+2Отличие в плотности цифр под слоем минимальная, рентген не возьмёт. тем более их как раз и разрабатывали в том числе и для того чтобы нельзя было просветить рентгеном. Структура билетиков многослойная и в каждом слое свой мусор, который не отличается материалом от самих цифр — на фоне этого мусора на рентгеновском снимке надпись просто не будет видно.
И ещё, я бы предположил что краска которую надо стирать меняет цвет после рентгена. Это было бы логичным шагом.nikitasius
17.12.2015 11:17Спасибо :(
Что ж… значит придется по старинке: запираться в сортире с светящим мобильником и пультом от теле. Долго-долго смотреть… обычно после женщины в красном что-то, да проглядывает через черную краску.
Nashev
18.12.2015 10:25Томограмма могла бы побороться с многослойностью, будь разрешение достаточным
Alexeyslav
18.12.2015 10:54У меня есть сомнения на этот счёт. Эти билетики делают специально чтобы было трудно просвечивать, куда не ткни, а луч вынужден проходить множество неоднородностей из такого же материала как и сама надпись. С какой стороны леса не посмотри, а на опушку в глубине леса никак не посмотришь.
Nashev
18.12.2015 11:02Эй, это же томограмма! Это и есть способ смотреть «со всех сторон леса на любую опушку в глубине леса»
Alexeyslav
18.12.2015 11:57А если опушка заросла такими же деревьями, другими по форме но точно такими же по проникающей способности для рентгена?
PapaBubaDiop
17.12.2015 00:45+6Решили подобную задачу (для большого томографа) в 2000 году. Отличие — пациент облучался не 360 раз, а 4. По 4 снимкам восстанавливали трехмерную картинку исходя из принципа максимальной энтропии. Решали методом Метрополиса. Изоповерхность строили так же, но вокселов не знали, картинки были статичные, строились медленно.
Nashev
18.12.2015 10:29Подробности есть? В продакшн пошло? За снижение лучевой нагрузки все томографостроители борятся же!
PapaBubaDiop
18.12.2015 15:22Заказчик — медицинская контора из Бостона. Проект в рамках МНТЦ. На нас было лишь R&D. Про конечный железный продукт ничего не могу сказать.
Nashev
18.12.2015 15:52Алгоритмы насколько выкуплены и защищены? Есть право / возможность выложить в open source? Хотя бы в виде статьи с формулами?
PapaBubaDiop
22.12.2015 23:20+1Source code принадлежит заказчику. Было два доклада и 1 публикация в ВАНТ (Вопросы Атомной Науки и Техники), пытался найти ссылку — не нашел. Нашел только выпуск про построение изоповерхности восстановленных данных от томографа, что тривиально в наше время. Спрошу у соавтора, если он еще жив.
AndreyDmitriev
18.12.2015 18:53Ну не все же борятся. Есть ведь медицина, а есть промышленность.
Ближайший аналог того томографа, что в статье, это, к примеру Phoenix Nanotom, он в серийном продакшене:
У него паспортное разрешение до 0,2 микрона, но, правда и стоимость соответствующая. Вот видео
На самом деле алгоритмы реконструирования сами по себе хорошо известны и относительно несложные, однако дьявол кроется в деталях — если применять их «в лоб», то мы получим артефакты преобразования (в виде концентрических окружностей, полос, и т.п.) и довольно большая часть не кода не столько сама реконструкция, сколько борьба с артефактами. Ну и стабильность железа важна — зачастую всё монтируют на монолитной гранитной плите, ставят подшипники на воздушной подушке, термостабилизированный детектор, встроенный кондиционер и всё такое.Nashev
18.12.2015 19:43Шикарно. Даже не думал, почему-то, что это ещё и 3D сканеры ведь, идеально подходящие для дальнейшей 3D печати!
AndreyDmitriev
18.12.2015 19:58А, кстати, мысль. Ведь воспроизвести можно не только внешний контур, но и всю внутреннюю структуру. Пожалуй 3D принтер и станет моей следующей игрушкой, надо только жабу задушить.
AndreyDmitriev
17.12.2015 18:47+2Здорово, просто молодцы. Хотя у меня по ходу чтения появилось немножко вопросов.
Камера, которую вы использовали (PIXIS-XF надо полагать) отдаёт картинки 2048х2048, а в статье вы пишете «до 8000х8000». Это вы как получили? Вы перемещаете также образец или камеру по вертикали и горизонтали и склеиваете потом картинки?
Демо изображения, которые в видео в конце статьи — они все получены из 360 проекций? Если так, то это хорошо, ведь 360 проекций с шагом в градус — довольно мало, обычно идут с шагом треть/четверть градуса, иначе будут артефакты реконструкции. Вроде есть формула оптимального количества проекций для заданного разрешения, но вот запамятовал.
Ещё я не очень понял про частоту камеры. По спецификации она двухчастотная на 100 килогерц или два мегагерца. Если у неё четыре мегапиксела — это значит, что она отдаёт кадр каждые две секунды на максимальной частоте?
Вы перемещаете манипулятор пошагово или непрерывно? Сколько времени занимает сканирование типичных образцов, что представлены на видео?
Ну и про 12 бит — очень любопытно. Ну, то что отрезать четыре бита и упаковать каждые два пиксела в три байта можно для хранения и передачи — это понятно. Но для реконструирования вам же придётся снова развернуть каждый пиксел как минимум в два байта? Или у вас вся математика на 12-ти битах — в этом случае как вы решили проблему того, что пикселы занимают полтора байта и не выравниваются на границу?
Спасибо.
Nashev
18.12.2015 10:41Смотрелку сохранённого воксельного массива не хотите выложить в исходниках и сделать расширяемой энтузиастами? А то все они убогие и закрытые…
viktorpanasiuk
Впечатляющая работа!!!