Цифровая обработка изображений (стр. 5 )

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

,

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

. (5.27)

Соответственно, интерполированный элемент изображения вычисляется как

(5.28)

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

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

Рис. 5.12. Результат компенсации дробного сдвига с применением линейного прогноза

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

и (5.29)

, ,

где - область, в которой изображения можно считать стационарными, а - количество точек в этой области.

Рис. 5.13 иллюстрирует эффективность применения линейного прогноза в задаче обнаружения малоразмерных объектов на медленно движущемся фоне. На рис. 5.13а и б показаны два изображения, отличающиеся небольшим сдвигом фона и тем, что на рис. 5.13б присутствует четыре одинаковых по форме, но разных по амплитуде объекта. Оба изображения сопровождаются некоррелированным случайным шумом, среднеквадратичное отклонение (СКО) которого существенно меньше СКО фона, но сравнимо с амплитудой объектов. Амплитуда верхнего правого объекта составляет 1.5 СКО шума, верхнего правого – 3.0, нижнего левого – 4.5, и нижнего правого – 6.0. Для обнаружения применяется изложенная в предыдущем разделе согласованная фильтрация, причем пороговое значение для принятия решения о наличии объекта выбирается так, чтобы ложных тревог не возникало. Предварительно перед фильтрацией осуществляется подавление фона путем вычитания первого изображения из второго.

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

На рис. 5.13 в, г и д показаны объекты, обнаруженные соответственно на простой разности, на разности, вычисленной после компенсации целочисленного сдвига, и на разности, полученной посредством линейного прогноза. Последний способ, как и следовало ожидать, дает наилучший результат.

Вопросы

Как можно представить искажения, вносимые в изображения линейной системой? В чем состоит задача восстановления изображения, искаженного линейной системой? Чем отличается фильтр Винера от инверсного фильтра? В чем проявляется при восстановлении конечность размеров изображения? Как можно ослабить влияние конечности размеров на результат восстановления? В чем заключается задача обнаружения объектов на изображении? На чем основан выбор порогового значения при обнаружении? Чем обусловлено повышение качества обнаружения при использовании согласованного фильтра? Для чего необходимо совмещать изображения? Что дает устранение дробного сдвига между дискретными изображениями? Какая модель сигнала лежит в основе метода линейного предсказания?

Литература

, и др. Цифровая обработка изображений в информационных системах. Учебное пособие. Новосибирск, НГТУ, 2002 Ван Теория обнаружения, оценок и модуляции, т.1., М., “Советское радио”, 1972

6. Компьютерная томография

Открытие в 1895 проникающего излучения, впоследствии названного по имени человека, его обнаружившего, рентгеновским, впервые позволило заглянуть внутрь непрозрачных в видимом свете объектов. Самое широкое практическое применение рентгеновское излучение нашло в медицине для получения изображений костных тканей и внутренних органов человека. В дальнейшем рентгеновская технология исследования внутренней структуры непрозрачных в видимом свете объектов распространилась в область материаловедения, диагностики промышленных изделий, строительных конструкций и т. д. Как и любое другое, рентгеновское излучение ослабляется средой, через которую оно проходит, и степень ослабления зависит от физических свойств среды и от длины пройденного в среде пути, что и является основой получения рентгеновских изображений. Важной его особенностью является способность распространяться в среде практически без рассеяния и преломления, обеспечивающая при использовании точечного источника излучения получение “теневых” изображений. При прохождении излучения через однородную поглощающую среду без рассеяния и преломления интенсивность его изменяется в соответствии с законом Бугера-Ламберта [6.1]

, (6.1)

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

,

где - координата вдоль направления распространения излучения, и - координаты входа в среду и выхода из нее. Полагая при , , это соотношение можно переписать в виде

. (6.2)

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

Соотношения (6.1) и (6.2), в частности, показывают, что четкие теневые изображения, отражающие структуру объекта, получаются лишь в том случае, когда поглощение среды неизменно вдоль направления распространения излучения, а ширина ее слоя постоянна. В противном случае изменения поглощения в изображении “смазываются”.

Преобразование Радона. Задача томографии

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

, (6.3)

где - угол между направлением оси OX и нормалью к лучу (рис. 6.1).

Рис. 6.1. Схема построения проекции в плоском параллельном пучке

Тогда интенсивность излучения на выходе из объекта, в соответствии с (6.2) может быть представлена в виде

, (6.4)

где - дельта-функция.

Интеграл

, (6.5)

