Результаты тестирования программы компьютерной томографии
vkUtility/Science/CTomo
Автор: Виктор Кон (http://kohnvict.ucoz.ru/main.htm). Только для посвященных.
Введение в проблему
Метод компьютерной томографии (КТ) состоит в том, чтобы получить распределение плотности в сечении объекта на определенной высоте по набору проекций, измеренных под разными углами, в каждой из которых записана интегральная вдоль лучей плотность. Метод является неразрушающим и получил исключительно широкое распространение в медицине. Авторы метода А. Кормак и Г. Хаунсфилд даже получили в 1979 году Нобелевскую премию по физиологии и медицине.
Справедливости ради следует отметить, что математический алгоритм для КТ разработал И. Радон еще в 1917 году, но тогда не было ни хороших рентгеновских источников, ни компьютеров. Метод КТ также широко используется в физике для диагностики структуры различных материалов. Существует несколько алгоритмов реализации метода. Наиболее распространены метод FBP (filtered back projection) и алгебраический метод. Все методы томографии давно и хорошо описаны, например, в книге
Avinash C. Kak, Malcolm Slaney "Principles of Computerized Tomographic Imaging", IEEE Press, 1988
C 1999 года книгу в электронном виде можно скачать по этой ссылке 
. Там представлено оглавление книги, и название каждой главы представляет собой ссылку на скачивание pdf файла. По данным сайта Google Академии эта книга на май 2016 года цитировалась 6832 раза.
В 2005 году меня попросили сделать программу томографии, удобную для чайников, и я, как примерный пионер, ее сделал. Дело в том, что коммерческие программы стоят дорого, а любительские, как правило, не удобные. Одновременно я записал небольшой pdf файл о программе, в котором кратко изложил основы метода FBP и детали программной реализации. Вот ссылка на этот файл . Файл необходимо переписать, но пока что есть, то есть.
В то время программа так и не была использована, и работа оказалась невостребованной. Несколько лет назад я переписал программу под другой набор данных в виде tiff файлов, и программа стала использоваться на Курчатовском источнике синхротронного излучения. Сейчас я решил включить ее в большой набор программ vkUtility, и заодно оптимизировать как методы работы, так и внешний вид, а также протестировать на разных моделях.
Описание метода тестирования и первый пример
Самый лучший способ проверки правильности работы программы -- это использование модельных объектов, структура которых известна. Этот способ предполагает на первом этапе расчет экспериментальных проекций от известного объекта. На втором этапе проекции используются для восстановления структуры объекта с последующим сравнением результата реконструкции с оригиналом. Именно сравнение с оригиналом позволяет понять насколько правильно работает программа.
Сложность в том, что необходимо очень точно приготовить экспериментальные проекции из оригинальной картинки. Если это делать не очень точно, то восстановленная картинка с оригиналом не совпадет. В качестве первого объекта в данной статье я рассмотрю фантом Shepp-Logan, который часто используется для тестирования алгоритмов КТ. Я впервые увидел его в учебной программе LearnCT, которую можно найти в интернете по названию. Родной сайт программы , на котором она стоит с 2006 года, сейчас почему-то не запускается.
Я нашел аналогичную картинку в интернете, привел ее к размеру 512*512, и затем приготовил проекции. Для этой цели я из картинки вытащил значения пикселей в интервале от 0 до 255, которые образуют матрицу размером 512*512. Эту матрицу я вращал на разные углы вокруг ее центра. При этом формировалась новая матрица точно такого же размера 512*512. Для каждой точки новой матрицы находился ее аналог в исходной матрице, то есть координаты, которые она бы получила после обратного вращения. И затем значение в этой точке получалось интерполяцией по 4-м ближайшим точкам исходной матрицы.
Все матрицы (для разных углов вращения) суммировались по вертикальной координате, и таким образом получались проекции размером 512 точек, которые записывались в файл формата tiff после умножения на коэффициент -M/(D*255), вычисления экспоненты и умножения на большое число, меньшее 65536, которое является максимальным для записи в tiff файлах. Здесь D=512 есть линейный размер проекции. При этом в тех точках проекции, где материала не было совсем, в сумме останется чистый ноль, экспонента будет равна 1, а в проекции будет максимальное число Fmax. А там, где материал был, будут числа меньше максимального. До нуля не дойдет, но могут быть и маленькие числа, если параметр M взять достаточно большим.
Так как детекторы обычно записывают свои данные в картинки формата tiff как двухбайтовые целые числа без знака, то максимальное число, которое можно записать равно 65535. То есть можно умножать на любое число меньше этого. Величины M и Fmax являются параметрами расчета, но результат реконструкции от них практически не зависит. Я выбрал M=1, Fmax=32500. Число проекций N=360, то есть угол менялся с шагом на половину градуса.
Программа расчета проекций каждую проекцию записывала в отдельный файл формата tiff так, как делают детекторы, то есть она записала 360 файлов с номерами от 000 до 359. На этом этап моделирования экспериментальных данных закончился. Программа, которая это делает, не является программой томографии. Она написана отдельно и не имеет нужного интерфейса с пользователем.
Затем в работу вступила основная программа томографии, которая тестируется. Она сделана так, что имеет три независимых модуля расчета. Первый модуль вычисляет синограммы для всех сечений трехмерного объекта по высоте, если объект имеет высоту в направлении оси вращения. Синограмма -- это матрица размерности D*N, в которой все проекции заданного сечения лежат как бы одна над другой. То есть программа выбирает из каждого файла двумерной проекции (если объект имеет высоту) одну строку на заданной высоте, и записывает ее в строку синограммы с вертикальной координатой, соответствующей углу проекции. Эту матрицу она записывает в файл, но одновременно определяет минимальное и максимальное значения матрицы. Их она записывает в другой текстовый файл с названием minmax.txt.
Эта операция выполняется относительно быстро, и ее необходимо выделить для того, чтобы определить самое максимальное значения на всех синограммах. На втором этапе расчета это максимальное значение надо указать в форме ввода входных данных. Правильное задание этого значения гарантирует получение минимального фона на томограмме в тех местах, где нет материала. Программа вычисляет логарифм максимального значения и потом она вычитает из этого значения логарифм всех чисел синограммы. Соответственно, для максимального значения разница будет равна нулю.
Итак, на первом этапе программа записывает синограмму, а затем из нее получает томограмму в виде картинки точно такого же размера и типа, как и оригинал. Дополнительная оптимизация состоит в том, что отрицательные значения в томограмме заменяются на ноль, потому что плотность не может быть отрицательной. Ниже слева показана картинка оригинала, а справа томограмма.
Рис. 1. Тест с медицинским фантомом Shepp-Logan для 360 проекций. Слева оригинал, справа томограмма.
После ряда оптимизаций кода время расчета томограммы равно 11 сек на ноутбуке Dell (процессор Intel Сore-i5-3210M 2.50 Hz, оперативная память 8 гб, ). Следует отметить, что размер оперативной памяти в задачу не входит, лишь бы ее хватило. Видно, что томограмма почти совпадает с оригиналом, как по форме объектов, так и по цвету. Единственное очень слабое отличие можно заметить в цвете верхнего светлого эллипса.
Программа позволяет по тому же набору данных провести расчет с меньшим числом проекций. Я задал число проекций 180 и 90, чтобы понять как портится картинка при уменьшении числа проекций. Результаты показаны на рисунке 2. Интересно, что однородные области становятся неоднородными, но цвет областей в целом сохраняется. И сохраняются резкие границы областей. Если провести обработку изображения разными сглаживающими фильтрами, то можно потерять резкость границ. Время расчета прямо пропорционально числу проекций.
Рис. 2. Топограммы медицинского фантома Shepp-Logan для 180 проекций (слева) и 90 проекций(справа).
Интересно посмотреть что происходит при меньшем числе проекций. Это показано на рисунке 3. Видно, что даже при очень малом числе проекций крупные объекты картинки все равно видны, а больше всего искажается фон за пределами светлого контура. Результаты можно признать более, чем удовлетворительными, но для получения совпадения весьма критическим является правильное приготовление проекций из картинки. Также важно, что фантом черепа с ярким контуром по границе объекта является очень устойчивым к восстановлению.
Рис. 3. Топограммы медицинского фантома Shepp-Logan для 45 проекций (слева) и 24 проекций(справа).
Тесты для различных моделей другого типа
В качестве более сложного модельного объекта я взял картинку острова в океане, точнее фотографию из космоса, которая показана на сайте go2life.net. Вот ссылка на картинку 
. С помощью программы FastStone Image Viewer я превратил ее в серую, вырезал фрагмент размером 512*512 пикселей, и затем уже с помощью собственной ACL программы сделал клип в виде круга, за пределами которого был просто черный фон.
В результате получилась картинка, показанная на рисунке 4 слева. Далее по описанной выше схеме, я записал для нее проекции и провел реконструкцию. Результат для 360 проекций показан на том же рисунке справа. Снова отличное совпадение. Потрясает не только отличное совпадение всех контуров деталей, но и полное соответствие цветов. Также можно сделать вывод, что даже для картинки столь сложной структуры вполне достаточно 360 проекций для полного восстановления всех деталей.
Рис. 4. Тест с фотографией острова в океане, слева оригинал, справа томограмма из 360 проекций.
Интересно проверить модель с минимумом цветов, но с большим числом резких границ мелких объектов. Такой моделью может быть текст в круге. Результаты тестирования такой модели показаны на рисунке 5. Слева, как и раньше, оригинал, справа томограмма на 360 проекций.
Рис. 5. Тест с фотографией печатного текста, слева оригинал, справа томограмма из 360 проекций.
На этот раз совпадение оставляет желать лучшего. Хотя текст отлично виден и легко читается., но цвет фона под ним не белый, а за пределами круга цвет не черный. То есть на двухцветной (черно-белой) модели объекта неточность FBP метода особенно видна. На рисунке 6 показаны результаты реконструкции из большего (720) и меньшего (180) числа проекций. При этом увеличение числа проекций не решает проблему, а уменьшение приводит к еще большим искажениям.
Рис. 6. Томограммы печатного текста, слева из 720 проекций, справа из 180 проекций.
Такой результат можно объяснить тем, что в методе FBP делается преобразование Фурье с помошью процедуры FFT (Fast Fourier Transformation). Эта процедура в местах резких скачков функции приводит к осцилляциям фурье-образа с высокой частотой. Эти частоты могут выходить за пределы конечной области расчетной сетки, и приводить к сглаживанию скачков. Более того, simple ramp фильтр, который используется при обратном проектировании, вероятно, только усиливает размазывание черного цвета по области рисунка.
Как улучшить качество реконструктции сечений такого типа я пока не знаю. Возможно, это самый сложный тип структуры для метода FBP, а возможно и для других методов. Я проводил расчет и для большего числа точек сечения, например, 1024. Ситуация остается такой же, то есть сохранить белый фон за текстом не удается, а цвет области за кругом тоже не воспроизводится.
Я также провел расчет на ультрабуке Самсунг (процессор Intel Сore-i5-2537M 1.40 Hz, оперативная память 4 гб) и с JRE-32. Расчеты на Dell проводились с JRE-64. Время счета было примерно таким же, лишь чуть чуть больше. Получается так, что время даже не зависит от версии JRE. Однако, вот интересный момент. Это время Самсунг показывает при зарядке от сети. А если работать на аккумуляторе, то время возрастает в полтора раза. Вероятно у процессора Intel Core-i5 есть разные режимы работы. Режим от аккумулятора на Dell я не проверял.
Рис. 7. Томограммы печатного текста из 720 проекций, с коррекцией синограмм с краев: 50 точек (слева) и 100 точек (справа).
У программы есть дополнительные возможности, к которым относится коррекция синограмм с краев области. Параметр Background correction. Предполагается, что объект имеет конечные размеры и по краям каждой проекции интенсивность пучка максимальна, а синограмма имеет нулевые значения. Если же это не так, то программа позволяет принудительно обнулить синограмму с краев области при ее наборе из проекций. Можно задать число точек, которые необходимо обнулить.
Результаты такого расчета показаны на рисунке 7. Слева обнуляется 50 точек с каждого края, справа -- 100 точек. Не очень понятно как интерпретировать результаты. С одной стороны, обнуление синограмм дает более черный фон за пределами информационной области. Но и внутри этой области фон тоже становится более черным. А хотелось бы получить обратный результат. С другой стороны, в этом случае граница информационной области получается очень белой, и картинку очень легко скорректировать выбирая контрасть не из всего интервала значений, а только из части.
Рис. 8. Томограмма рисунка 7 с коррекцией контраста. Черный цвет -- минимум, белый -- 30% максимума (слева) и 20% (справа)
На рисунке 8 показаны два варианта такого представления. Видно, что в этом случае контраст сильно уменьшается, а текст виден не только в пределах информационного круга, но и за его пределами, хотя и в искаженной форме. То есть обнуление краев синограммы не убивает полностью текст за пределами круга. Более того, он оказывается даже более черным, чем фон. Можно еще и снизу обрезать интервал значений. И получить чистую черно-белую картинку. То же самое можно сделать и с исходной картинкой, то есть при нулевом обрезании синограммы.
Вообще говоря, такая коррекция возможно только тогда, когда знаешь что нужно получить. А если структура объекта заранее не известна, то такие эксперименты только показывают, что возможно восстановление части большого объекта даже в том случае, когда размер проекции меньше размеров объекта. Такая ситуация может иметь место в большом числе случаев. С другой стороны объекты типа черно-белого текста редко встречаются в природе.
Есть и еще один способ увеличения контраста в том случае, если размер объекта меньше размеров проекции. Он состоит в том, что программа позволяет вырезать часть числовой матрицы из полной матрицы томограммы и нарисовать ее с максимальным контрастом. Пример использования программы в таком режиме показан на рисунке слева от этого текста. Небольшой фрагмент вырезан из левой томограммы рисунка 7. Видно, что контраст текста значительно выше, чем на исходном рисунке.
Так обычно поступают для того, чтобы не показывать части томограммы, которые не содержат полезной информации, в частности границу круга. Не очень понятно как это согласуется с трехмерной реконструкцией из сечений. Но с точки зрения увеличения контраста -- это хороший способ.