СИСТЕМА АВТОМАТИЧЕСКОГО ПОИСКА И КОРРЕКЦИИ ОСИ ВРАЩЕНИЯ ОБЪЕКТА Российский патент 2024 года по МПК G01N23/46 

Изобретение относится к области рентгенотехники и может быть использовано как в медицине, так и в области технического неразрушающего контроля.

Из уровня техники известен проекционный способ позиционирования центра вращения системы компьютерной томографии, раскрытый в патентном документе CN 102539460 А, опубликованном 04.07.2012. Изобретение раскрывает проекционный способ позиционирования центра вращения системы компьютерной томографии. Способ включает этапы: получение данных проекции на 360 градусов для получения проекционной синусоидальной диаграммы; установка надлежащего порога, сегментация проекционной синусоидальной диаграммы, наложение и усреднение данных проекции всех углов проекции в соответствии с соответствующими пикселями, т.е. наложение всех линий проекционной синусоидальной диаграммы и усреднение для получения одномерного сигнала y=f(s), с использованием квадратичного многочлена y=as2+bs+c и подбора методом наименьших квадратов y=f(s); и вычисление значения горизонтальных координат s=b/(-2a) симметричного центра параболы y=as2+bs+c, где значениями горизонтальных координат являются значения координат центра вращения проекции.

Недостатком прототипа является невозможность использования его для конусной модели распространения рентгеновских лучей, представляющей на текущий момент наиболее распространенную модель.

Задачей изобретения является устранение недостатков прототипа.

Технический результат изобретения заключается в повышении качества реконструкции внутренней структуры исследуемого объекта.

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

a) получения проекционных данных,

b) вычисления усредненного проекционного изображения,

c) предобработки усредненного проекционного изображения с помощью линейной нормализации и масштабирования,

d) двухэтапного поиска величин угла наклона и сдвига оси вращения,

e) уточнения положения оси вращения объекта.

Краткое описание чертежей.

На фиг. 1 представлен общий вид системы автоматического поиска и коррекции оси вращения объекта.

На фиг. 1 обозначено: 1 - источник рентгеновского излучения, 2 - изучаемый объект, 3 - устройство поворота и позиционирования изучаемого объекта (пунктирной линией обозначена ось вращения), 4 - плоскопанельный детектор, 5 - вычислительное устройство для сбора и обработки рентгеновских изображений.

На фиг. 2 представлена блок-схема обработки данных.

На фиг. 3 представлена блок-схема двухэтапного определения угла наклона оси вращения.

На фиг. 4 представлена блок-схема двухэтапного определения величины сдвига оси при произведенном повороте усредненного проекционного изображения.

На фиг. 5 представлен график функции оценки несимметрии повернутого предобработанного усредненного проекционного изображения l(Δх) при различных пробных значениях сдвига положения оси Δх для усредненного проекционного изображения фантома.

На фиг. 6 представлен график функции оценки несимметрии повернутого предобработанного усредненного проекционного изображения l(Δх) при различных пробных значениях сдвига положения оси Δх для усредненного проекционного изображения фантома в окрестности точки минимума.

На фиг. 7 представлен график функции оценки несимметрии повернутого предобработанного усредненного проекционного изображения l(Δх) при различных пробных значениях сдвига положения оси Δх для усредненного проекционного изображения тубуса с шахматами.

На фиг. 8 представлен график функции оценки несимметрии повернутого предобработанного усредненного проекционного изображения l(Δх) при различных пробных значениях сдвига положения оси Δх для усредненного проекционного изображения тубуса с шахматами в окрестности точки минимума.

На фиг. 9 представлена блок схема уточнения величины сдвига оси вращения по изначальному усредненному проекционному изображению.

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

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

На фиг. 12 представлен реконструированный центральный слой синтетических данных (фантом с геометрическими фигурами): а.) идеальное изображение, б.) реконструкция центрального слоя без коррекции оси вращения, в.) реконструкция центрального слоя со скорректированной осью вращения.

