[Numpy-discussion] Detecting phase windings

gruben@bigpond.n... gruben@bigpond.n...
Mon Jun 16 11:30:24 CDT 2008


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?

Gary Ruben


More information about the Numpy-discussion mailing list