Added pure-python implementation of scipy.ndimage.map_coordinates

This commit is contained in:
Luke Campagnola 2014-03-11 02:25:05 -04:00
parent 1eac666d02
commit 72a3902a18
3 changed files with 76 additions and 104 deletions

View File

@ -372,7 +372,7 @@ def affineSlice(data, shape, origin, vectors, axes, order=1, returnCoords=False,
"""
Take a slice of any orientation through an array. This is useful for extracting sections of multi-dimensional arrays such as MRI images for viewing as 1D or 2D data.
The slicing axes are aribtrary; they do not need to be orthogonal to the original data or even to each other. It is possible to use this function to extract arbitrary linear, rectangular, or parallelepiped shapes from within larger datasets. The original data is interpolated onto a new array of coordinates using scipy.ndimage.map_coordinates (see the scipy documentation for more information about this).
The slicing axes are aribtrary; they do not need to be orthogonal to the original data or even to each other. It is possible to use this function to extract arbitrary linear, rectangular, or parallelepiped shapes from within larger datasets. The original data is interpolated onto a new array of coordinates using scipy.ndimage.map_coordinates if it is available (see the scipy documentation for more information about this). If scipy is not available, then a slower implementation of map_coordinates is used.
For a graphical interface to this function, see :func:`ROI.getArrayRegion <pyqtgraph.ROI.getArrayRegion>`
@ -413,8 +413,9 @@ def affineSlice(data, shape, origin, vectors, axes, order=1, returnCoords=False,
"""
try:
import scipy.ndimage
have_scipy = True
except ImportError:
raise Exception("This function requires the scipy library, but it does not appear to be importable.")
have_scipy = False
# sanity check
if len(shape) != len(vectors):
@ -436,7 +437,6 @@ def affineSlice(data, shape, origin, vectors, axes, order=1, returnCoords=False,
#print "tr1:", tr1
## dims are now [(slice axes), (other axes)]
## make sure vectors are arrays
if not isinstance(vectors, np.ndarray):
vectors = np.array(vectors)
@ -456,8 +456,10 @@ def affineSlice(data, shape, origin, vectors, axes, order=1, returnCoords=False,
output = np.empty(tuple(shape) + extraShape, dtype=data.dtype)
for inds in np.ndindex(*extraShape):
ind = (Ellipsis,) + inds
#print data[ind].shape, x.shape, output[ind].shape, output.shape
output[ind] = scipy.ndimage.map_coordinates(data[ind], x, order=order, **kargs)
if have_scipy:
output[ind] = scipy.ndimage.map_coordinates(data[ind], x, order=order, **kargs)
else:
output[ind] = mapCoordinates(data[ind], x.T)
tr = list(range(output.ndim))
trb = []
@ -474,6 +476,54 @@ def affineSlice(data, shape, origin, vectors, axes, order=1, returnCoords=False,
else:
return output
def mapCoordinates(data, x, default=0.0):
"""
Pure-python alternative to scipy.ndimage.map_coordinates
*data* is an array of any shape.
*x* is an array with shape[-1] == data.ndim
Returns array of shape (x.shape[:-1] + data.shape)
"""
result = np.empty(x.shape[:-1] + data.shape, dtype=data.dtype)
nd = data.ndim
# First we generate arrays of indexes that are needed to
# extract the data surrounding each point
fields = np.mgrid[(slice(0,2),) * nd]
xmin = np.floor(x).astype(int)
xmax = xmin + 1
indexes = np.concatenate([xmin[np.newaxis, ...], xmax[np.newaxis, ...]])
fieldInds = []
totalMask = np.ones(x.shape[:-1], dtype=bool) # keep track of out-of-bound indexes
for ax in range(nd):
mask = (xmin[...,ax] >= 0) & (xmax[...,ax] < data.shape[ax])
totalMask &= mask
axisIndex = indexes[...,ax][fields[ax]]
axisMask = mask.astype(np.ubyte).reshape((1,)*(fields.ndim-1) + mask.shape)
axisIndex *= axisMask
fieldInds.append(axisIndex)
# Get data values surrounding each requested point
# fieldData[..., i] contains all 2**nd values needed to interpolate x[i]
fieldData = data[tuple(fieldInds)]
## Interpolate
s = np.empty((nd,) + fieldData.shape, dtype=float)
dx = x - xmin
# reshape fields for arithmetic against dx
for ax in range(nd):
f1 = fields[ax].reshape(fields[ax].shape + (1,)*(dx.ndim-1))
s[ax] = f1 * dx[:,ax] + (1-f1) * (1-dx[:,ax])
s = np.product(s, axis=0)
result = fieldData * s
for i in range(nd):
result = result.sum(axis=0)
result[~totalMask] = default
return result
def transformToArray(tr):
"""
Given a QTransform, return a 3x3 numpy array.

View File

@ -1083,105 +1083,6 @@ class ROI(GraphicsObject):
mapped = fn.transformCoordinates(img.transform(), coords)
return result, mapped
### transpose data so x and y are the first 2 axes
#trAx = range(0, data.ndim)
#trAx.remove(axes[0])
#trAx.remove(axes[1])
#tr1 = tuple(axes) + tuple(trAx)
#arr = data.transpose(tr1)
### Determine the minimal area of the data we will need
#(dataBounds, roiDataTransform) = self.getArraySlice(data, img, returnSlice=False, axes=axes)
### Pad data boundaries by 1px if possible
#dataBounds = (
#(max(dataBounds[0][0]-1, 0), min(dataBounds[0][1]+1, arr.shape[0])),
#(max(dataBounds[1][0]-1, 0), min(dataBounds[1][1]+1, arr.shape[1]))
#)
### Extract minimal data from array
#arr1 = arr[dataBounds[0][0]:dataBounds[0][1], dataBounds[1][0]:dataBounds[1][1]]
### Update roiDataTransform to reflect this extraction
#roiDataTransform *= QtGui.QTransform().translate(-dataBounds[0][0], -dataBounds[1][0])
#### (roiDataTransform now maps from ROI coords to extracted data coords)
### Rotate array
#if abs(self.state['angle']) > 1e-5:
#arr2 = ndimage.rotate(arr1, self.state['angle'] * 180 / np.pi, order=1)
### update data transforms to reflect this rotation
#rot = QtGui.QTransform().rotate(self.state['angle'] * 180 / np.pi)
#roiDataTransform *= rot
### The rotation also causes a shift which must be accounted for:
#dataBound = QtCore.QRectF(0, 0, arr1.shape[0], arr1.shape[1])
#rotBound = rot.mapRect(dataBound)
#roiDataTransform *= QtGui.QTransform().translate(-rotBound.left(), -rotBound.top())
#else:
#arr2 = arr1
#### Shift off partial pixels
## 1. map ROI into current data space
#roiBounds = roiDataTransform.mapRect(self.boundingRect())
## 2. Determine amount to shift data
#shift = (int(roiBounds.left()) - roiBounds.left(), int(roiBounds.bottom()) - roiBounds.bottom())
#if abs(shift[0]) > 1e-6 or abs(shift[1]) > 1e-6:
## 3. pad array with 0s before shifting
#arr2a = np.zeros((arr2.shape[0]+2, arr2.shape[1]+2) + arr2.shape[2:], dtype=arr2.dtype)
#arr2a[1:-1, 1:-1] = arr2
## 4. shift array and udpate transforms
#arr3 = ndimage.shift(arr2a, shift + (0,)*(arr2.ndim-2), order=1)
#roiDataTransform *= QtGui.QTransform().translate(1+shift[0], 1+shift[1])
#else:
#arr3 = arr2
#### Extract needed region from rotated/shifted array
## 1. map ROI into current data space (round these values off--they should be exact integer values at this point)
#roiBounds = roiDataTransform.mapRect(self.boundingRect())
##print self, roiBounds.height()
##import traceback
##traceback.print_stack()
#roiBounds = QtCore.QRect(round(roiBounds.left()), round(roiBounds.top()), round(roiBounds.width()), round(roiBounds.height()))
##2. intersect ROI with data bounds
#dataBounds = roiBounds.intersect(QtCore.QRect(0, 0, arr3.shape[0], arr3.shape[1]))
##3. Extract data from array
#db = dataBounds
#bounds = (
#(db.left(), db.right()+1),
#(db.top(), db.bottom()+1)
#)
#arr4 = arr3[bounds[0][0]:bounds[0][1], bounds[1][0]:bounds[1][1]]
#### Create zero array in size of ROI
#arr5 = np.zeros((roiBounds.width(), roiBounds.height()) + arr4.shape[2:], dtype=arr4.dtype)
### Fill array with ROI data
#orig = Point(dataBounds.topLeft() - roiBounds.topLeft())
#subArr = arr5[orig[0]:orig[0]+arr4.shape[0], orig[1]:orig[1]+arr4.shape[1]]
#subArr[:] = arr4[:subArr.shape[0], :subArr.shape[1]]
### figure out the reverse transpose order
#tr2 = np.array(tr1)
#for i in range(0, len(tr2)):
#tr2[tr1[i]] = i
#tr2 = tuple(tr2)
### Untranspose array before returning
#return arr5.transpose(tr2)
def getAffineSliceParams(self, data, img, axes=(0,1)):
"""
Returns the parameters needed to use :func:`affineSlice <pyqtgraph.affineSlice>` to

View File

@ -19,3 +19,24 @@ def testSolve3D():
tr2 = pg.solve3DTransform(p1, p2)
assert_array_almost_equal(tr[:3], tr2[:3])
def test_mapCoordinates():
data = np.array([[ 0., 1., 2.],
[ 2., 3., 5.],
[ 7., 7., 4.]])
x = np.array([[ 0.3, 0.6],
[ 1. , 1. ],
[ 0.5, 1. ],
[ 0.5, 2.5],
[ 10. , 10. ]])
result = pg.mapCoordinates(data, x)
import scipy.ndimage
spresult = scipy.ndimage.map_coordinates(data, x.T, order=1)
assert_array_almost_equal(result, spresult)