Повышаем быстродействие расчётов на примере медианного фильтра
Рассмотрим некоторые приёмы повышения быстродействия вычислительных программ на примере алгоритма медианного фильтра.
Медианный фильтр — алгоритм, применяемый к изображению или подобной ему матрице, который значение каждого элемента заменяет на медиану среди его соседей. Сам по себе он может применяться для грубого сглаживания изображения, однако более полезен для устранения битых и горячих пикселей. Идея тут такая, что, если значение исходного элемента очень сильно отличается от значения медианного фильтра в той же точке, то значение элемента заменяется значением медианного фильтра.
В качестве исходного медианного фильтра используем свободно распространяемый исходный код профессора факультета математики Калифорнийского технологического университета Джона Буркардта. Этот код (gray_median_news) входит в модуль image_denoise и представлен на языках C, C++ и Fortran 90. Мы здесь рассмотрим оптимизацию кода на Фортране, однако практически всё то же самое применимо и к сишным версиям, различия чисто синтаксические.
Исходная процедура gray_median_news отлично работает (хотя содержит небольшую ошибку при обработке правого нижнего пикселя), однако довольно медленна. Посмотрим, что тут можно сделать.
Взглянем на главный цикл обработки данных в фильтре:
do i = 2, m - 1
do j = 2, n - 1
p(1) = gray(i-1,j )
p(2) = gray(i+1,j )
p(3) = gray(i ,j+1)
p(4) = gray(i ,j-1)
p(5) = gray(i ,j )
call i4vec_median ( 5, p, gray2(i,j) )
end do
end do
Во-первых, тут напрашивается единственная чисто фортрановская оптимизация. Данные в фортрановских массивах хранятся по столбцам, а циклы обходят их по строкам, поэтому для более эффективного регулярного доступа к памяти поменяем внутренний и внешний циклы местами:
do j = 2, n - 1
do i = 2, m - 1
p(1) = gray(i-1,j )
p(2) = gray(i+1,j )
p(3) = gray(i ,j+1)
p(4) = gray(i ,j-1)
p(5) = gray(i ,j )
call i4vec_median ( 5, p, gray2(i,j) )
end do
end do
Далее, обратим внимание, что обработка здесь заключается в пересчитывании данных из массива gray в массив gray2 и асинхронна по отношению к порядку элементов. Поэтому хотелось бы распараллелить и заасинхронить эти циклы. Но нам здесь будет мешать массив p. Конечно, можно его сделать приватным для разных веток, но вместо этого давайте сначала разберёмся с подпрограммой i4vec_median.
Что же профессор делает в i4vec_median? А вот что:
subroutine i4vec_median ( n, a, median )
implicit none
integer ( kind = 4 ) n
integer ( kind = 4 ) a(n)
integer ( kind = 4 ) k
integer ( kind = 4 ) median
k = ( n + 1 ) / 2
call i4vec_frac ( n, a, k, median )
return
end
Здесь он вызывает другую подпрограмму i4vec_frac, предназначенную для поиска k-го элемента отсортированного массива a длиной n, полагая k равным середине массива. Сама подпрограмма i4vec_frac выглядит так:
subroutine i4vec_frac ( n, a, k, frac )
implicit none integer ( kind = 4 ) n
integer ( kind = 4 ) a(n)
integer ( kind = 4 ) frac
integer ( kind = 4 ) i
integer ( kind = 4 ) iryt
integer ( kind = 4 ) ix
integer ( kind = 4 ) j
integer ( kind = 4 ) k
integer ( kind = 4 ) left
integer ( kind = 4 ) t
if ( n <= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'I4VEC_FRAC - Fatal error!'
write ( *, '(a,i8)' ) ' Illegal nonpositive value of N = ', n
stop
end if
if ( k <= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'I4VEC_FRAC - Fatal error!'
write ( *, '(a,i8)' ) ' Illegal nonpositive value of K = ', k
stop
end if
if ( n < k ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'I4VEC_FRAC - Fatal error!'
write ( *, '(a,i8)' ) ' Illegal N < K, K = ', k
stop
end if
left = 1
iryt = n
do
if ( iryt <= left ) then
frac = a(k)
exit
end if
ix = a(k)
i = left
j = iryt
do
if ( j < i ) then
if ( j < k ) then
left = i
end if
if ( k < i ) then
iryt = j
end if
exit
end if
!
! Find I so that IX <= A(I).
!
do while ( a(i) < ix )
i = i + 1
end do
!
! Find J so that A(J) <= IX.
!
do while ( ix < a(j) )
j = j - 1
end do
if ( i <= j ) then
t = a(i)
a(i) = a(j)
a(j) = t
i = i + 1
j = j - 1
end if
end do
end do
return
end
Разумеется, учинять сортировку массива произвольного размера — это совсем не то, что бы нам хотелось иметь для решения вполне конкретной задачи нахождения медианы из пяти чисел (то есть такого числа, что два других больше его, а ещё два других не больше). Дальше по ходу алгоритма профессора ещё используются медианы тройки чисел на краях матрицы.
Попробуем подумать, как нам найти медиану пятёрки и медиану тройки наиболее эффективным способом.
Первое, что приходит в голову — это тупо расписать все альтернативы упорядоченности чисел вложенными условными операторами. Например, для тройки это может выглядеть так:
pure integer function median3 (a1, a2, a3)
implicit none
integer, intent(in) :: a1, a2, a3
if (a1 > a2) then
if (a1 > a3) then
if (a2 > a3) then
median3 = a2
else
median3 = a3
end if
else
median3 = a1
end if
else if (a2 > a3) then
median3 = a3
else
median3 = a2
end if
end function median3
К сожалению, для пятёрки таким путём получится монструозная конструкция из 120 веток вложенных условных операторов, которую написать представляется возможным разве что путём автоматической генерации. Но надо ли это делать?
Вспомним, что современные процессоры включают наборы векторных инструкций SSE и подобные им, которые обрабатывают небольшие массивы гораздо эффективнее, чем процессор применял бы последовательные ветвления. Беглый анализ существующей практики показывает, что наиболее эффективными процедурами, использующими SSE, в нашем случае будут такие:
pure integer function median3 (a, b, c)
implicit none
integer, intent(in) :: a, b, c
median3 = max(min(a,b),min(c,max(a,b)))
end function median3
pure integer function median5 (a, b, c, d, e)
implicit none
integer, intent (in) :: a, b, c, d, e
integer: f, g
f = max (min (a, b), min (c, d))
g = min (max (a, b), max (c, d))
median5 = median3 (e, f, g)
end function median5
Можно практически убедиться, что функция median3, использующая max и min, действительно примерно вдвое быстрее, чем использующая if.
Перепишем наши циклы с использованием процедур median3 и median5:
do j = 2, n - 1
do i = 2, m - 1
gray2(i,j) = median5 (gray(i-1,j ), &
gray(i+1,j ), &
gray(i ,j+1), &
gray(i ,j-1), &
gray(i ,j ))
end do
end do
Обратите внимание, что мы избавились от массива p, так как можем передать параметры просто в 5 аргументах функции median5. И теперь нам ничего не мешает, действуя чисто формально, распараллелить и заасинхронить циклы обработки:
!$omp parallel do private (i)
do j = 2, n - 1
do concurrent (i=2:m-1)
gray2(i,j) = median5 (gray(i-1,j ), &
gray(i+1,j ), &
gray(i ,j+1), &
gray(i ,j-1), &
gray(i ,j ))
end do
end do
!$omp end parallel do
Аналогичным образом мы можем исправить и менее нагруженные циклы обработки точек на краях матрицы.
Почему мы параллелим внешний цикл и не параллелим внутренний? Потому что параллельности внешнего цикла вполне достаточно при разумном для обычного компьютера предположении о том, что размер массива по высоте превышает количество ядер процессоров. А накладные расходы на организацию потоков вычислений лучше нести однократно, для внешнего цикла.
Каков же итог этих оптимизаций? На стареньком тестовом компьютере автора работу фильтра удалось ускорить более чем в 130 раз, откомпилировав программу с использованием OpenMP. Для матрицы в 4 миллиона элементов это соответствует прогрессу от десятых до тысячных долей секунды. Фильтр стал пригоден, например, для обработки кадров в реальном времени.
Разумеется, быстрые алгоритмы фильтрации изображения имеются в специализированных библиотеках, так что эту статью не следует рассматривать как руководство по обработке изображений. Идея статьи в том, чтобы показать некоторые приёмы оптимизации и проиллюстрировать важность периодического подключения к процессу программирования мозга, а не только скачанных из интернета библиотек.
