trjhtr

I have a results of simulations, which are discrete complex numbers representing wave function on a NxN grid. I calculate phase of the wavefunction. For every point on that grid I have to perform this set of operations. Example:

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

Sum of difference between neighbours around element e:

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

But every substraction δ have is subjected to this conditions

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

I wrote a code in Octave, which is extremaly inefficient and slow, but I was hoping there is a way to move from buteforce calculations to faster calculations using some of the functions in image package. I know I can perfom sumation using convolution but I still have problem with other operations.

Code in 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

The idea is to rewrite this code in OpenCV and use some of that librarys methods.

First some simple things:

• I would use end instead of endfor and endif, to keep the code compatible with MATLAB.




• It is not necessary to end a function with end (though it's not harmful either). The end is only necessary when writing nested functions that have access to the parent function's workspace.




• You should break up very long lines of code using ... (see an example below).




• When looping through a large matrix, always make the first index your inner loop. Matrix elements M(1,j) and M(2,j) are adjacent in memory, M(i,1) and M(i,2) are not. If you switch your two nested loops, your code will be faster for large arrays, as you use the cache better.

Now for the hard part: changing your logic.

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

Sum of difference between neighbours around element e:

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

Note that you compute the phase difference between two neighboring cells several times: for example, (a-b) will be computed when determining your value for e, but also in the previous loop iteration, when d was the center element; it will be computed two more times when that pair appears at the bottom of such a neighborhood.

Thus, you should pre-compute these differences and store them:

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

For this to work, phaseDifference must be rewritten, see below.

Now, vdiff has one fewer rows, and hdiff has one fewer columns, than M. We need to take this into account when indexing. Think of these as the elements in between the rows/columns of M.

The sum for one element at (i,j) now becomes:

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);

The full function now becomes:

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

Since Octave is slow with loops, it is likely that this vectorized version will be faster:

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);

I'm testing on MATLAB R2017a, and the vectorized form is actually slower than the version with loops, because indexing is still quite slow, whereas MATLAB has made huge strides in the last year speeding up code with loops.

Timings for an input of 2000x3000 elements:

• 
  checkPhase(M) (the OP's version) runs in 1.3s,


• and swapping the loop order it runs in 1.0s,


• 
  checkPhase2(M) runs in 0.17s,


• 
  checkPhase3(M) runs in 0.25s.

Your timings on Octave will be vastly different (much slower for the code with loops).

For phaseDifference to work with an array input, it cannot use if the way that the original function does. It is possible to rewrite the same logic using logical indexing to work on a full matrix, but this version is simpler:

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

Note I replaced the two inputs with a single one, as diff already computes the difference between neighbors.