стоящий в показателе экспоненты выражения (6.4), называется проекцией объекта в направлении . Делая замену , в выражении (6.5), нетрудно убедиться, что . Кроме того, очевидно, что периодична по с периодом .

Отметим, что здесь просто фиксирует плоскость, в которой происходит построение проекции, и фактически является параметром, определяющим вид распределения показателя поглощения в выбранной плоскости. Т. о., формально выражение (6.5) определяет преобразование исходной двумерной функции в двумерную функцию путем интегрирования вдоль секущих, задаваемых уравнением (6.3). Это преобразование в литературе известно как прямое преобразование Радона, а функция называется радоновским образом функции .

Нетрудно вычислить радоновский образ от точечного объекта, заданного функцией

:

.

Следовательно, точечный объект в плоскости () превращается в плоскости () в синусоиду с фазой, определяемой координатами .

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

.

Изображение этой функции приведено на рис. 6.2а. Здесь , , , , , , , . Ее радоновский образ, представленный на рис. 6.2б, имеет вид

.

На рис. 6.2б горизонтальная ось соответствует координате (), а вертикальная – координате ().

а	б
Рис. 6.2. Функция (а) и ее радоновский образ (б)

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

Теорема о центральном сечении

Вычислим одномерный спектр Фурье от проекции по координате :

(6.6)

Как можно заметить, представляет собой двумерный спектр Фурье от функции , вычисленный вдоль прямой , , проходящей через начало координат в частотной плоскости (). Тем самым доказана теорема, известная как теорема о центральном сечении [6.2]:

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

С другой стороны, функция выражается через свой двумерный спектр как

.

Вводя в частотной области полярную систему координат (), связанную с системой координат () соотношениями

, ,

и делая замену переменных в предыдущем выражении, его можно представить в виде

,

или, учитывая (6.6),

. (6.7а)

Из того, что , следует, что . Учитывая это, получим окончательно

. (6.7б)

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

Фурье-алгоритм восстановления томограммы

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

Рис. 6.3. Одномерные Фурье-образы проекций в двумерном частотном пространстве

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

Для обеспечения хорошего качества интерполяции необходимо, чтобы расстояния между соседними лучами в полярной решетке на всех частотах не превышало шага декартовой решетки. Следовательно, должно выполняться условие , где - количество проекций (направлений пучка, при которых регистрируются проекции), - шаг декартовой решетки по , . Из этого вытекает требование к количеству проекций: .

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

. (6.8а)

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

. (6.8б)

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

Восстановление томограммы методом свертки и обратного проецирования

Подставим в (6.7б) вместо его представление через проекцию :

.

Изменяя порядок интегрирования, получим

(6.9)

где введено обозначение

. (6.10)

Выражение (6.9) определяет способ восстановления непосредственно через проекции без вычисления спектров. Если внимательно посмотреть на внутренний интеграл в (6.9), то можно обнаружить, что он имеет вид одномерной свертки:

,

причем зависит как от угла, под которым получена проекция, так и от координат восстанавливаемой точки: . Здесь подстрочный индекс показывает, что результат свертки относится к проекции, зарегистрированной под углом . Следовательно, является результатом фильтрации проекции одномерным линейным пространственно инвариантным фильтром с импульсным откликом . Будем называть фильтрованной проекцией.

Представим (6.9) в виде

.

Чтобы лучше понять смысл этого выражения, вычислим приближенное значение интеграла, заменив его суммой

, , .

Обозначим вклад, вносимый в сумму проекцией, полученной под углом , как , и вычислим вдоль прямой , направление которой совпадает с направлением просвечивания при построении выбранной проекции. Вдоль этой прямой , т. е. во все точки, лежащие на ней, данная проекция вносит одинаковый вклад. Она как бы “размазывается” по плоскости () в направлении просвечивания. Восстановление происходит за счет суммирования всех “размазанных” проекций. Такая процедура называется обратным проецированием. На рис. 6.4 показаны две “размазанных” фильтрованных проекции объекта, приведенного на рис. 6.2а (а – полученная при , б – при ) и результат их накопления (в).

а	б	в
Рис. 6.4. Восстановление посредством обратного проецирования фильтрованных проекций

Определенные трудности возникают с фильтрацией проекций. Дело в том, что импульсный отклик фильтра (6.10) определен через расходящийся интеграл. Обойти эту трудность при реализации фильтра можно, если учесть, что спектр Фурье проекций ограничен частотой . Тогда выражение (6.7б) можно модифицировать следующим образом