CC BY
22EFFECTIVE IMPLEMENTATION OF USFFT ALGORITHM
A. S. Matveev, V.V. Nikitin*, A. A. Romanenko**, A. A. Duchkov
IPGG SB RAS, 630090, Novosibirsk, Russia *NSU,
630090, Novosibirsk, Russia ** Center for Mathematical Sciences, Lund University, 22100, Lund, Sweden
This article is devoted to the Unequispaeed Fast Fourier Transform (USFFT), which is a popular analytical tool for solving physics and engineering problems. The most common applications of the transform include seismology, optics, computed tomography, crystallography, etc. Despite the favorable computational complexity of the USFFT algorithm (0(N logN)), the execution time remains rather high due to the algorithm structure and large input data sizes. There are two main types of USFFT: Fourier transform from equispaeed grid to unequispaeed grid and Fourier transform from unequispaeed grid to equispaeed grid. Corresponding computational algorithms consist of three main steps: convolution, Fast Fourier Transform (FFT) and deeonvolution. Profiling shows that up to 95% of execution time is spent on the convolution step. In this paper, we propose a parallel USFFT algorithm and its effective cache-optimized implementation on CPU for one-, two- and three-dimensional cases. Cache performance optimization is based on the sorting of unequispaeed grid points. The constructed sorting procedure sufficiently reduces the number of cache misses. For instance, for the two-dimensional case the number of cache misses is reduced by 36 times, which results in 2x speed-up of the transform evaluation. Next, we propose a parallel block algorithm for the convolution step and implement it by making use of OpenMP, a popular extension for the C programming language supporting multi-platform shared memory parallel programming. The obtained parallel implementation was optimized in terms of optimal block sizes and type of scheduling for the convolution step. Numerical tests show high parallel efficiency: speed-up on 16 processors compared to the seciuential implementation is approximately eciual to 13. The tests also show that the performance is several times higher than the performance of the commonly-used library for the fast Fourier transform at noneciuispaeed nodes (NFFT 3.0). USFFT is commonly used for fast evaluation of the Radon transform operator which is one of the main mathematical tools in computed tomography. In this paper, we consider a standard reconstruction of tomography data by inversion of the Radon transform, and an iterative reconstruction by using the Expectation-maximization algorithm. The iterative reconstruction is well-suited for processing data with [21] or irregularly-structured data. Since iterative schemes assume applying the forward and adjoint Radon operators several times, computational times for preprocessing procedures such as sorting of grid points and allocation memory can be diminished. The obtained program for evaluating iterative schemes was tested for synthetic Radon data containing Poisson noise. The program outperforms the implementation via NFFT by 4.4 times for the same accuracy level.
Key words: fast Fourier transform, uneciuispaeed grids, parallel algorithm, optimization, high performance computing.
References
1. Coolev J. W., Tukev J. W. An algorithm for the machine calculation of complex Fourier series // Mathematics of computation. 1965. T. 19. N 90. P. 297-301.
2. Intel Math Kernel Library (Intel MKL) [El. res.] // https://software.intel.com/en-us/ intel-mkl
3. cuFFT | NVIDIA Developer [El. res.] // https://developer.nvidia.com/cufft
4. FFTW Home page [El. res.] // http://www.fftw.org
5. Bracewell R. N. Strip integration in radio astronomy // Australian Journal of Physics. 1956. T. 9. N 2. P. 198-217.
6. Duchkov A. A., Andersson F., De Hoop M. V. Discrete almost-symmetric wave packets and multiscale geometrical representation of ([8]) waves // Geoscience and Remote Sensing, IEEE Transactions on. 2010. T. 48. N 9. P. 3408-3423.
7. Bevlkin G., Burridge R. Linearized inverse scattering problems in acoustics and elasticity // Wave motion. 1990. T. 12. N 1. P. 15-52.
8. Zwartjes P. M.. Sacchi M. D. Fourier reconstruction of nonuniformlv sampled, aliased seismic data // Geophysics. 2006. T. 72. N 1. P. V21-V32.
9. Dutt A., Rokhlin V. Fast Fourier transforms for nonequispaced data // SIAM Journal on Scientific computing. 1993. T. 14. N 6. P. 1368-1393.
10. Bevlkin G. On the fast Fourier transform of functions with singularities // Applied and Computational Harmonic Analysis. 1995. T. 2. N 4. P. 363-381.
11. Greengard L., Lee J. Y. Accelerating the nonuniform fast Fourier transform // SIAM review. 2004. T. 46. N 3. P. 443-454.
12. Fessler J. A., Sutton B. P. Nonuniform fast Fourier transforms using min-max interpolation // Signal Processing, IEEE Transactions on. 2003. T. 51. N 2. P. 560-574.
13. NUFFT page [El. res.] // http://www.cims.nyu.edu/cmcl/Clll/nufft.html
14. NFFT — TU Chemnitz [El. res.] // https://www-user.tu-chemnitz.de/$\sim$potts/nfft/
15. Andersson F. Algorithms for unequally spaced fast Laplace transforms // Applied and Computational Harmonic Analysis. 2013. T. 35. N 3. P. 419-432.
16. Herman G. T., Louis A. K., Natterer F. (ed.). Mathematical methods in tomography: proceedings of a conference held in Oberwolfach, Germany, 5-11 June, 1990. Springer, 2006.
17. Yilmaz O. Seismic data analysis. Tulsa : Society of exploration geophvsicists, 2001. T. 1. P. 74170-2740.
18. Tretiak O., Metz C. The exponential Radon transform // SIAM Journal on Applied Mathematics. 1980. T. 39. N 2. P. 341-354.
19. Natterer F. Inversion of the attenuated Radon transform // Inverse problems. 2001. T. 17. N 1. P. 113.
20. Shepp L. A., Logan B. F. The Fourier reconstruction of a head section // Nuclear Science, IEEE Transactions on. 1974. T. 21. N 3. P. 21-43.
21. Barrett H. H., Wilson D. W., Tsui B. M. W. Noise properties of the EM algorithm. I. Theory // Physics in medicine and biology. 1994. T. 39. N 5. P. 833.
22. Yan M., Vese L. A. Expectation maximization and total variation-based model for computed tomography reconstruction from undersampled data // SPIE Medical Imaging. International Society for Optics and Photonics, 2011. P. 79612X-79612X-8.
23. Champlev K. SPECT reconstruction using the expectation maximization algorithm and an exact inversion formula: /pic. MS Thesis, Oregon State University, 2004.
24. Dempster A. P., Laird N. M., Rubin D. B. Maximum likelihood from incomplete data via the EM algorithm // Journal of the royal statistical society. Series B (methodological). 1977. P. 1-38.
25. Miqueles E. X., Helou E. S., De Pierro A. R. Generalized Backprojection Operator: Fast Calculation // Journal of Physics: Conference Series. IOP Publishing, 2014. T. 490. N 1. P. 012148.
ЭФФЕКТИВНАЯ РЕАЛИЗАЦИЯ АЛГОРИТМА БЫСТРОГО ПРЕОБРАЗОВАНИЯ ФУРЬЕ НА НЕРЕГУЛЯРНЫХ СЕТКАХ
A.C. Матвеев, В. В. Никитин*, A.A. Романенко**, A.A. Дучков
Институт нефтегазовой геологии и геофизики СО РАН, 630090, Новосибирск, Россия * Новосибирский национальный исследовательский государственный университет,
630090, Новосибирск, Россия
**
22100, Лупд, Швеция
УДК 621
Статья посвящена преобразованию Фурье на нерегулярных сетках (USFFT), популярному средству анализа во многих естественнонаучных задачах. Большинство практических задач, использующих USFFT, имеют большой объем данных, что приводит к значительным вычислительным затратам. В данной работе предложена реализация алгоритма USFFT, использующая такие особенности современных центральных процессоров как параллелизм и наличие большого кэша данных. Оптимизация последовательной программы позволила сократить время выполнения наиболее трудоемкого этана преобразования в два раза, а последующее распараллеливание дало тринадцатикратное ускорение на вычислительном узле с 16 ядрами.
Ключевые слова: быстрое преобразование Фурье, нерегулярные сетки, оптимизация, параллелизм.
Введение. Преобразование Фурье является популярным средством анализа сигналов и изображений. В настоящее время в связи с повсеместным использованием вычислительной техники наиболее часто используется дискретный аналог данного преобразования — дискретное преобразование Фурье (ДПФ), Отчасти это связано с существованием алгоритма быстрого преобразования Фурье (БПФ или FFT) |1|, Алгоритм БПФ предполагает равномерную дискретизацию данных и имеет вычислительную сложность O(N log N) вместо O(N2) в случае ДПФ (N — размер входных данных). Данный алгоритм реализован во многих библиотеках, например Intel MKL |2|, cnFFT |3| и FFTW |4|,
Однако существует ряд задач, дня которых условие регулярности входных данных не выполняется. Например, возникает необходимость вычисления ДПФ в полярных координатах в задачах радиоастрономии |5|. В геофизике преобразование Фурье на нерегулярных сетках нашло применение дня сжатия данных по волновым пакетам |6| и миграции Стол та |7|. Также существует задача сиектралыюго анализа сейсмических данных, полученных с нерегулярно расположенных приемников |8|.
Учитывая большой объем данных, характерных дня таких задач, отсутствие быстрого алгоритма дня нерегулярных данных представляло собой проблему. Это послужило причиной того, что был разработан вычислительный алгоритм дня выполнения быстрого преобразования Фурье на нерегулярных сетках |9, 10| (USFFT, в некоторых источниках
NUFFT [11, 12]) вычислительная сложность которого совпадает со сложностью FFT и равна O(N log N),
С вычислительной точки зрения возникает задача эффективной реализации данного алгоритма, которая бы максимально использовала возможности современных центральных процессоров, многообразие и возможности которых постоянно растут и изменяются. Одним из направлений, благодаря которым растет их производительность, является увеличение количества ядер. Совокупность этих факторов приводит к задаче создания и оптимизации параллельных алгоритмов USFFT,
В настоящее время появляются первые библиотеки, реализующие данное преобразование, Самыми известными из них являются NUFFT, разработанная в Кураiпоиском институте математических наук (подразделение Нью-Йоркского университета) [13], и NFFT [14], разработанная в Хемницком техническом университете. Библиотека NUFFT сделана для образовательных целей и является последовательной, в свою очередь для NFFT имеется параллельная версия, В данной статье будет произведено сравнение производительности библиотеки NFFT с полученной реализацией алгоритма USFFT,
1. Алгоритмы быстрого преобразования Фурье на нерегулярных сетках. Рассмотрим ДПФ на нерегулярных сетках, заданные следующим образом:
N-1
Fk = £ fje-2nikxj, (1)
j=0
N/2-1
/з = £ Рке-2п1кх;. (2)
к=^/ 2
Здесь рациональные числа {х3- Для К0Т0Рых |х3-1 < 1/2, задают нерегулярную сетку, а целые числа к = -N/2,... ,N/2 — 1 задают регулярную сетку, 3начения / принадлежат множеству комплексных чисел. Вычисление сумм в формулах (1) и (2) напрямую выполняется за О^2) операций, в случаях двумерного и трехмерного аналогов данного преобразования вычислительная сложность равна О^4) и О^6) соответственно.
Рассмотрим эвристическое описание работы алгоритмов иЯГГТ для быстрого вычисления преобразований (1) и (2), более подробная формулировка и оценки точности представлены в [9], Начнем с того, что экспоненциальные множители внутри сумм могут быть представлены в виде свертки с дельта-функцией,
/те
е-2ткх 8(х — х3 )йх. (3)
-те
Тогда сумма в (1) может быть представлена в следующем виде:
N-1 »те N-1
£ ¡1 е-2ткх = £ ¡16(х — х3)в-2ткх¿х. (4)
3=0 J-те 1=0
Также введем гауссиан в форме ф(х) = в-Хх\ для некоторого пара метра А, контролирующего ширину гауссиана, Домножим выражение (4) с обеих сторон на преобразование Фурье от гауссиана, обозначенное как ф(к), тогда
N-1 N-1
fje-2mkXj = Е fj- Xj)e-2nikxdx (5)
j=0 j=0
или
N-1 - ^ N-1
\ Л .f „—2nikxj _ 1 I \ Л f JJ„„ ™ \„-2nikx,
j=о ф(к) J-^ j=o
^ fje-2mkxi = ^- ^ fjф(х - Xj)e dx (6)
отметить тот факт, что имеется возможность вычислять значение функции ф (^ — х только для таких значений аргумента — хЛ, которые находятся в небольшой (отно-
Интеграл в (6) может быть аппроксимирован при помощи правила трапеции, потому что получаемая под интегралом функция является гладкой за счет свертки с гауссианом. Поэтому достаточно точные результаты аппроксимации интеграла при помощи суммы
x
l 1
ответствии с фактором v, называемым oversampling facto г): введем {——:} _2 vN в качестве
vN l=-VN
x
Таким образом, алгоритм преобразования Фурье на нерегулярных сетках может быть представлен как последовательное применение операции дискретной свертки с гауссианом, БПФ и операции деконволюции (операция, обратная свертке, необходима для получения
спектра исходного сигнала, который был модифицирован в результате свертки), Стоит
" — - х-
vN X3
yvN xjJ, которые находятся в неболыно xj
рассматривать как имеющую ограниченный численный носитель. Причем размер окрестности зависит только от наперед заданной точности вычислений е [9], В таком случае операция дискретной свертки становится эквивалентной операции взвешенного рассеивания,
В итоге формула аппроксимации для (1) выглядит следующим образом
n-1 1 1 *2г-1 /N-1 / , N\
Fk = gfje_2nikXj-¿¿Д jfj*(VN-x-))e"2ni- (7)
где сумма по переменной l состоит из L = 2M + 1 ^ N (параметр M будет введен в Алгоритме 1 ненулевых элементов в виду ограниченного численного носителя функции ф(х).
На практике используется oversampling factor v = 2, так как он является оптимальным с точки зрения вычислительных ресурсов и требуемой точности.
Алгоритмы 1 и 2 составлены в соответствии с формулой (7), Отметим, что операторы
xj - xj
для преобразования (2) может быть построен аналогичным образом, заменив операцию взвешенного рассеивания на взвешенный сбор и выполняя шаги алгоритма 1 в обратном порядке.
Алгоритм 1. USFFT с нерегулярной сетки на регулярную
1) Взвешенное рассеивание значений точек нерегулярной сетки на точки регулярной сетки в шаре радиуса M = \yj2,25/пlogе][15], где \х] — округление х к наименьшему большему целому, е — требуемая точность вычислений с центром в точке Vj = [2Nxj], где [х] — округление х к ближайшему целому. Весовая функция — гауссиан ф(х) = e-Ax , где
д = - 2n2N2 " " " "
log £ "
Рис. 1. Возможное взаимное расположение регулярной и нерегулярной сетки
2) Дискретное преобразование Фурье на регулярный сетке размера 2 N, на которую произвелось взвешенное рассеивание. Может быть выполнено с помощью алгоритма БПФ.
3) Корректирующий фильтр к шагу 1, т.е. умножение на 1/ф(к).
Алгоритм 2. USFFT с регулярной сетки на нерегулярную.
1) Корректирующий фильтр к шагу 3, значение А вычисляется так же, как в преобразовании с нерегулярной сетки на регулярную
2) Дискретное преобразование Фурье на исходной регулярный сетке размера 2N. Может быть выполнено с помощью алгоритма БПФ.
3) Взвешенный сбор значений точек регулярной сетки, находящихся в шаре радиуса М с центром в точке Vj = [2Nxj], весовая функция идентична аналогу из операции взвешенного рассеивания.
На рис. 1 представлено возможное взаимное расположение регулярной и нерегулярной сетки. Квадратами изображены точки регулярной сетки, крестиками — нерегулярной. Красным цветом выделена точка нерегулярной сетки для которой производится операция взвешенного сбора или рассеивания, красный квадрат — ближайшая к ней точка регулярной сетки. Пунктиром обведены точки регулярной сетки, принадлежащие шару радиуса М = 2.
Операции взвешенного рассеивания и сбора имеют вычислительную сложность O(MN), преобразование Фурье имеет сложность O(NlogN), и корректирующий фильтр может быть выполнен за O(N) операций. В итоге вычислительная сложность каждого из алгоритмов равна 0(N log N + MN). что для М N аппроксимируется как 0(N log N).
2. Оптимизация последовательной программы. В ходе профилирования последовательной версии преобразований обоих типов было установлено, что до 95 % времени занимают этапы взвешенного рассеивания и взвешенного сбора. Эти этапы характеризуются малым количеством вычислений на каждое обращение в память, что приводит к большой зависимости времени вычисления от эффективности использования кэша.
Так как в общем случае точки нерегулярной сетки распределены случайным образом, доступ к точкам регулярной сетки, участвующим в сборе или рассеивании, является хаотическим. Следствием этого является большое количество кэш-промахов. Для того чтобы их уменьшить, необходимо организовать обход точек нерегулярной сетки таким образом, чтобы каждая последующая точка нерегулярной сетки имела как можно больше общих точек регулярной сетки совместно с предыдущей. Для достижения этого точки нерегу-
Таблица 1
Количество промахов в кэш последнего уровня
Размер сетки Без сортировки С сортировкой Идеальный случай
224 4,6е+06 2,9е+07 1,1е+07
212 х 212 9,5е+08 :]. !е 1)7 2,6е+07
28 х 28 х 28 7,6е+09 1,4е+08 •"). !с 1)7
лярной сетки были отсортированы в порядке возрастания координат. При многомерном преобразовании происходит последовательная сортировка сначала по одной координате, затем по остальным.
Данная оптимизация позволила уменьшить время взвешенного сбора в 2 раза при двумерном преобразовании и в 5 раз в трехмерном, В случае одномерного преобразования накладные расходы по сортировке оказались больше выигрыша, полученного в результате оптимизации попаданий в кэш.
Также было проверено то, насколько может возрасти производительность при такой организации доступа в память, при которой количество кеш-промахов было бы минимальным, Тестовый пример, в котором все точки нерегулярной сетки имели одну координату, что гарантировало присутствие в кэше задействованных в преобразовании данных, показал, что в идеальном случае количество кэш-промахов для преобразования может быть уменьшено в 36 раз (табл. 1) (данные приведены для двумерной версии преобразования с регулярной на нерегулярную сетку размера 212 х 212 для процессора Intel i7-4770),
3. Параллельный алгоритм USFFT, Аналогично этапу оптимизации последовательной версии программы, основное внимание было уделено этапам взвешенного сбора и рассеивания.
При выполнении взвешенного сбора цикл по точкам нерегулярной сетки возможно естественным образом распараллелить, используя технологию ОрепМР, Каждый поток использует точки нерегулярной сетки только для чтения, поэтому распараллеливание в этом случае не представляет сложности.
При этом установлено, что на количество попаданий в кэш последнего уровня прямым образом влияет политика планирования выполнения потоков. Это особенно актуально для трехмерного сбора, который характеризуется обработкой большого количества точек регулярной сетки, приходящихся на каждую точку нерегулярной.
Планирование выполнения необходимо организовать таким образом, чтобы в каждый момент времени различные потоки обрабатывали как можно более локальные данные, для предотвращения конфликта в кэше последнего уровня, который является общим для всех ядер в рамках одного процессора. Для этого необходимо разбить итерации цикла на блоки небольшого размера и установить динамическое распределение блоков между потоками. Оптимальный размер такого блока напрямую зависит от характеристик конкретного вычислительного устройства. На графике такой зависимости (рис, 2) для преобразования размером 128 х 128 х 128 для системы с двумя процессорами Intel Xeon Е5-2690 видно, что оптимальными являются блоки размера от 4 до 32 элементов.
Тестирование показало, что такой подход к распараллеливанию является достаточно эффективным — при использовании максимального количества доступных вычислитель-
Рис. 2. Зависимость времени выполнения от размера блока
Таблица 2
распараллеливания этапа взвешенного сбора
Размер сетки Количество потоков
2 4 8 16
224 0,97 0,93 0,89 0,85
228 0,95 0,89 0,86 0,85
212 х 212 0,95 0,90 0,87 0,87
213 х 213 0,95 0,89 0,87 0,87
28 х 28 х 28 0,95 0,90 0,87 0,83
ных ядер эффективность распараллеливания (отношение ускорения к числу используемых потоков) достигает 0,87 (ускорение в 13,9 раз на 16 потоках — табл. 2).
Однако, такое распараллеливание цикла по точкам нерегулярной сетки невозможно при выполнении взвешенного рассеивания. Это связано с тем, что возможна ситуация, когда два потока будут производить запись в один и тот же участок памяти, что приведет к ошибке (race condition). Пример такой ситуации изображен на рис. 3. Видно, что существуют точки регулярной сетки, принадлежащие к окрестностям (области ограниченные красным пунктиром) двух одновременно обрабатываемых точек нерегулярной сетки (красные крестики).
Для решения данной проблемы было предложено разбить множество точек регулярной сетки на группы блоков. Так, чтобы блоки, принадлежащие к одной группе, были расположены на расстоянии 2M + 1 друг от друга, где M — радиус рассеивания.
х ~ х
Тх1.1. ЛД * хл " - - • " * -
• я ■ >■ а ■ ■ ■ а,
I ЕЙ Щ Ш
! х ;
Й Ш
Я -В 9 л , 6 я. в
Т—т—Т—Т—Т—Т—Т—^— * - - • Ж---'
К X
ттт
X
—-ш-X
Рис. 3. Пример пересечения областей взвешенного рассеивания
х х X х х * х * X Л х V • V XX .1 х * X* X 1 * у Х>* х X Л XX •* XV* *** X ХхкХ*
X а' ^ ^С * X X X X *** ^V * X г X
' X х . 'X X Су" х ■ у X X
Xх X Xх Xх* X х х X * X* х (' х Xх XX ^ V
1 у X I1 у X х* 1 * X X 1 * , V V К * X # V " X
Рис. 4. Пример разбиения двумерной сетки
Тогда для корректного и эффективного выполнения взвешенного рассеивания необходимо распараллелить цикл обработки блоков, принадлежащих одной группе. Здесь под обработкой блока понимается проведение операции взвешенного рассеивания для всех точек нерегулярной сетки, принадлежащих этому блоку.
На рис. 4 изображен пример такого разбиения. Цветными крестиками обозначены точки нерегулярной сетки, принадлежащие одной группе блоков, различные цвета этих крестиков обозначают принадлежность к различным блокам. Цветные пунктирные линии обозначают участки памяти, куда может быть произведена запись при обработке блока соответствующего цвета. Видно, что эти области не пересекаются, а значит различные блоки могут быть обработаны независимо.
При равномерной плотности распределения точек нерегулярной сетки данный подход достигает эффективность распараллеливания, равную 0,86 (ускорение в 13,7 раза на 16 потоках — табл. 3). Меньшая по сравнению с операций взвешенного рассеивания эффективность распараллеливания связана в первую очередь с необходимостью барьерной синхронизации после обработки группы блоков.
4. Сравнение с существующими аналогами. Важным этапом при разработке параллельного алгоритма и реализации на его основе библиотеки (далее обозначена как библиотека ТТЯГРТ), нацеленной на высокие показатели производительности, является ее сравнение с существующими аналогами. В настоящее время наиболее распространенной реализацией преобразования Фурье на нерегулярных сетках является библиотека КП-ТТ. Сравнение производилось на двусокетном сервере с установленными на нем процессора-
Таблица 3
распараллеливания этапа взвешенного рассеивания
Размер сетки Количество потоков
2 4 8 16
224 0,96 0,90 0,86 0,86
228 0,95 0,90 0,86 0,86
212 х 212 0,95 0,89 0,85 0,85
213 х 213 0,94 0,89 0,85 0,85
28 х 28 х 28 0,95 0,87 0,81 0,70
Рис. 5. Время вьшолнения двумерного преобразования Фурье с нерегулярной сетки на регулярную библиотек иЗГЕТ и МГГТ при различных размерностях
ми Intel Xeon Е5-2690. В качестве тестов использовалась двумерная задача со случайно распределенными точками нерегулярной сетки, количество регулярных и нерегулярных точек совпадает.
Стоит отметить, что в библиотеке NFFT используется фиксированное значение параметра M = 12 для вычислений с двойной точностью. Учитывая высокую зависимость времени вычисления этапов взвешенного суммирования и рассеивания от этого параметра (график зависимости представлен на рис. 7), все тесты производительности выполнены с таким значением.
Как видно из графиков на рис. б, преобразования обоих типов из библиотеки USFFT до трех раз быстрее своих аналогов из библиотеки NFFT для некоторых размерностей сеток.
горитмов USFFT была протестирована для решения задачи обращения преобразования
Рис. б. Время вьшолиения двумерного преобразования Фурье с регулярной сетки на нерегулярную библиотек иЗГЕТ и NFFT при различных размерностях
Радона. Данное преобразование применяется во многих областях, например в компьютерной томографии [16], радиоастрономии [5] или обработке сейсмических данных [17].
Преобразование Радона задается следующим образом: пусть ¡(х,у) — функция двух действительных переменных, определенная на всей плоскости и достаточно быстро убывающей на бесконечности. Тогда преобразованием Радона функции /(.г.у ) называется функ-
К! (в,в)
! (в сое в г вт в, в вт в + г сое в)в,г.
(8)
Производится вычисление интеграла вдоль прямых, отстоящих на расстоянии з от центра координат, с углом наклона в. Схематически преобразование Радона для задачи томографии представлено на рис. 8. Лучи проходят через объект, затухают (функция затухания /(х,у)) и записываются приемником для разных значений 0 и
5.1. Метод фильтрованных обратных проекций. На практике большое значение имеет обратная задача — по заданной функции Радона определить исходную функцию, для задачи томографии — внутреннюю структуру объекта. Одним из наиболее популярных методов восстановления функции ¡(х,у) является метод фильтрованных обратных проекций [18], [19]. В таком случае формула для восстановления функции выглядит следующим образом:
! = К* ^ К!, (9)
где сопряженный оператор Н* к преобразованию Радона задан следующим образом:
Зависимость времени вычисления 20 преобразований от
значения параметра М
0,7
0,6
0,5
" 0,4 ОС Ol ä о,з
0,2 ^---
ОД ----"
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 М С регулярной на нерегулярную С нерегулярной на регулярную
Рис. 7. Зависимость времени вычисления от параметра M
[■ п/2
R*g(x,y) = g(x cos в + y sin в,в)йв (10)
J п/2
Оператор свертки W действует только по переменной s и задан в частотной области через функцию w(a) = |<т|, также называемую инверсионным фильтром.
Решение задачи восстановления томографических данных при помощи метода фильтрованных обратных проекций кратко описывается следующим образом. Для данных Rf необходимо из одномерных Фурье-образов Д/(<тД) сформировать двумерный спектр Rf(a,9) на полярной сетке, применить фильтр (умножение в частотной области на функ-
а)
Рис. 9. Радон-с
б)
и восстановленное изображение (б) фантома Шеппа-Логана
цию w(a) для каждого в), а затем выполнить обратное двумерное преобразование Фурье в полярной системе координат.
Программа, реализующая данное преобразование, была реализована и протестирована на изображении, которое является стандартом для проверки алгоритмов восстановления изображения: Фантом Шеппа-Логана [20].
5.2. Восстановление данных с шумом. Томографические данные могут быть неполными, содержать шум или другие артефакты, возникающие при выполнении измерений. В этих случаях необходимо использовать итерационные схемы, так как метод фильтрованных обратных проекций не гарантирует стабильность решения. В большей части проблема нестабильности решения связана с инверсионным фильтром (IT) из формулы (9), так как его применение (умножение на w(a) = |<т|) увеличивает высокочастотные компоненты, которые обычно представляют шум и другие резкие колебания в данных, обусловленные ошибками измерений.
Итерационные схемы восстановления томографических данных в большинстве случаев основываются на многократном применении прямого и сопряженного операторов преобразования Радона и не подразумевают применения инверсионного фильтра. Итерационные алгоритмы достаточно ресурсозатратны в случаях, где для сходимости необходимо очень большое число итераций. Поэтому разработка быстрых алгоритмов для прямого и сопряженного операторов Радона является очень важной задачей.
Рассмотрим задачу восстановления томографических данных с шумом. Более обоснованным будет рассмотреть шум с распределеним Пуассона [21], [22]. Одним из наиболее простых методов для решения задачи восстановления данных, содержащих такого рода шум, является ЕМ-алгоритм (англ. Expectation-maximization (ЕМ) algorithm) [23], [24], [25]. Для томографических данных g итерационная схема может быть описана следующим образом:
Рис. 10. Выполнение преобразования Радона для данных с шумом: а) входные данные, б) результат метода фильтрованных обратных проекций, в) результат 50 итераций ЕМ алгоритма
/
к+1
/
Я* I
к Я I я/к
Я*Хс
(11)
где функция х(8$) равна 1, если линия, параметризованная при помощи (в, в), проходит через единичный круг (носитель /), в противном случае функция равна 0. По формуле видно, что на каждом итерационном шаге нужно вычислять прямое преобразование Я/к
и сопряженное Я* (я/х
На рис. 10, а, представлены данные с шумом, имеющим распределение Пуассона, рис. 10, б, показывает результат применения метода фильтрованных обратных проекций. Отметим, что уровень шума в восстановленном изображении гораздо выше. В свою очередь, восстановление при помощи 50 итераций ЕМ-алгоритма (рис. 10, в) практически полностью удалило шум.
Результаты производительности итерационного алгоритма, реализованного с помощью библиотек МРРТ и 118РРТ, представлены на рис. 11. Из графиков видно, что реализация с помощью библиотеки иЭРРТ быстрее в 4,4 раза для различного количества итераций алгоритма.
Заключение. Учет особенностей современных микропроцессоров, таких как наличие большого кэша данных и нескольких вычислительных ядер в рамках одного узла, позволил эффективно реализовать библиотеку быстрого преобразования Фурье на нерегулярных сетках.
Основное внимание при оптимизации было уделено наиболее вычислительно-сложным этапам преобразования — взвешенному рассеиванию и сбору. Были предложены подходы к параллельной реализации данных этапов. Благодаря этому удалось достичь ускорения до 13,7 и 13,9 раз соответственно на шестнадцатиядерном узле. Производительность полученной реализации оказалась до 3 раз выше по сравнению с реализацией из библиотеки ЫРРТ при сравнении на синтетических тестах.
Факт наличия эффективной реализации алгоритма ЦБРРТ особенно важен при реализации итеративных методов. Примером такой задачи является процедура обра-
Зависимость времени выполнения EM алгоритма от
25 20 числа итераций
20,67
15,5
и 15 10
I s Oi a 10,33
CO 5,17 3,53
5 0 1,18 2,35
10 20 30 40 Количество итераций -USFFT -NFFT
Рис. 11. Зависимость времени выполнения ЕМ алгоритма от числа итераций
щения преобразования Радона с помощью ЕМ-алгоритма. Тестирование показало, что в данной задаче использование полученной реализации позволяет добиться более чем четырехкратного прироста производительности по сравнению с реализацией на основе библиотеки NFFT.
Список литературы
1. Cooley J. W., Tukey J. W. An algorithm for the machine calculation of complex Fourier series // Mathematics of computation. 1965. V. 19. N 90. P. 297-301.
2. Intel Math Kernel Library (Intel MKL). [Electron res.], https://software.intel.com/en-us/ intel-mkl
3. cuFFT | NVIDIA Developer. [Electron res.], https://developer.nvid.ia.com/cufft
4. FFTW Home page. [Electron res.], http://www.fftw.org
5. Bracewell R. N. Strip integration in radio astronomy // Australian Journal of Physics. 1956. V. 9. N 2. P. 198-217.
6. Duchkov A. A., Andersson F., De Hoop M. V. Discrete almost-symmetric wave packets and multiscale geometrical representation of ([8]) waves // Geoscience and Remote Sensing, IEEE Transactions on. 2010. V. 48. N 9. P. 3408-3423.
7. Beylkin G., Burridge R. Linearized inverse scattering problems in acoustics and elasticity // Wave motion. 1990. V. 12. N 1. P. 15-52.
8. Zwartjes P. M.. Sacchi M. D. Fourier reconstruction of nonuniformly sampled, aliased seismic data // Geophysics. 2006. V. 72. N 1. P. V21-V32.
9. Dutt A., Rokhlin V. Fast Fourier transforms for nonequispaced data // SIAM Journal on Scientific computing. 1993. V. 14. N 6. P. 1368-1393.
10. Beylkin G. On the fast Fourier transform of functions with singularities / / Applied and Computational Harmonic Analysis. 1995. V. 2. N 4. P. 363-381.
11. Greengard L., Lee J. Y. Accelerating the nonuniform fast Fourier transform // SIAM review. 2004. V. 46. N 3. P. 443-454.
12. Fessler .J. A., Sutton B. P. Nonuniform fast Fourier transforms using min-max interpolation /7 Signal Processing, IEEE Transactions on. 2003. V. 51. N 2. P. 560 574.
13. NUFFT page. [Electron res.]. http://www.cims.nyu.edu/ciiicl/Clll/nufft.htiiil
14. NFFT TU Chemnitz. [Electron res.], https://www-user.tu-chemnitz.de/~potts/nfft/
15. Andersson F. Algorithms for unequally spaced fast Laplace transforms /7 Applied and Computational Harmonic Analysis. 2013. V. 35. N 3. P. 419 432.
16. Herman G. T., Louis A. K., Nattcrcr F. (cd.). Mathematical methods in tomography: proceedings of a conference held in Oberwolfaeh, Germany, 5 11 .June, 1990. Springer, 2006.
17. Yilmaz O. Seismic data analysis. Tulsa : Society of exploration geophysicists, 2001. V. 1. P. 74170 2740.
18. Trctiak O., Metz C. The exponential Radon transform /7 SIAM .Journal on Applied Mathematics. 1980. V. 39. N 2. P. 341 354.
19. Nattcrcr F. Inversion of the attenuated Radon transform /7 Inverse problems. 2001. V. 17. N 1. P. 113.
20. Shepp L. A., Logan B. F. The Fourier reconstruction of a head section /7 Nuclear Science, IEEE Transactions on. 1974. V. 21. N 3. P. 21 43.
21. Barrett H. H., Wilson D. W., Tsui B. M. W. Noise properties of the EM algorithm. I. Theory /7 Physics in medicine and biology. 1994. V. 39. N 5. P. 833.
22. Yan M., Vese L. A. Expectation maximization and total variation-based model for computed tomography reconstruction from undcrsamplcd data /7 SPIE Medical Imaging. International Society for Optics and Photonics, 2011. P. 79612X 79612X-8.
23. Champley K. SPECT reconstruction using the expectation maximization algorithm and an exact inversion formula: diss. MS Thesis, Oregon State University, 2004.
24. Dempster A. P., Laird N. M., Rubin D. B. Maximum likelihood from incomplete data via the EM algorithm /7 .Journal of the royal statistical society. Scries B (methodological). 1977. P. 1 38.
25. Miqueles E. X., Helou E. S., De Pierre A. R. Generalized Backprojcction Operator: Fast Calculation /7 .Journal of Physics: Conference Scries. IOP Publishing, 2014. V. 490. N 1. P. 012148.
Матвеев Алексей Сергеевич получил степень бакалавра в 2014 году и степень магистра в 2016 году в Новосибирском государственном университете но специальности „Информационные технологии". С 2014 года но настоящее время работает в Институте нефтегазовой геологии и геофизики СО РАН, работал в лаборатории динамических проблем сейсмики, возглавляемой Антоном Альбертовичем Дучковым. Его научными интересами являются высокопроизводительные вычисления, параллельные алгоритмы и обработка геофизических данных.
Aleksey Matveev received the В. S. degree in 2014 and the M.S. degree in 2016 in Information Technology, Novosibirsk State University, Novosibirsk, Russia. From 2014 to present, he has been working at the Trofimuk
Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, Russia, where he was a member of the laboratory of Dynamic Problems in Seismic led by prof. Anton Duchkov. His research interests include high-performance computing, parallel algorithms and seismic data processing.
Никитин Виктор Валерьевич получил степень бакалавра в 2011 году и степень магистра в 2013 году в Новосибирском государственном университете, обе но специальности „Информационные технологии". С 2009 по 2013 годы он был инженером в Институте нефтегазовой геологии и геофизики СО РАН и работал в лаборатории динамических проблем сейсмики, возглавляемой Антоном Альбертовичем Дучковым. В настоящее время учится в аспирантуре Лундского университета в Швеции, где также
является членом исследовательской группы по решению обратных задач, возглавляемой Фред-риком Андерссоном. Список ei'o научных интересов включает математические методы обработки сигналов, восстановление сейсмических и томографических данных, высокопроизводительные вычисления и оптимизация кода.
Viktor V. Nikitin received the B.S. degree in 2011 and the M.S. degree in 2013, both in Information Technology from Novosibirsk State University, Novosibirsk, Russia. From 2009 to 2013 he was an engineer at the Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, Russia, where he was a member of the laboratory of Dynamic Problems in Seismic-led by prof. Anton Duehkov. He has worked as a Ph.D. student in Mathematics at Lund University, Sweden, since then. He is also a member of the Inverse Problems group led by prof. Frcdrik Andersson from the Centre for Mathematical Sciences, Lund University. His research interests include image processing, reconstruction techniques in seismology and tomography, high-performance computing, and code optimization.
Романенко Алексей Анатольевич, зав. отделом компьютерной техники ФИТ ИГУ. Доцент кафедры систем информатики ФИТ ИГУ Канд. ^ | техн. наук, 2004. Интересы: ма-
тематическое моделирование, высокопроизводительные системы, адаптация программ и алгоритмов под современные вычислительные архитектуры.
Alexey Anatolievich Romanenko. Head of IT division of Information technologies department, Novosibirsk State University. Associate Professor, Novosibirsk State University. PhD in Computer science, 2004. Interests: mathematical modeling, HPC, programs and algorithms adaptation to modern computer systems.
Дучков Антон Альбертович, зав. лаборатории динамических проблем еейемики ИНГГ СО РАН. Доцент кафедры геофизики ГГФ ИГУ Канд. физ.-мат. наук, 2004. Интересы: исследования в области эффективных а;и'оритмов обработки сейсмических данных, сейсмической миграции, рмуляризации сейсмических данных, микроеейеми ческохх) мониторинга.
Anton Albertovich Duehkov. Head of Laboratory of Dynamic Problems in Seismic, Trofimuk Institute of Petroleum Geology and Geophysics SB RAS. Associate Professor, Novosibirsk State University. PhD in Geophysics, 2004. Interests: microseismic monitoring, seismic-data processing, regularization, and imaging.
Дата поступления 23.06.2016
