Быстрый способ подсчета соседних точек в 3D-массиве

У меня есть код, в котором я хочу перебрать все точки в сетке и для каждой точки проверить, выполняется ли данное условие для достаточного количества соседних точек. Кроме того, у меня есть периодические границы на сетке.

Задача очень похожа на игру «Жизнь».

Мой текущий код выглядит примерно так

do k=1,ksize; do j=1,jsize; do i=1,isize ! Loop over all points
  ncount = 0
  kkloop: do kk=k-1,k+1 ! Loop over neighbours
    ktmp = kk
    if(kk>ksize) ktmp = 1 ! Handle periodic boundary
    if(kk<1) ktmp = ksize
    do jj=j-1,j+1
      jtmp = jj
      if(jj>jsize) jtmp = 1
      if(jj<1) jtmp = jsize
      do ii=i-1,i+1
        if(ii == 0 .and. jj == 0 .and. kk == 0) cycle ! Skip self
        itmp = ii
        if(ii>isize) itmp = 1
        if(ii<1) itmp = isize

        if(grid(itmp,jtmp,ktmp)) ncount = ncount + 1 ! Check condition for neighbour
        if(ncount > threshold) then ! Enough neigbours with condition?
          do_stuff(i,j,k)
          break kkloop
        end if
      end do
    end do
  end do
end do; end do; end do

Это не элегантно и, вероятно, не очень эффективно. Есть лучший способ сделать это? Этот код будет повторяться много раз, поэтому я хотел бы сделать его как можно быстрее.


optimization fortran
person pafcu    schedule 19.06.2012    source источник

Ответы (3)


arrow_upward
5
arrow_downward

Я сделаю это в 2D, предоставлю вам раздувать до 3D.

Первое, что я бы сделал, это дополнил массив ореолом глубины, равным глубине интересующего вас района. Итак, если ваш массив объявлен как, скажем,

real, dimension(100,100) :: my_array

и вас интересуют 8 непосредственных соседей каждой ячейки,

real, dimension(0:101,0:101) :: halo_array
.
.
.
halo_array(1:100,1:100) = my_array
halo_array(0,:) = my_array(100,:)
! repeat for each border, mutatis mutandis

Это сэкономит много времени на проверку границы, и это стоит сделать независимо от того, следуете ли вы следующему предложению или нет. Вы можете сделать это «на месте», если хотите, я имею в виду просто развернуть my_array, а не копировать его.

Для элегантного решения вы можете написать что-то вроде этого

forall (i=1:100,j=1:100)
   if (logical_function_of(my_array(i-1,j),my_array(i+1,j),my_array(i,j-1),my_array(i,j+1),...) then
      do_stuff(my_array(i,j))
   end if
end forall

Здесь logical_function_of() возвращает true, когда окрестности my_array(i,j) удовлетворяют вашим критериям. Я устал после перечисления соседей N, S, E, W, и для производственного кода я, вероятно, все равно напишу это как функцию индексов. По моему опыту, forall элегантен (для некоторых), но не так эффективен, как вложенные циклы.

person High Performance Mark    schedule 19.06.2012
comment
+1 за чистый код и за mutatis mutandis, которые я никогда раньше не видел в вычислительном контексте. - person Jonathan Dursi; 19.06.2012

arrow_upward
1
arrow_downward

Вы можете использовать kd-дерево или октодерево для разделения трехмерного массива. Кривая заполнения пространства, такая как кривая порядка z Мортона, полезна в 3D для создания ключа для восьми кубов. Но лучше всего это работает с мощностью 2 3d-массивов.

person Gigamegs    schedule 19.06.2012

arrow_upward
1
arrow_downward

Вы вызываете функцию grid 27 раз для каждой точки сетки. Если это дорогой звонок, вы хотели бы звонить по нему меньше раз.

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

В противном случае вы можете использовать битовую матрицу размерности isize*jsize*3 для кэширования значений сетки для всех точек с kk=k-1,k+1. Если это слишком много, вы можете выбрать промежуточное решение, используя битовую матрицу размера isize*3*3 для кэширования.

person salva    schedule 19.06.2012
comment
Каждый раз, когда я читаю код, grid(itmp,jtmp,ktmp) выглядит для меня оператором доступа к массиву, а не вызовом функции. - person High Performance Mark; 19.06.2012
comment
Особенно, когда вы перебираете ktmp в самом внешнем цикле и itmp в самом внутреннем, это выглядит как оператор доступа к массиву. - person Hristo Iliev; 20.06.2012
comment
Что ж, мы все только предполагаем. Чтение вопроса (проверка каждой точки при заданном условии) заставило меня подумать, что grid это функция. Посмотрим, прояснит ли это ОП. - person salva; 20.06.2012