На фиг. 13 представлен реконструированный центральный слой реальных данных (тубус с шахматами): а.) идеал, реконструкция центрального слоя с подобранными вручную параметрами положения оси вращения, б.) реконструкция центрального слоя без коррекции оси вращения, в.) реконструкция центрального слоя со скорректированной осью вращения.

Схема работы заявленной системы представлена на фиг. 1.

Система состоит из источника рентгеновского излучения 1, плоскопанельного детектора 4, между которыми установлено устройство поворота и позиционирования 3 изучаемого объекта 2 и вычислительного устройства 5, подключенное к детектору 4.

Схема сканирования является круговой. Детектор - плоскопанельный. Геометрическая модель распространения излучения - параллельная, веерная или конусная, размеры детектора - от 100×100 до 5000×5000 пикселей, угловой размер объекта - до 30 градусов. Источник - синхротрон или лабораторный.

В время сканирования данные, полученные плоскопанельным детектором, передаются в вычислительное устройство, в котором осуществляются следующие этапы обработки данных:

- получение проекционных данных,

- вычисление усредненного проекционного изображения,

- предобработка усредненного проекционного изображения с помощью линейной нормализации и масштабирования,

- двухэтапный поиск величин угла наклона и сдвига оси вращения,

- уточнение положения оси вращения объекта.

Получение проекционных данных

На первом этапе вычислительное устройство получает проекционные данные от плоскопанельного детектора.

Вычисление усредненного проекционного изображения.

Вычисление арифметически усредненного проекционного изображения :

где p_ϕ=p_ϕ(х, z) - проекционный снимок объекта (одноканальное изображение), отвечающий углу проецирования ϕ ∈ Ф, Ф - множество обозреваемых углов проецирования в количестве |Ф|.

Предположим, что усредненное проекционное изображение - это осесимметричное изображение, причем ось симметрии совпадает с проекцией оси вращения объекта на плоскость окна детектора при просвечивании источником рентгеновского излучения. В случае наличия сдвигов усредненное проекционное изображение можно считать приближенно осесимметричным изображением: несимметричность усредненного проекционного изображения становится заметной при сдвигах, превышающих нескольких десятков-сотен пикселей детектора. В параллельной геометрической схеме при рассмотрении множества допустимых движений оси вращения усредненного проекционного изображения всегда строго осесимметрична.

Усредненное проекционное изображение складывается из множества эллипсов - проекций круговых траекторий движения отдельных материальных точек исследуемого объекта при вращении вокруг оси. Последнее служит обоснованием наличия осевой симметрии усредненного проекционного изображения.

В предположении осесимметричности ось симметрии усредненного проекционного изображения совпадает с проекцией оси вращения объекта на плоскость окна детектора: величина угла наклона оси в плоскости, параллельной плоскости окна детектора, равна величине угла между осью симметрии усредненного проекционного изображения и вертикалью, а сдвиг оси вращения отличается от сдвига оси симметрии усредненного проекционного изображения от центра изображения в количество раз, задаваемое коэффициентом увеличения конусной геометрической схемы, или совпадает со сдвигом оси симметрии усредненного проекционного изображения в параллельной геометрической схеме. Таким образом, задача определения скалярных параметров положения оси вращения объекта -величин сдвига и наклона в плоскости, параллельной плоскости окна детектора, - сведена к задаче нахождения уравнения оси симметрии усредненного проекционного изображения, в частности, ее угла наклона и сдвига в горизонтальном направлении (от центра изображения).

Предобработка усредненного проекционного изображения с помощью линейной нормализации и масштабирования.

После вычисления усредненного проекционного изображения производится его предобработка. Этап предобработки изображения включает линейную нормализацию и масштабирование.

Сначала осуществляется линейная минимаксная нормализация усредненного проекционного изображения, при которой значения интенсивности пикселей усредненного проекционного изображения приводятся к шкале [0,1].

