Новая версия моей программы компьютерной томографии
Виктор Кон, . . . 1 апреля 2022 года, . . . http://kohnvict.narod.ru
ВВЕДЕНИЕ.
Метод компьютерной томографии (КТ) состоит в том, чтобы получить распределение плотности некоторого параметра (например, массы) в сечениях объекта на разных уровнях высоты по набору двумерных проекций, измеренных под разными углами, в каждой из которых записана интегральная вдоль лучей плотность в горизонтальном направлении. Метод является неразрушающим и получил исключительно широкое распространение в медицине. Авторы метода А. Кормак и Г. Хаунсфилд даже получили в 1979 году Нобелевскую премию по физиологии и медицине.
Справедливости ради следует отметить, что математический алгоритм для КТ разработал И. Радон еще в 1917 году, но тогда не было ни хороших рентгеновских источников, ни компьютеров. Метод КТ также широко используется в физике для диагностики структуры различных материалов. Существует несколько алгоритмов реализации метода. Наиболее распространены метод FBP (filtered back projection) и алгебраический метод. Все методы томографии давно и хорошо описаны, например, в книге
Avinash C. Kak, Malcolm Slaney "Principles of Computerized Tomographic Imaging", IEEE Press, 1988
Книгу в электронном виде можно скачать по этой ссылке 
. Там представлено оглавление книги, и название каждой главы представляет собой ссылку на скачивание pdf файла. По данным сайта Google Академии эта книга на дату написания этой статья цитировалась 10352 раза.
В 2005 году знакомый экспериментатор из старшего поколения попросил меня сделать программу томографии, удобную для чайников, и я, как примерный пионер, ее сделал. Дело в том, что коммерческие программы стоят дорого, а любительские, как правило, не удобные и не пригодны для серьезной работы. Одновременно я записал небольшой pdf файл о программе, в котором кратко изложил основы метода FBP и детали программной реализации. Вот ссылка на этот файл . Файл было бы неплохо переписать, но пока что есть, то есть.
В то время программа так и не была использована, и работа оказалась невостребованной. Несколько лет назад я переписал программу под другой набор данных в виде tiff файлов, и программа некоторое время использовалась на Курчатовском источнике синхротронного излучения. Одной из наиболее широко известных коммерческих программ КТ является программа Октопус (Octopus). Возможно, одной из первых научных публикаций о ней была эта . Статья вышла в 2004 году. Программа с тех пор развивалась, потом стала продаваться за деньги и до сих пор продается, но уже не так удобно как раньще. Учитывая сложившуюся к 2022 году ситуацию получить ее из России стало весьма проблематично, хотя и раньше было не просто.
С другой стороны, я как раз получил возможность обработать кое какие реальные экспериментальные данные с помощью своей программы и сравнить свои результаты с результатами, полученными с помошью программы Октопус. Оказалось, что в мою программу надо добавить методы обработки данных как до, так и после применения программы. Я также добавил кое-какие изменения и в саму программу. Так появилась новая версия программы, которая и рассматривается в данной статье.
НЕМНОГО СЛОВ О САМОЙ ПРОГРАММЕ.
Программа написана на языке ACL и для ее работы необходим интерпретатор этого языка. Интерпретатор сделан таким образом, что он, фактически, сам является программой vkACL.jar с внешним видом и большими возможностями для работы. Он не только дает средства разработки новых программ на этом языке, но и является удобным навигатором по большому набору уже написанных автором программ. Автором как языка, так и программы интерпретатора, является автор этой статьи. Программа томографии является одной из большого списка таких программ. Программа vkACL.jar имеет свой сайт , откуда ее можно скачать. Там же дается архив большого числа базовых программ, представляющих интерес для всех. Научные программы в этот архив не входят, их можно скачать из депозитария, каждую программу отдельно, кому что надо. В этом депозитарии находится и программа томографии, а которой идет речь в этой статье. Она там называется (Com Tom). Сама программа vkACL.jar написана на языке Java, и для работы также нужен интерпретатор. Но это уже широко известная техника, на сайте написано что и как надо делать.
Описание программы томографии сейчас вынесено в интернет и его можно прочитать даже без скачивания самой программы вот по этой ссылке . Разумеется, что и новая версия еще далека от совершенства и наверняка я буду развивать ее и дальше, просто невозможно делать то, о чем не знаешь. Всегда и любая разработка основана на том опыте. который удается получить. При возникновении новой информации появляются новые идеи и становится ясно, что надо делать дальше. Но надо понимать, что данный вид деятельности весьма широко развит, очень широко используется в медицине и сделано уже много. Однако меня и моих коллег интересуют микро и нано объекты. И вот как раз для них еще сделано не так много. И проблемы существуют.
ДЕМОНСТРАЦИЯ И ТЕСТИРОВАНИЕ ПРОГРАММЫ.
Любая приличная программа должна включать в себя возможность для пользователя посмотреть как она работает (демонстрация), средства тестирования, а также средства максимального легкого обучения работы с ней. Это особенно касается коммерческих программ, но и просто научные программы, авторы которой заботятся о пользователях, тоже должны иметь такой набор. В случае с программой томографии ситуация усложняется тем, что для работы программы надо иметь очень много файлов входных данных, которые к тому же имеют специальный формат tiff, который совсем не так прост, как хотелось бы.
С другой стороны, эти данные получены от неизвестного объекта и сравнить результаты не с чем. Самый лучший способ проверки правильности работы программы -- это использование модельных объектов, структура которых известна с самого начала. Этот способ предполагает на первом этапе расчет "экспериментальных" проекций от известного объекта. Этот этап работы не является необходимым для самого метода. Он является дополнением именно с целью демонстрации, проверки и обучения. Сам метод и основная функция программы используются на втором этапе для восстановления структуры объекта с последующим сравнением результата реконструкции с оригиналом. Именно сравнение с оригиналом позволяет понять насколько правильно работает программа.
Модельным объектом в данной программе является картинка в формате png. Сложность в том, что необходимо очень точно приготовить "экспериментальные" проекции из оригинальной картинки. Если это делать не очень точно, то восстановленная картинка с оригиналом не совпадет. Картинки могут быть разными и их может быть много. В данной статье я рассмотрю только одну, а именно, фантом Shepp-Logan, который часто используется для тестирования алгоритмов КТ. Я впервые увидел его в учебной программе (LearnCT), которую можно найти в интернете по названию. К сожалению многие статьи в интернете меняют адреса. На момент написания данной статьи программу можно скачать тут . Она тоже написана на языке Java.
Программа (LearnCT) даже судя по ее названию является учебной, но у нее версия 0.3 не меняется много лет, нет описания и использовать ее очень сложно. С другой стороны, она использует оба метода, как FBP, так и алгебраический, и не выглядит слишком примитивной. Вероятно есть и другие бесплатные программы, но я пока ничего не нашел. Да особенно и не искал.
Я нашел аналогичную картинку в интернете, привел ее к размеру 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=180, то есть угол менялся с шагом на один градус. Как следует из описания основного расчетного модуля программы, у нее есть 4 способа работы. Номер 1 -- это расчет синограмм, номер 2 -- это расчет томограмм, номер 3 -- вырезания части картинки из целой картинки, номер 4 -- расчет проекций из картинки.
Для работы по номеру 4 надо открыть новую папку внутри папки программы (я открыл папку var000) и записать туда картинку с названием source.png. Она должна быть квадратной и иметь число пикселей в виде целой степени 2. У меня было 512. Это наиболее удобный размер. При этом программа работает относительно быстро, но все хорошо видно на экране. Программа расчета проекций каждую проекцию записывает в отдельный файл формата tiff так, как делают детекторы, то есть она записала 180 файлов с номерами от 000 до 179. На этом этап моделирования экспериментальных данных закончился.
Затем в работу вступили основные способы работы программы томографии, которые и тестируются. Первый способ вычисляет синограммы для всех сечений трехмерного объекта по высоте, если объект имеет высоту в направлении оси вращения. Синограмма -- это матрица размерности D*N, в которой все проекции заданного сечения лежат как бы одна над другой. То есть программа выбирает из каждого файла двумерной проекции (если объект имеет высоту) одну строку на заданной высоте, и записывает ее в строку синограммы с вертикальной координатой, соответствующей углу проекции. Эту матрицу она записывает в файл, но одновременно определяет минимальное и максимальное значения матрицы. Их она записывает в другой текстовый файл с названием minmax.txt.
Эта операция выполняется относительно быстро, и ее необходимо выделить для того, чтобы определить самое максимальное и самое минимальное значения для всех синограммах. На втором этапе расчета это значения надо указать в форме ввода входных данных. Правильное задание этих значений влияет на результат. Но можно указывать не только сами эти значения и и любые другие внутри интервала, который они задают. Нельзя лишь выходить за этот интервал. Программа все числа, которые выходят за указанный интервал ставит на ближайшую границу. То есть меньшие на минимум, большие на максимум. Программа вычисляет логарифм максимального значения и потом она вычитает из этого значения логарифм всех чисел синограммы. Соответственно, для максимального значения разница будет равна нулю.
Итак, на первом этапе программа записывает синограмму, а затем из нее получает томограмму в виде картинки точно такого же размера и типа, как и оригинал. Ниже слева показана картинка оригинала, а справа томограмма. Картинки даются в натуральную величину, то есть одна точка экрана -- один пиксель, если размер показа в браузере выставлен без масштабирования. Время расчета томограммы на моем компьютере 9 сек, но это зависит от параметров компьютера. На слабом компьютере может быть и больше. В новой версии программы нужно задавать сдвиг записи картинки от начала в tiff файле. Раньше предполагалось, что картинка находится в конце файла и сдвиг определялся автоматически. Но это не всегда так бывает. Для tiff файлов, которые делает сама программа, это именно так, и сдвиг равен 168.
Рис. 1. Тест с медицинским фантомом Shepp-Logan для 180 проекций. Слева оригинал, справа томограмма.
Видно, что томограмма почти совпадает с оригиналом, как по форме объектов, так и по цвету. Единственное отличие состоит в том, что правая картинка имеет мелкую структуру, как мелкая рябь на воде при слабом ветре. Причина в том, что 180 проекций не хватает, чтобы очень точно восстановить объект. Расчет с 360 проекциями дает такую же картинку, но рябь становится существенно меньше. Записывая проекции с разным их числом, то есть с разным шагом по углу вращения легко изучить как это влияет на качество томограмм. Я не стану здесь это показывать, так как такой анализ я уже публиковал раньше тут .
ОБУЧЕНИЕ. ПЕРВЫЕ ШАГИ.
Рисунок 1 является одновременно демонстрацией и тестированием. Для обучения необходимо указать первые шаги. В данном случае использовался только основной расчетный модуль программы. Для ввода входных данных он использует вертикальную форму подписанных окон ввода. Всего их 22. Сама программа аккумулирует результаты ввода в текстовой строке, в которой один результат от другого отделен символом вертикальной черты. В файле с названием (сt.par) имеется 26 таких строк. В самой первой строке записываются текущие значения, которые программа показывает, а остальные 25 строк содержат значения, которые запоминаются и потом используются по кнопкам формы (Save var) и (Choose var). Для указанного на рисунке 1 расчета использовались следующие значения
var000|z000.tif|0|1|180|512|1|512|1|0|0|1|0 0|0|0|1|1|32500|25019|168|2|
Папка с названием (var000) в дистрибутиве программы уже присутствует. В этой папке уже есть файл с названием (source.png) и внутрення папка с названием (r). Начиная новый проект папки такого типа надо открывать вручную средствами операционной системы. И туда надо записывать tiff файлы экспериментальных проекций. В демонстрационно-тестовом режиме проекций нет. Чтобы их создать надо сначала в 22-м окне 2 заменить 4. Эта операция использует не все поля ввода, а только 1-е (имя папки), 5-е (число проекций), 8-е (размер проекций), 18-е (степень поглощения) и 19-е (максимальное значение). При этом 18-е поле ввода используется в другом смысле чем в основном расчете. После работы программы в указанной папке появляются 180 файлов проекций. Они имеют ширину 512 и высоту 1.
Теперь все готово для основной работы программы. Для расчета синограммы нужно в 22-м поле ввода поставить значение 1. Все остальные значения уже указаны выше. Тут важно отметить, что 21-е поле требует знать какой сдвиг в байтах картинка имеет относительно начала файла. Это таг с номером 273 в tiff файле. Узнать его значение можно только с помощью специальных программ. Такие программы в интернете есть. Но важно, что такая программа есть среди готовых ACL программ интерпретатора ACL. Для это надо нажать клавиши Alt+B, выбрать кнопку (Img pro), затем (Special operations) и наконец (Find tags in tiff files). Максимальное и минимальное значения на этом этапе задавать не надо, они как раз вычисляются.
После выполнения операции в папке (var000/r) появится два файла (sin0000.dat) и (sin0000.jpg). Второй файл можно сразу посмотреть. Его покажет даже сама операционная система. В программе vkACL тоже есть много способов посмотреть файл. В том числе и среди кнопок главного меню программы. Для этого надо прочитать описание программы на сайте. Это совсем не проблема для тех, кто смотрит фотографии на компьютере. Убедившись, что все в порядке и файлы есть, можно приступать к основной работе.
Просто меняем значение в 22-м окне на 2 и вводим правильные значения для минимального и максимального значений на синограмме. И снова запускаем программу в работу. Это будет уже не так быстро. И в результате в той же папке появится еще два файла. Как раз файл jpg из тех, которые появляются на этом этапе, и показан на рисунке 1. Расчет можно повторить с другим числом прекций, или меняя другие значения параметров. Или с другой картинкой. Этого достаточно для учебы. А как делать более сложные расчеты с реальными экспериментальными данными написано в описании работы программы и будет более детально написано в других документах.
.