Повышаем быстродействие расчётов на примере медианного фильтра

35eca0f0230c519cc0c0390034b8ffe7

Рассмотрим некоторые приёмы повышения быстродействия вычислительных программ на примере алгоритма медианного фильтра.
Медианный фильтр — алгоритм, применяемый к изображению или подобной ему матрице, который значение каждого элемента заменяет на медиану среди его соседей. Сам по себе он может применяться для грубого сглаживания изображения, однако более полезен для устранения битых и горячих пикселей. Идея тут такая, что, если значение исходного элемента очень сильно отличается от значения медианного фильтра в той же точке, то значение элемента заменяется значением медианного фильтра.

В качестве исходного медианного фильтра используем свободно распространяемый исходный код профессора факультета математики Калифорнийского технологического университета Джона Буркардта. Этот код (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 миллиона элементов это соответствует прогрессу от десятых до тысячных долей секунды. Фильтр стал пригоден, например, для обработки кадров в реальном времени.

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

© Habrahabr.ru