Затем усредненное проекционное изображение масштабируется в горизонтальном и вертикальном направлениях с одним и тем же коэффициентом масштабирования При масштабировании размерности усредненного проекционного изображения уменьшаются так, что их отношение остается постоянным, причем наименьшая размерность масштабированного (с усреднением) изображения составляет 128 пикселей. В рамках шага масштабирования усредненного проекционного изображения производится интерполяция.

Масштабирование применяется с целью ускорения последующих шагов обработки данных, включающих перебор различных значений сдвига оси вращения - по всей ширине изображения.

Будем различать предобработанное и непредобработанное усредненное проекционное изображение и использовать различные обозначения для этих объектов: исходное непредобработанное изображение и предобработанное (нормализованное и масштабированное) изображение

Двухэтапный поиск величин угла наклона и сдвига оси вращения

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

Этап 1 заключается в поиске угла наклона оси вращения по разреженной сетке значений, в рамках этапа 2 найденная на этапе 1 величина угла наклона оси уточняется по сетке с уменьшенным шагом.

Опишем детально каждый из двух этапов текущего шага метода. Этапы визуализированы посредством блок-схемы, представленной на фиг. 3.

Этап 1: поиск значения угла наклона оси вращения по разреженной сетке значений.

Поворот предобработанного усредненного проекционного изображения. Имеем предобработанное усредненное проекционное изображение Будем искать первое приближение к искомой величине угла наклона оси вращения в промежутке от start_angle до end_angle градусов с шагом first_step.

Параметры first_step, start_angle и end_angle разреженной равномерной сетки рассматриваемых значений могут быть выбраны равными 1.0, -15.0 и 15.0. Для поиска первого приближения величины угла наклона оси шаг сетки выбирается не превосходящим 1 градуса.

Для каждого значения angle из рассматриваемого промежутка производится поворот предобработанного изображения на угол angle градусов. Поворот осуществляется относительно центра изображения с использованием интерполяции.

Повернутые на angle градусов предобработанное изображение и исходное усредненное проекционное изображение будем обозначать и соответственно. В таких обозначениях, в частности,

Заметим, что в силу симметрии относительно проекции оси вращения на плоскость окна детектора, если angle - правильно найденная величина угла наклона оси, имеет в точности вертикальную ось симметрии. Следовательно, наилучшее приближение angle истинной величины угла наклона оси соответствует той ситуации, когда определяемая далее функция оценки вертикальной несимметрии изображения минимальна.

Целевой функционал. Функция оценки вертикальной несимметрии усредненного проекционного изображения. Для измерения функции оценки вертикальной несимметрии произвольного изображения Р определим l(Δх) - функцию оценки несимметрии усредненного проекционного изображения Р относительно вертикали, сдвинутой на Δх пикселей относительно серединной вертикали изображения в горизонтальном направлении (целевой функционал):

где F - оператор зеркального отражения изображения в горизонтальном направлении, S_Δx - оператор сдвига изображения на Δх в горизонтальном направлении, ° - операция композиции отображений, l_q - стандартная q норма Минковского. Вместо нормы Минковского может быть вычислена иная матричная норма.

Вертикальная ось, сдвинутая горизонтально относительно центра изображения на Δх пикселей, является осью симметрии изображения Р, если при горизонтальном отзеркаливании Р и последующем сдвиге на 2⋅Δх пикселей также получим исходное изображение Р. Функция оценки l(Δх) несимметрии изображения относительно оси, сдвинутой от центра на Δx пикселей, есть l_q - норма разности между исходным изображением Р и его зеркально отраженной и сдвинутой на 2⋅Δх пикселей копией (S_2⋅Δх°F)[P]. Чем меньше значение l(Δх), тем более уверенными мы можем быть в том, что х=Δх - вертикальная ось симметрии изображения (отсчитывая от центра изображения).

