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
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
