Набор матричных операций для суммирования вокруг элемента матрицы


У меня есть результаты моделирования, которые являются дискретными комплексными числами, представляющие собой волновой функции на сетке. Рассчитать фазу волновой функции. Для каждой точки На этой сетке я должен выполнять эти операции. Пример:

|a b c|
|d e f|
|g h i|

Сумма разницы между соседями вокруг элемента е:

 (a-b)+(b-c)+(c-f)+(f-i)+(i-h)+(h-g)+(g-d)+(d-a)

Но каждый отнимание δ иметь подвергают эти условия

if δ > π than δ = δ - 2π 
if δ < -π than δ = δ + 2π 

Я написал код в Octave, который extremaly неэффективным и медленным, но я надеялся, что есть способ, чтобы перейти от расчетов buteforce ускорить расчеты, используя некоторые функции в пакет изображения. Я знаю, что я могу перфом sumation используя свертку, но я все еще есть проблемы с другими операциями.

Код в Octave:

function deltaPhi = phaseDifference(phi1, phi2)
 deltaPhi = phi1 - phi2;
  if(deltaPhi > pi)
    deltaPhi = deltaPhi - 2*pi;
  endif

  if(deltaPhi < -pi)
    deltaPhi = deltaPhi + 2*pi;
  endif;
end

function [phase] = checkPhase(M)

phase = zeros(size(M)-2);

  for i = 2:size(M,1)-1
    for j = 2:size(M,2)-1      
      phase(i-1,j-1) = phaseDifference(M(i-1,j-1),M(i,j-1)) + phaseDifference(M(i,j-1),M(i+1,j-1)) + phaseDifference(M(i+1,j-1),M(i+1,j)) + phaseDifference(M(i+1,j),M(i+1,j+1)) + phaseDifference(M(i+1,j+1),M(i,j+1)) + phaseDifference(M(i,j+1), M(i-1,j+1)) + phaseDifference(M(i-1,j+1), M(i-1,j)) + phaseDifference(M(i-1,j), M(i-1,j-1));
    endfor
  endfor

end

Идея в том, чтобы переписать этот код в OpenCV и использовать некоторые из этих методов библиотеки.



122
1
задан 29 января 2018 в 02:01 Источник Поделиться
Комментарии
1 ответ

Во-первых, некоторые простые вещи:


  • Я хотел бы использовать end вместо endfor и endif, чтобы сохранить код, совместимый с MATLAB.

  • Не стоит в конце функции end (хотя это и не вредны). В end необходимо только при написании вложенных функций , которые имеют доступ к родительской функции рабочей области.

  • Вы должны разбить очень длинные строки кода, используя ... (см. пример ниже).

  • После перехода через большой матрицы, всегда сделать первый показатель вашего внутреннего цикла. Матричные элементы M(1,j) и M(2,j) являются смежными в памяти, M(i,1) и M(i,2) нет. Если вы переключите два вложенных цикла, ваш код будет быстрее на больших массивах, как вы используете кэш лучше.


Теперь самое сложное: изменение вашей логике.


|a b c|
|d e f|
|g h i|


Сумма разницы между соседями вокруг элемента e:

(a-b)+(b-c)+(c-f)+(f-i)+(i-h)+(h-g)+(g-d)+(d-a)

Обратите внимание, что для вычисления разности фаз между двумя соседними клетками несколько раз: например, (a-b) будет вычислено при определении стоимости для e, но и в предыдущем цикле итерации, когда d был центром элемента; это будут вычислены еще два раза, когда эта пара появляется в нижней части такого соседства.

Таким образом, вы должны предварительно вычислить эти различия и хранить их:

vdiff = phaseDifference(diff(M,1,1));
hdiff = phaseDifference(diff(M,1,2));

Для этого, чтобы работать, phaseDifference должны быть переписаны, см. ниже.

Теперь, vdiff имеет один меньше строк, и hdiff имеет один меньше столбцов, чем M. Мы должны это учитывать при индексации. Воспринимайте их как элементы между строк/столбцов M.

Сумма для одного элемента (i,j) теперь становится:

phase(i-1,j-1) = hdiff(i-1,j-1) + hdiff(i-1,j) ...
- hdiff(i+1,j-1) - hdiff(i+1,j) ...
- vdiff(i-1,j-1) - vdiff(i,j-1) ...
+ vdiff(i-1,j+1) + vdiff(i,j+1);


Полная функция теперь будет:

function phase = checkPhase2(M)
vdiff = phaseDifference(diff(M,1,1));
hdiff = phaseDifference(diff(M,1,2));
phase = zeros(size(M)-2);
for j = 2:size(M,2)-1
for i = 2:size(M,1)-1
phase(i-1,j-1) = hdiff(i-1,j-1) + hdiff(i-1,j) ...
- hdiff(i+1,j-1) - hdiff(i+1,j) ...
- vdiff(i-1,j-1) - vdiff(i,j-1) ...
+ vdiff(i-1,j+1) + vdiff(i,j+1);
end
end

Поскольку Октав медленно с петель, вполне вероятно, что это векторизованных версия будет быстрее:

function phase = checkPhase3(M)
vdiff = phaseDifference(diff(M,1,1));
hdiff = phaseDifference(diff(M,1,2));
j = 2:size(M,2)-1;
i = 2:size(M,1)-1;
phase = hdiff(i-1,j-1) + hdiff(i-1,j) ...
- hdiff(i+1,j-1) - hdiff(i+1,j) ...
- vdiff(i-1,j-1) - vdiff(i,j-1) ...
+ vdiff(i-1,j+1) + vdiff(i,j+1);

Я тестирую на R2017a MATLAB и векторные формы, это на самом деле медленнее, чем версия с петель, так как индексация по-прежнему довольно медленно, тогда как Матлаб добилась огромного успеха в прошлом году ускоряя код с петель.

Тайминги для ввода 2000х3000 элементы:


  • checkPhase(M) (ФП версии) работает в 1.3 с

  • и меняя порядок циклов он работает в 1.0 с

  • checkPhase2(M) работает в 0,17 с,

  • checkPhase3(M) работает в 0,25 сек.

Ваши тайминги на октаву будет сильно отличается (гораздо медленнее кода с петлями).


Для phaseDifference для работы с массивом входных данных, его нельзя использовать if кстати, что исходная функция делает. Можно переписывать той же логике, используя логическое индексирование для работы на полной матрице, но эта версия проще:

function deltaPhi = phaseDifference(deltaPhi)
deltaPhi = mod(deltaPhi + pi, 2*pi) - pi;

Обратите внимание, Я заменил двумя входами с одной, как diff уже вычисляет разницу между соседями.

1
ответ дан 1 февраля 2018 в 07:02 Источник Поделиться