Функция оценки вертикальной несимметрии l_P изображения Р показывает, насколько вероятно наличие у изображения Р какой-либо вертикальной оси симметрии:

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

Значение l_P=0 (отметим, l_P≥0) функции оценки вертикальной несимметрии говорит о существовании у усредненного проекционного изображения Р вертикальной оси симметрии. Чем больше значение 1_Р, тем менее симметрично изображение (в предположении - достаточно мелкая сетка).

Вертикальной осью симметрии усредненного проекционного изображения Р мы считаем такую вертикальную сдвинутую относительно центра изображения на Δх пикселей ось, на которой функционал l_P(Δх) достигает значения l_P:l_P(Δх)=l_P - условие того, что рассматриваемая вертикальная ось является вертикальной осью симметрии усредненного проекционного изображения Р.

Двухэтапная минимизация целевого функционала. Искомая величина угла наклона оси вращения такова, что повернутое предобработанное усредненное проекционное изображение имеет вертикальную ось симметрии, то есть оптимальное значение angle подбирается из условия минимизации функции оценки вертикальной несимметрии усредненного проекционного изображения

Для каждого значения angle разреженной сетки (от start_angle до end_angle градусов с шагом first_step) значение вычисляется в ходе двухэтапной минимизации целевого функционала (Δх) (см. фиг. 4).

Этап 1: минимизация целевого функционала по разреженной сетке значений Δх от -width до width пикселей с большим шагом first_step, где width - ширина изображения , а шаг first step может быть принят равным 5 (не менее нескольких пикселей). В конце этапа 1 получаем начальное приближение минимального значения функции оценки вертикальной несимметрии и значения Δх, задающего положение минимизирующей вертикальной оси симметрии. Графики целевого функционала (Δх) (при величине угла angle, совпадающей с истинной величиной угла наклона оси вращения) приведены на фиг. 5, 6, 7, 8.

Этап 2: минимизация целевого функционала по измельченной сетке значений сдвига Δх, в ходе которой последовательно для каждого p ∈ {1,0, -1} производится минимизация целевого функционала (Δх) в окрестности ранее найденной точки минимума Δх=best_shift радиуса 10^р с шагом 10^р-1 (шаг равен линейному размеру пикселя при р=1, десятой и сотой доле линейного размера пикселя при р=0 и р=-1 соответственно, минимизация прямым перебором). Примеры графиков целевого функционала в окрестности точки минимума можно увидеть на фиг. 6 и 8.

Результат этапа 2 - вычисленное значение и Δх, при котором целевой функционал (Δх) равен (при предварительно произведенном повороте предобработанного изображения на угол величины angle относительно центра усредненного проекционного изображения).

Далее предобработанное усредненное проекционное изображение поворачиваем на другой угол (см. фиг. 3), вычисляем для произведенного угла поворота angle функцию оценки вертикальной несимметрии Если найденное значение norm функции оценки вертикальной несимметрии повернутого усредненного проекционного изображения меньше функции оценки при предыдущем повороте (то есть min_norm > norm), обновляем значения угла наклона и сдвига оси angle_best=angle, shift_best=shift и min_norm=norm. Перебираются величины угла angle поворота от start_angle до end_angle градусов с шагом first_step.

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

Этап 2: уточнение величины наклона оси вращения.

В рамках этапа 2 шага вычисления значения угла наклона оси вращения рассматриваются измельченные сетки значений угла наклона в окрестности найденного во время исполнения этапа 1 оптимального значения angle_best. По измельченным сеткам производится минимизация функции оценки вертикальной несимметрии повернутого предобработанного усредненного проекционного изображения

