2008/6/16 gruben@bigpond.net.au <gruben@bigpond.net.au>:

I have a speed problem with the approach I'm using to detect phase wrappings in a 3D data set. In my application, phaseField is a 3D array containing the phase values of a field. In order to detect the vortices/phase windings at each point, I check for windings on each of 3 faces of a 2x2 cube with the following code:

def HasVortex(cube): ''' cube is a 2x2x2 array containing the phase values returns True if any of 3 orthogonal faces contains a phase wrapping ''' # extract 3 cube faces - flatten to make indexing easier c1 = cube[0,:,:].flatten() c3 = cube[:,0,:].flatten() c5 = cube[:,:,0].flatten()

# now order the phases in a circuit around each face, finishing in the same corner as we began # Use numpy's phase unwrapping code, which has a default threshold of pi u1 = unwrap(c1[cwi]) u3 = unwrap(c3[cwi]) u5 = unwrap(c5[cwi])

# check whether the final unwrapped value coincides with the value in the cell we started at return not allclose(r_[u1[0],u3[0],u5[0]], r_[u1[4],u3[4],u5[4]])

vortexArray = array([int(HasVortex(phaseField[x:x+2,y:y+2,z:z+2])) for x in range(phaseField.shape[0]-1) for y in range(phaseField.shape[1]-1) for z in range(phaseField.shape[2]-1)]\ ).reshape((xs-1,ys-1,zs-1))

Whilst this does what I want, it's incredibly slow. I'm wondering whether anyone has done this before, or can think of a better approach.

My thoughts about a better approach are to stack the values along 3 new dimensions making a 6D array and use unwrap along the 3 new dimensions. Something like the following may give you an idea of what I mean, but it's a toy example trying to extend 2D into 3D, rather than 3D into 6D, because I haven't come to grips with enough of numpy's axis reshaping and stacking tricks to avoid getting a severe headache when trying to generalize it:

a = array([[1.,3.], [6.,5.]]) b = np.dstack((a, roll(a,-1,axis=1), roll(roll(a,-1,axis=0),-1,axis=1), roll(a,-1,axis=0), a)) print np.unwrap(b, axis=2)

A problem with this approach is that the memory requirements for the stacked version will be much greater, so some version using views would be preferable.

Any ideas?

I'd start by looking at the problem one face at a time: def find_vortices(X, axis=0): XX = np.rollaxis(X,axis) loop = np.concatenate((XX[np.newaxis,:-1,:-1,...], XX[np.newaxis,1:,:-1,...], XX[np.newaxis,1:,1:,...], XX[np.newaxis,:-1,1:,...], XX[np.newaxis,:-1,:-1,...]),axis=0) loop = np.unwrap(loop) r = np.abs(loop[0,...]-loop[-1,...])<np.pi/2 return np.rollaxis(r,0,axis+1) This does require temporary space some five or ten times your original space required, which is not good. You might be able to do better by reimplementing the phase wrapping code; your condition is equivalent to "the number of upward steps of more than pi is different from the number of downward steps of more than pi". But every operation on the whole array at once is going to produce a temporary roughly the size of the array, as in fact is the result boolean array (though you will presumably use nonzero() to convert it to a list of locations, which will presumably be fairly small until you hit a totally decorrelated image). If your array is large, as sounds like the case, you can try a standard trick for cutting down on temporary sizes: use a for loop along the smallest dimension and a vector along the rest. In fact in this case I'd use a 2D vortex finder and iterate along the remaining axis; do this three times and you're done. (This can bit you if you have an extremely long skinny array, but otherwise the python overhead will be pretty small.) Anne