CC BY
25
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2012 Управление, вычислительная техника и информатика № 4(21)
УДК 532.57, 519.6
В.А. Дедок, М.П. Токарев, А.Н. Бондаренко
ЗАДАЧА ВОССТАНОВЛЕНИЯ ПОЛОЖЕНИЙ ДИСКРЕТНЫХ РАССЕИВАЮЩИХ ЦЕНТРОВ ПО НАБОРУ ПРОЕКЦИЙ1
Рассматривается актуальная задача определения максимального подмножества множества точек в Я3по ограниченному набору их двумерных проекций. Данная задача естественно возникает в приложениях физической гидроаэродинамики для оптической диагностики реальных течений жидкости и газа путем измерения мгновенных полей скорости в объеме потока. Предложены методы восстановления указанного множества и определения достаточности измерений для однозначного решения обратной задачи.
Ключевые слова: анемометрия по изображениям частиц, оптическая томография, обратные задачи.
В настоящее время томографические методы широко применяются для исследования внутренней структуры объектов. В основном они используются в медицинской диагностике, в промышленности и строительстве для контроля качества всевозможных изделий и инженерных сооружений. Математический аппарат вычислительной томографии изложен в монографиях [1, 2]. В последнее время все чаще стали появляться примеры использования томографии и в других областях науки и техники. Например, оптическая томография рассеивающих сред сравнительно недавно стала применяться в хорошо зарекомендовавшем себя, методе цифровой трассерной визуализации - Particle Image Velocimetry (PIV) [3].
Благодаря возможности регистрировать мгновенные пространственные распределения скорости метод цифровой трассерной визуализации используется в фундаментальных научных исследованиях для изучения потоков, содержащих вихревые структуры, информация о которых частично теряется при применении одноточечных методов диагностики. К подобным течениям можно отнести турбулентные течения, включая струи, следы, слои смешения и многое другое. Применение полевых методов дает возможность получения информации о динамике структур, их масштабах, расчета дифференциальных характеристик, пространственных и пространственно-временных корреляций, а также статистических характеристик потока. Кроме того, измерительные системы на основе данного метода применяются в прикладных исследованиях для авиастроительной индустрии, автомобилестроения, энергетики, химической промышленности, биологии и медицины.
Основным отличием томографической цифровой трассерной визуализации (Tomo PIV) для экспериментальной гидроаэродинамики от медицинских томографических исследований является малое количество проекций, используемое для восстановления объемного изображения измерительной области потока. Вследствие записи небольшого числа проекций задача получения приближения
1 Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00105), междисциплинарным
интеграционным проектом СО РАН № 14 «Обратные задачи и их приложения: теория, алгоритмы,
программы».
исходного распределения интенсивности рассеянного света на трассерах усложняется, однако форма и разреженность полезного сигнала (локализованные в пространстве узкие пики интенсивности), в конечном итоге, позволяют реконструировать объемное распределение интенсивности требуемого качества, подходящее для дальнейшего корреляционного анализа и оценки смещения трассеров в потоке. В данном случае для томографической реконструкции измерительной области потока применяются итерационные алгебраические методы [4, 5]. На рис. 1 представлена схема расчета поля скорости в измерительном объеме.
Альтернативными методами получения трехкомпонентной трехмерной картины поля течения являются голографический, сканирующий и 3D Particle Tracking Velocimetry (3D PTV) [6]. Последний метод, 3D PTV, по принципу действия наиболее близок к томографической цифровой трассерной визуализации. Однако он требует поддержания сравнительно низкой концентрации трассеров для успешной идентификации отдельных частиц на всех проекциях путем триангуляции и отслеживания их перемещения для оценки скорости.
Перечислим некоторые известные сложности Tomo PIV-измерений. Во-первых, концентрация частиц в области измерения должна быть достаточной для получения хорошего пространственного разрешения оценки скорости. C другой стороны, она не должна превышать критического значения ~0,2 частиц/пикс, выше которого томографическая реконструкция существенно затрудняется вследствие затенения частицами друг друга. При этом стандартные алгоритмы реконструкции ведут к появлению ложных образов трассеров в восстанавливаемом объеме вследствие возрастания неопределенности при недостатке проекционных данных [7, 8]. Таким образом, в эксперименте необходимо поддерживать определенную концентрацию трассеров на проекции, которая бы позволила увеличить вероятность получения истинной конфигурации рассеивающих центров и, в итоге, надежной оценки поля скорости потока.
В данной работе предложен и математически обоснован алгоритм получения максимального множества точек, проекции которого в плоскости, заданной светочувствительными матрицами камер, совпадают с зарегистрированными проекциями. Обсуждаются вопросы необходимости и достаточности проекционных данных для надежного определения внутренней структуры среды, состоящей из набора дискретных рассеивающих центров.
трассерные проекции объемное изображение для расчетной области
Рис. 1. Иллюстрация принципов работы алгоритма оценки поля скорости в методе Tomo PIV
1. Метод реконструкции
Сформулируем задачу в ее математической постановке.
Пусть X = р,х2,...,хп} - множество попарно различных точек в Я3. Обозначим через Рп={р1,Р2,...,рт} - проекцию точек из множества Xна плоскость п. Очевидным образом выполняется неравенство т < п и, в общем случае, т = т(Х,п). Здесь и далее, будем считать, что любая точка из исходного множества не совпадает с любой ее проекцией.
Пусть далее задан набор плоскостей п1,п2,...,пк и набор проекций точек из множества X на каждую из плоскостей Р Рп ,■■■,РПк. Сформулируем следующую обратную задачу: для заданных множеств точек проекций Р Р ,...,РПк на плоскости п1, п2,..., пк восстановить исходное множество X = р, х2,..., хп}.
Очевидным образом, исходное множество точек X для конкретного набора плоскостей и проекций не всегда может быть восстановимо. В частности, задача не имеет однозначного решения, если количество плоскостей равно 2. Без потери общности можно считать, что плоскости ппопарно непараллельны, так как в
случае параллельных плоскостей множества проекций совпадают с точностью до параллельного переноса.
Рассмотрим плоскость пи проекции исходного множества точек на нее.
Тогда
' Х - ҐіЛп1 = С1 і,1 Р},1 + С1і,2Рі,2 + . + С1 і^Ріті,
_ Х2 - ^і,2пі = С2і ,1 Рі,1 + С2і,2 Рі,2 + . + С2т Ріті , (1)
, Хп - іппі = Спі,1 Рі ,1 + Спі,2 Р і,2 + ••• + Спі,т]Рі,т] ■
Здесь пі - вектор нормали к плоскости пі , Р^к - проекция некоторой точки х{
исходного множества X на плоскость п і. Неизвестные коэффициенты ї^п є Я , а Су к таковы, что для каждой строки I все, кроме одного, равны 0, а один равен 1.
Система (1) представляет собой однородную систему линейных уравнений из 3п уравнений и 3(2п + ті) неизвестных х, ґ, с (каждый со своим набором индексов).
Записав систему (1) к раз (по числу проекций), получим систему уравнений (2), решения которой будут определять решение обратной задачи по восстановлению множества X через набор проекций на плоскости п1,п2,...,пк. Система (2) совместна по построению, т.е. решение существует, но, в общем случае, не единственно.
Как показано на рис. 2, даже в двумерном случае (аналогичные примеры можно построить и для Я3) без обладания какой-либо дополнительной информацией об исходном множестве X достаточно сложно, да и не всегда возможно, однозначно восстановить множество по его проекциям.
Добавление третьей плоскости, на которую выполняется проекция исходного множества, расположенной под углом 45° к каждой плоскости (как показано на рис. 3), не уменьшает неопределенности - допустимы 2 конфигурации точек ис-
ходного множества. Изменив угол (сделав его отличным от 45 и 90°) между плоскостями, получим возможность однозначного восстановления множества по его проекциям.
Р 4
Рз
Х2 Хз
Хз
Рі
Х1
Р 2
Рис. 2. Набор проекций р1,р2,р3,р4 на 2 координатные плоскости, определяющий 2 набора исходных точек х1, х2 , х3 и Х1 , Х2 ,
Рис. 3. Сохранение неопределенности
Таким образом, мы доказали следующую теорему:
Теорема 1. Пусть X = р, х2, ..., хп } - множество точек в Я3, Р , Рп,..., Р -
проекции X на набор плоскостей п1,п2,.,(Рп={р1,р2,...,рт}). Пусть также
система (2) имеет единственное решение, такое, что у є Я , все, кроме одного,
Су к для каждой строки I равны 0, а один равен 1. Тогда исходное множество X
единственным образом определяется по множеству проекций.
Замечание 1. Отметим, что теорема 1 явным образом не дает конструктивного алгоритма для вычисления исходного множества. Кроме того, априори требуется знание количества точек в исходном множестве.
Практически более интересной является задача определения характеристик исходного множества без знания какой-либо дополнительной информации о нем.
Как было уже замечено, обратная задача по восстановлению множества по набору проекций в общем случае не имеет единственного решения, поэтому переформулируем ее в несколько ином контексте.
Определение. В условиях предыдущих обозначений назовем множество У = {р, у2,..., р} подходящим, если проекция этого множества на плоскость пу
совпадает с множеством Ру для любого у и любое расширение множества У нарушает это условие.
Тем самым множество У является максимальным множеством, проекции которого на заданные плоскости пу совпадают с множествами Ру. Очевидно, что
X с У . Дадим конструктивный алгоритм восстановления подходящего множества У.
Рассмотрим xq - некоторую точку из исходного множества X. Очевидным образом найдутся минимум 2 точки проекции p, и pj такие, что
xq = p, + nlt1 = p} + П2^ (2)
где n1 и n2 - нормали к плоскостям п1 и п2; t1 и t2 - некоторые неизвестные.
Запишем теперь эти уравнения на t1 и t2 вида (2) для всех возможных p, и Pj. Решения этих уравнений определят множество Z = {p, + n1t} (t - всевозможные решения вышеописанных уравнений). Верно следующее включение X с Y с Z .
Для получения множества Y из множества Z необходимо исключить те точки, которые не являются проекциями на все оставшиеся плоскости п3,..., nk .
Чтобы упростить задачу исключения точек, попадающих в Z \ Y, достаточно применить аффинное преобразование, переводящее плоскости п1, п2 и nq
(3 < q < k) во взаимно ортогональные (координатные) плоскости. Тогда точка (x, у, z) будет принадлежать множеству Y, если точка (х, у, 0) принадлежит множеству проекций на п1, точка (х,0, z) - проекций на п2, а точка (0, у, z) - проекций на nq. В противном случае, ее нужно будет исключить из множества Z и продолжить ту же схему до перебора всех возможных вариантов.
Замечание 1. Для построения исходного множества Z также можно использовать аффинное преобразование, приводящее плоскости п1 и п2 во взаимно ортогональные.
Замечание 2. Последовательное исключение лишних точек из Z обеспечит большую скорость сходимости алгоритма, чем полный перебор. Компьютерная реализация алгоритма поиска также может быть хорошо оптимизирована.
Тем самым мы предложили алгоритм восстановления подходящего множества для множества проекций. Более точное приближение множества Y к множеству X можно выполнять за счет увеличения количества проекций, что очень трудоемко в условиях реальных практических измерений.
Заключение
Рассмотрена обратная задача восстановления максимального множества точек
в R по набору его двумерных проекций. Доказана теорема, позволяющая дать ответ на вопрос о достаточности изменений для однозначного определения множества по набору его двумерных проекций. Предлагается итеративный алгоритм для восстановления подходящего (максимального) множества по его набору проекций.
ЛИТЕРАТУРА
1. Хермен Г. Т. Восстановление изображений по проекциям. Основы реконструктивной томографии. М.: Мир, 1983.
2. Лаврентьев М.М., Зеркаль С.М., Трофимов О.Е. Численное моделирование в томографии и условно-корректные задачи. Новосибирск: Изд-во ИДМИ НГУ, 1999.
3. Raffel M., Willet C.E., Wereley S.T., Kompenhans J. Particle Image Velocimetry. Practical Guide. Berlin: Springer, 2007.
4. Elsinga G.E., Scarano F., Wieneke B., van Oudheusden B.W. Tomagraphic particle image ve-locimetry // Exp. Fluids. 2006. V. 41. P. 933-947.
5. Бильский А.В., Ложкин В.А., Маркович Д.М. и др. Оптимизация и тестирование томографического метода измерения скорости в объеме потока // Теплофизика и аэромеханика. 2011. Т. 18. № 4. C. 555-566.
6. MaasH.G., Gruen A.,Papantoniou,D. Particle tracking velocimetry in tree-dimentional flows // Exp. Fluids. 1993. V. 15. P. 133-146.
7. Novara M., Scarano F. Performances of motion tracking enhanced Tomo-PIV on turbulent shear flows // Exp. Fluids. 2012. V. 52. P. 1027-1041.
8. Wieneke B. Iterative reconstruction of volumetric particle distribution // Proc. 9th International Symposium on Particle Image Velocimetry, July 21-23, 2011, Kobe, Japan.
Бондаренко Анатолий Николаевич Дедок Василий Александрович
Институт математики СО РАН им. С.Л. Соболева (г. Новосибирск)
Токарев Михаил Петрович
Институт теплофизики СО РАН (г. Новосибирск)
E-mail: [email protected];
[email protected]; [email protected] Поступила в редакцию 30 апреля 2012г.
Bondarenko Anatoly N., Dedok Vasily A., Tokarev Mikhail P. (Sobolev Institute of Mathematics, Institute of Thermophysics, Siberian Branch of the Russian Academy of Sciences). Discrete scattering centers reconstruction problem on a limited set of projections.
Keywords: particle image velocimetry (PIV), optical tomography, inverse problems.
We consider the actual problem of reconstruction the maximum subset of points in R3 on a limited set of two-dimensional projections. This problem arises naturally in applications of physical hydrodynamics to real optical diagnostics of liquid and gas by measuring the instantaneous velocity fields in the bulk flow. During the experimental data processing is solved the problem of reconstruction the internal structure of the scattering medium which is a set of discrete small particles, moving with the speed closely like to the speed of a continuous medium. The instantaneous flow rate is determined by the change in the position of the particles over time. We prove a theorem that allows us to answer the question of the sufficiency of the changes to uniquely identify a set of two-dimensional set of its projections. Iterative algorithm is proposed to restore the maximum set by its set of projections.
