[Numpy-discussion] Detecting phase windings

Anne Archibald peridot.faceted at gmail.com
Fri Jun 20 00:13:07 EDT 2008


2008/6/16 gruben at bigpond.net.au <gruben at 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



More information about the NumPy-Discussion mailing list