[Numpy-discussion] Detecting phase windings

[Gary Ruben]

I had a chance to look at Anne's suggestion from this thread
<http://www.mail-archive.com/numpy-discussion@scipy.org/msg10091.html>
and I thought I should post my phase winding finder solution, which is 
slightly modified from her idea. Thanks Anne. This is a vast improvement 
over my original slow code, and is useful to me now, but I will probably 
have to rewrite it in C, weave or Cython when I start generating large 
data sets.

import numpy as np
from pyvtk import *

def find_vortices(x, axis=0):
     xx = np.rollaxis(x, axis)
     r = np.empty_like(xx).astype(np.bool)
     for i in range(xx.shape[0]):
         print i,
         xxx = xx[i,...]
         loop = np.concatenate(([xxx],
                                [np.roll(xxx,1,0)],
                                [np.roll(np.roll(xxx,1,0),1,1)],
                                [np.roll(xxx,1,1)],
                                [xxx]), axis=0)
         loop = np.unwrap(loop, axis=0)
         r[i,...] = np.abs(loop[0,...]-loop[-1,...])>pi/2

     return np.rollaxis(r, 0, axis+1)[1:-1,1:-1,1:-1]

and call it like so on the 3D phaseField array, which is a float32 array 
containing the phase angle at each point:

# Detect the nodal lines
b0 = find_vortices(phaseField, axis=0)
b0 |= find_vortices(phaseField, axis=1)
b0 |= find_vortices(phaseField, axis=2)

# output vortices to vtk
indices = np.transpose(np.nonzero(b0)).tolist()
vtk = VtkData(UnstructuredGrid(indices))
vtk.tofile('%s_vol'%sys.argv[0][:-3],'binary')
del vtk

--
Gary R.