Опишем этап 2 более подробно. Для каждого р ∈ {0, -1} рассматривается сетка значений angle от angle best -2⋅10^р до angle_best +2⋅10^р с шагом 10^р-1 (десятая доля градуса при р=0 и сотая доля градуса при р=-1), где angle_best - ранее определенное оптимальное значение угла наклона оси. Для каждого angle производится вычисление функции оценки вертикальной несимметрии посредством уже описанной двухэтапной минимизации целевого функционала. Если новое значение функции оценки вертикальной несимметрии меньше ранее полученных при других значениях angle, обновляем значения angle_best=angle, shift_best=shift и min_norm=norm.

Результат этапа 2 - уточненное по измельченным сеткам оптимальное значение угла наклона оси вращения angle_best. Для вычисленного значения angle_best также находится значение сдвига оси вращения shift_best, при котором найденное минимальное значение функции оценки несимметрии усредненного проекционного изображения достигается. Значение angle_best является ответом на вопрос о значении угла наклона оси вращения и далее не уточняется.

Значение shift_best служит начальным приближением искомой величины сдвига оси вращения, которое далее подлежит уточнению.

Уточнение величины сдвига оси вращения с помощью усредненного проекционного изображения.

На предыдущем шаге были вычислены значения угла наклона оси вращения и сдвига оси вращения. На шаге уточнения величины сдвига оси вращения ранее найденное значение shift_best уточняется.

Уточнение значения сдвига оси производится при помощи исходного непредобработанного усредненного проекционного изображения (см. фиг. 9).

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

1) Производится линейная минимаксная нормализация изображения (значения интенсивности пикселей усредненного проекционного изображения приводятся к шкале [0,1]). Для единообразия результат нормализации исходного усредненного проекционного изображения обозначим Подчеркнем, что усредненное проекционное изображение не масштабируется.

2) Нормализованное усредненное проекционное изображение поворачивается вокруг центра изображения на найденную величину угла наклона оси вращения, равную angle_best. Результат поворота будем обозначать Повернутое нормализованное усредненное проекционное изображение обладает вертикальной осью симметрии, которую будем искать в окрестности положения вертикальной оси Поскольку начальное приближение значения сдвига оси было получено с помощью масштабированного усредненного проекционного изображения, shift_best делится на масштабирующий множитель scaling_factor.

3) Рассмотрим сетку возможных значений сдвига оси вращения от до пикселей с маленьким шагом, равным сотой доле линейного размера пикселя. По этой сетке прямым перебором значений сдвигов Δх минимизируется значение целевого функционала Точка минимума shift_best - итоговый ответ метода на вопрос о значении сдвига оси вращения. Примеры графиков целевого функционала, значения которого вычислены для различных значений сдвига оси вращения Δх по нормализованному немасштабированному изображению приведены на фиг. 10 и 11.

После вычисления shift_best следует учесть специфику модели распространения рентгеновских лучей, а именно в конусной модели необходимо разделить shift_best на коэффициент увеличения конуса (то же касается веерной схемы). Параллельная модель распространения рентгеновских лучей не требует нормировки найденной по усредненному проекционному изображению величины сдвига оси вращения объекта.

Итак, значения угла наклона и сдвига оси вращения найдены.

Коррекция положения оси вращения

После нахождения параметров положения оси вращения (величин угла наклона и сдвига оси вращения) каждый проекционный снимок р_ф, ϕ ∈ Ф, поворачивается относительно центра изображения на величину угла наклона оси вращения angle_best. Затем производится сдвиг каждого повернутого проекционного снимка на shift_best пикселей. При вычислении поворотов и сдвигов проекционных снимков используется интерполяция.

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

Пример осуществления представленного способа

Результаты работы изобретения на синтетических и реальных данных представлены на фиг. 12 и 13. На рисунках изображены реконструкции центральных томографических слоев синтетических данных (фантом с геометрическими фигурами) и реальных данных (тубус с шахматами) с коррекцией и без коррекции оси вращения предложенным автоматическим способом выравнивания оси вращения (фиг. 12б, 12в и 13б, 13в). На фиг. 12а, 13а представлены идеальные изображения центрального слоя фантома и тубуса с шахматами, с которыми производится сравнение полученных реконструированных центральных слоев. Для измерения изменения качества полученных реконструкций проведено вычисление L₂ нормы разности реконструкции и идеального изображения: если реконструированное изображение представляется матрицей а идеальное изображение - матрицей то значение l₂ нормы разности изображений равно Для реконструкции использовался алгоритм Feldkamp-Davis-Kress. Значение l₂ нормы разности реконструкций и идеальных изображений оформлены в виде таблицы 1.

Синтетические проекционные снимки (пакет 1024×1024×720 вокселей) симулированы с углом наклона оси вращения величины 9.20 градуса в плоскости, параллельной плоскости окна детектора, и сдвигом величины 2.83 пикселя. Найденное предложенным методом значение угла наклона оси составляет 9.20 градуса, найденное значение сдвига оси равно 2.81 пикселя. Результат коррекции положения оси вращения - на рис. 12в. Реконструкция центрального слоя при коррекции оси вращения визуально схожа с идеальным изображением центрального слоя фантома (ср. 12а и 12в). Меньшее значение L₂ нормы разности реконструкции и идеального изображения в случае использования коррекции оси вращения подтверждает улучшение качества реконструкции при использовании предложенного средства коррекции оси (таблица 1, ср. фиг. 12б и 12в).

Применительно к пакету реальных данных тубуса с шахматами (пакет 2064×1548×1024 вокселей) визуально очевидно улучшение качества реконструкции (рис. 13б и 13в). Найденная величина угла наклона оси вращения равна 0.33 градуса, найденная величина сдвига оси составляет -7.09. Ручной подбор параметров положения оси вращения подтверждает точность результатов метода: оптимальное значение сдвига между -7.00 и -8.00 пикселями, в то время как оптимальный угол наклона - между 0.00 и 0.50 градусов. На фиг. 13а представлена реконструкция центрального слоя с оптимальными параметрами оси, подобранными вручную и равными -7.00 пикселей и 0.30 градусов.

Приведенные примеры показывают, что использование описанного средства автоматического определения положения оси вращения позволяют получать томографические реконструкции более близкие идеальному изображению как при визуальном восприятии, так и по измеренной L₂ норме (см. табл. 1)., чем в случае без коррекции оси.

Использование предлагаемого изобретения предваряет шаг реконструкции внутренней структуры объекта. Коррекция положения оси вращения перед последующей реконструкцией дает возможность получить реконструированное изображение более высокого качества.

При этом описанное изобретение может быть использовано и в качестве средства предварительной калибровки измерительного оборудования при помощи специальных фантомов. Такое применение изобретения позволяет улучшить качество результатов в последующих целевых измерениях.

Использование: для автоматического поиска и коррекции оси вращения объекта при выполнении рентгеновской томографии. Сущность изобретения заключается в том, что система автоматического поиска и коррекции оси вращения объекта при выполнении рентгеновской томографии содержит источник рентгеновского излучения, устройство поворота и позиционирования объекта, плоскопанельный детектор и вычислительное устройство, при этом плоскопанельный детектор выполнен с возможностью формирования проекционных данных объекта, а вычислительное устройство выполнено с возможностью: а) получения проекционных данных, b) вычисления усредненного проекционного изображения, с) предобработки усредненного проекционного изображения с помощью линейной нормализации и масштабирования, d) двухэтапного поиска величин угла наклона и сдвига оси вращения, е) уточнения положения оси вращения объекта. Технический результат: обеспечение возможности повышения качества реконструкции внутренней структуры исследуемого объекта. 13 ил., 1 табл.

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

a) получения проекционных данных,

b) вычисления усредненного проекционного изображения,

c) предобработки усредненного проекционного изображения с помощью линейной нормализации и масштабирования,

d) двухэтапного поиска величин угла наклона и сдвига оси вращения,

e) уточнения положения оси вращения объекта.