From 7f0556b05f98b44e87d0bb14d91a8385ffd8a01f Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Tue, 13 Sep 2016 18:08:11 -0700 Subject: [PATCH 1/5] Implement order=0 in functions.interpolateArray; use scipy only for order>1. --- pyqtgraph/functions.py | 131 ++++++++++++++++++++++++----------------- 1 file changed, 76 insertions(+), 55 deletions(-) diff --git a/pyqtgraph/functions.py b/pyqtgraph/functions.py index d79c350f..187a4717 100644 --- a/pyqtgraph/functions.py +++ b/pyqtgraph/functions.py @@ -409,12 +409,45 @@ def eq(a, b): else: raise Exception("== operator returned type %s" % str(type(e))) + +def affineSliceCoords(shape, origin, vectors, axes): + """Return the array of coordinates used to sample data arrays in affineSlice(). + """ + # sanity check + if len(shape) != len(vectors): + raise Exception("shape and vectors must have same length.") + if len(origin) != len(axes): + raise Exception("origin and axes must have same length.") + for v in vectors: + if len(v) != len(axes): + raise Exception("each vector must be same length as axes.") + + shape = list(map(np.ceil, shape)) + + ## make sure vectors are arrays + if not isinstance(vectors, np.ndarray): + vectors = np.array(vectors) + if not isinstance(origin, np.ndarray): + origin = np.array(origin) + origin.shape = (len(axes),) + (1,)*len(shape) + + ## Build array of sample locations. + grid = np.mgrid[tuple([slice(0,x) for x in shape])] ## mesh grid of indexes + x = (grid[np.newaxis,...] * vectors.transpose()[(Ellipsis,) + (np.newaxis,)*len(shape)]).sum(axis=1) ## magic + x += origin + + return x + def affineSlice(data, shape, origin, vectors, axes, order=1, returnCoords=False, **kargs): """ - 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. + 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 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. + 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 either interpolateArray if order<2 + or scipy.ndimage.map_coordinates otherwise. For a graphical interface to this function, see :func:`ROI.getArrayRegion ` @@ -453,47 +486,24 @@ def affineSlice(data, shape, origin, vectors, axes, order=1, returnCoords=False, affineSlice(data, shape=(20,20), origin=(40,0,0), vectors=((-1, 1, 0), (-1, 0, 1)), axes=(1,2,3)) """ - try: - import scipy.ndimage - have_scipy = True - except ImportError: - have_scipy = False - have_scipy = False - - # sanity check - if len(shape) != len(vectors): - raise Exception("shape and vectors must have same length.") - if len(origin) != len(axes): - raise Exception("origin and axes must have same length.") - for v in vectors: - if len(v) != len(axes): - raise Exception("each vector must be same length as axes.") - - shape = list(map(np.ceil, shape)) + x = affineSliceCoords(shape, origin, vectors, axes) ## transpose data so slice axes come first trAx = list(range(data.ndim)) - for x in axes: - trAx.remove(x) + for ax in axes: + trAx.remove(ax) tr1 = tuple(axes) + tuple(trAx) data = data.transpose(tr1) #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) - if not isinstance(origin, np.ndarray): - origin = np.array(origin) - origin.shape = (len(axes),) + (1,)*len(shape) - - ## Build array of sample locations. - grid = np.mgrid[tuple([slice(0,x) for x in shape])] ## mesh grid of indexes - x = (grid[np.newaxis,...] * vectors.transpose()[(Ellipsis,) + (np.newaxis,)*len(shape)]).sum(axis=1) ## magic - x += origin - ## iterate manually over unused axes since map_coordinates won't do it for us - if have_scipy: + if order > 1: + try: + import scipy.ndimage + except ImportError: + raise ImportError("Interpolating with order > 1 requires the scipy.ndimage module, but it could not be imported.") + + # iterate manually over unused axes since map_coordinates won't do it for us extraShape = data.shape[len(axes):] output = np.empty(tuple(shape) + extraShape, dtype=data.dtype) for inds in np.ndindex(*extraShape): @@ -502,8 +512,8 @@ def affineSlice(data, shape, origin, vectors, axes, order=1, returnCoords=False, else: # map_coordinates expects the indexes as the first axis, whereas # interpolateArray expects indexes at the last axis. - tr = tuple(range(1,x.ndim)) + (0,) - output = interpolateArray(data, x.transpose(tr)) + tr = tuple(range(1, x.ndim)) + (0,) + output = interpolateArray(data, x.transpose(tr), order=order) tr = list(range(output.ndim)) trb = [] @@ -520,16 +530,21 @@ def affineSlice(data, shape, origin, vectors, axes, order=1, returnCoords=False, else: return output -def interpolateArray(data, x, default=0.0): + +def interpolateArray(data, x, default=0.0, order=1): """ N-dimensional interpolation similar to scipy.ndimage.map_coordinates. This function returns linearly-interpolated values sampled from a regular grid of data. - *data* is an array of any shape containing the values to be interpolated. - *x* is an array with (shape[-1] <= data.ndim) containing the locations - within *data* to interpolate. + ============== =========================================================================================== + **Arguments:** + *data* Array of any shape containing the values to be interpolated. + *x* Array with (shape[-1] <= data.ndim) containing the locations within *data* to interpolate. + *default* Value to return for locations in *x* that are outside the bounds of *data*. + *order* Order of interpolation: 0=nearest, 1=linear. + ============== =========================================================================================== Returns array of shape (x.shape[:-1] + data.shape[x.shape[-1]:]) @@ -574,8 +589,11 @@ def interpolateArray(data, x, default=0.0): This is useful for interpolating from arrays of colors, vertexes, etc. """ + if order not in (0, 1): + raise ValueError("interpolateArray requires order=0 or 1 (got %s)" % order) + prof = debug.Profiler() - + nd = data.ndim md = x.shape[-1] if md > nd: @@ -583,7 +601,7 @@ def interpolateArray(data, x, default=0.0): # First we generate arrays of indexes that are needed to # extract the data surrounding each point - fields = np.mgrid[(slice(0,2),) * md] + fields = np.mgrid[(slice(0,order+1),) * md] xmin = np.floor(x).astype(int) xmax = xmin + 1 indexes = np.concatenate([xmin[np.newaxis, ...], xmax[np.newaxis, ...]]) @@ -609,18 +627,21 @@ def interpolateArray(data, x, default=0.0): prof() ## Interpolate - s = np.empty((md,) + fieldData.shape, dtype=float) - dx = x - xmin - # reshape fields for arithmetic against dx - for ax in range(md): - f1 = fields[ax].reshape(fields[ax].shape + (1,)*(dx.ndim-1)) - sax = f1 * dx[...,ax] + (1-f1) * (1-dx[...,ax]) - sax = sax.reshape(sax.shape + (1,) * (s.ndim-1-sax.ndim)) - s[ax] = sax - s = np.product(s, axis=0) - result = fieldData * s - for i in range(md): - result = result.sum(axis=0) + if order == 0: + result = fieldData[0,0] + else: + s = np.empty((md,) + fieldData.shape, dtype=float) + dx = x - xmin + # reshape fields for arithmetic against dx + for ax in range(md): + f1 = fields[ax].reshape(fields[ax].shape + (1,)*(dx.ndim-1)) + sax = f1 * dx[...,ax] + (1-f1) * (1-dx[...,ax]) + sax = sax.reshape(sax.shape + (1,) * (s.ndim-1-sax.ndim)) + s[ax] = sax + s = np.product(s, axis=0) + result = fieldData * s + for i in range(md): + result = result.sum(axis=0) prof() From be07979b3904b84e584e992b7b545d654f7c1c58 Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Fri, 16 Sep 2016 17:16:16 -0700 Subject: [PATCH 2/5] Add returnMappedCoords option to LineSegmentROI.getArrayRegion --- pyqtgraph/graphicsItems/ROI.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/pyqtgraph/graphicsItems/ROI.py b/pyqtgraph/graphicsItems/ROI.py index 81a4e651..963ecb05 100644 --- a/pyqtgraph/graphicsItems/ROI.py +++ b/pyqtgraph/graphicsItems/ROI.py @@ -2070,9 +2070,9 @@ class LineSegmentROI(ROI): if len(positions) > 2: raise Exception("LineSegmentROI must be defined by exactly 2 positions. For more points, use PolyLineROI.") + self.endpoints = [] for i, p in enumerate(positions): - self.addFreeHandle(p, item=handles[i]) - + self.endpoints.append(self.addFreeHandle(p, item=handles[i])) def listPoints(self): return [p['item'].pos() for p in self.handles] @@ -2080,8 +2080,8 @@ class LineSegmentROI(ROI): def paint(self, p, *args): p.setRenderHint(QtGui.QPainter.Antialiasing) p.setPen(self.currentPen) - h1 = self.handles[0]['item'].pos() - h2 = self.handles[1]['item'].pos() + h1 = self.endpoints[0].pos() + h2 = self.endpoints[1].pos() p.drawLine(h1, h2) def boundingRect(self): @@ -2090,8 +2090,8 @@ class LineSegmentROI(ROI): def shape(self): p = QtGui.QPainterPath() - h1 = self.handles[0]['item'].pos() - h2 = self.handles[1]['item'].pos() + h1 = self.endpoints[0].pos() + h2 = self.endpoints[1].pos() dh = h2-h1 if dh.length() == 0: return p @@ -2109,7 +2109,7 @@ class LineSegmentROI(ROI): return p - def getArrayRegion(self, data, img, axes=(0,1), order=1, **kwds): + def getArrayRegion(self, data, img, axes=(0,1), order=1, returnMappedCoords=False, **kwds): """ Use the position of this ROI relative to an imageItem to pull a slice from an array. @@ -2120,15 +2120,15 @@ class LineSegmentROI(ROI): See ROI.getArrayRegion() for a description of the arguments. """ - imgPts = [self.mapToItem(img, h['item'].pos()) for h in self.handles] + imgPts = [self.mapToItem(img, h.pos()) for h in self.endpoints] rgns = [] - for i in range(len(imgPts)-1): - d = Point(imgPts[i+1] - imgPts[i]) - o = Point(imgPts[i]) - r = fn.affineSlice(data, shape=(int(d.length()),), vectors=[Point(d.norm())], origin=o, axes=axes, order=order, **kwds) - rgns.append(r) - - return np.concatenate(rgns, axis=axes[0]) + coords = [] + + d = Point(imgPts[1] - imgPts[0]) + o = Point(imgPts[0]) + rgn = fn.affineSlice(data, shape=(int(d.length()),), vectors=[Point(d.norm())], origin=o, axes=axes, order=order, returnCoords=returnMappedCoords, **kwds) + + return rgn class _PolyLineSegment(LineSegmentROI): From 92fc9dbe2f0e97f4cac04454877ea65232b1a19b Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Wed, 12 Oct 2016 09:58:03 -0700 Subject: [PATCH 3/5] Add unit test for interpolateArray with order=0 docstring update --- pyqtgraph/functions.py | 7 ++++-- pyqtgraph/tests/test_functions.py | 37 ++++++++++++++----------------- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/pyqtgraph/functions.py b/pyqtgraph/functions.py index 187a4717..32d9f2bf 100644 --- a/pyqtgraph/functions.py +++ b/pyqtgraph/functions.py @@ -536,12 +536,15 @@ def interpolateArray(data, x, default=0.0, order=1): N-dimensional interpolation similar to scipy.ndimage.map_coordinates. This function returns linearly-interpolated values sampled from a regular - grid of data. + grid of data. It differs from `ndimage.map_coordinates` by allowing broadcasting + within the input array. ============== =========================================================================================== **Arguments:** *data* Array of any shape containing the values to be interpolated. - *x* Array with (shape[-1] <= data.ndim) containing the locations within *data* to interpolate. + *x* Array with (shape[-1] <= data.ndim) containing the locations within *data* to interpolate. + (note: the axes for this argument are transposed relative to the same argument for + `ndimage.map_coordinates`). *default* Value to return for locations in *x* that are outside the bounds of *data*. *order* Order of interpolation: 0=nearest, 1=linear. ============== =========================================================================================== diff --git a/pyqtgraph/tests/test_functions.py b/pyqtgraph/tests/test_functions.py index bfa7e0ea..4c9cabfe 100644 --- a/pyqtgraph/tests/test_functions.py +++ b/pyqtgraph/tests/test_functions.py @@ -22,9 +22,17 @@ def testSolve3D(): assert_array_almost_equal(tr[:3], tr2[:3]) -def test_interpolateArray(): +def test_interpolateArray_order0(): + check_interpolateArray(order=0) + + +def test_interpolateArray_order1(): + check_interpolateArray(order=1) + + +def check_interpolateArray(order): def interpolateArray(data, x): - result = pg.interpolateArray(data, x) + result = pg.interpolateArray(data, x, order=order) assert result.shape == x.shape[:-1] + data.shape[x.shape[-1]:] return result @@ -48,7 +56,6 @@ def test_interpolateArray(): with pytest.raises(TypeError): interpolateArray(data, np.ones((5, 5, 3,))) - x = np.array([[ 0.3, 0.6], [ 1. , 1. ], [ 0.5, 1. ], @@ -56,9 +63,10 @@ def test_interpolateArray(): [ 10. , 10. ]]) result = interpolateArray(data, x) - #import scipy.ndimage - #spresult = scipy.ndimage.map_coordinates(data, x.T, order=1) - spresult = np.array([ 5.92, 20. , 11. , 0. , 0. ]) # generated with the above line + # make sure results match ndimage.map_coordinates + import scipy.ndimage + spresult = scipy.ndimage.map_coordinates(data, x.T, order=order) + #spresult = np.array([ 5.92, 20. , 11. , 0. , 0. ]) # generated with the above line assert_array_almost_equal(result, spresult) @@ -78,24 +86,13 @@ def test_interpolateArray(): [[1.5, 0.5], [1.5, 1.0], [1.5, 1.5]]]) r1 = interpolateArray(data, x) - #r2 = scipy.ndimage.map_coordinates(data, x.transpose(2,0,1), order=1) - r2 = np.array([[ 8.25, 11. , 16.5 ], # generated with the above line - [ 82.5 , 110. , 165. ]]) + r2 = scipy.ndimage.map_coordinates(data, x.transpose(2,0,1), order=order) + #r2 = np.array([[ 8.25, 11. , 16.5 ], # generated with the above line + #[ 82.5 , 110. , 165. ]]) assert_array_almost_equal(r1, r2) - # test interpolate where data.ndim > x.shape[1] - - data = np.array([[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]]) # 2x2x3 - x = np.array([[1, 1], [0, 0.5], [5, 5]]) - - r1 = interpolateArray(data, x) - assert np.all(r1[0] == data[1, 1]) - assert np.all(r1[1] == 0.5 * (data[0, 0] + data[0, 1])) - assert np.all(r1[2] == 0) - - def test_subArray(): a = np.array([0, 0, 111, 112, 113, 0, 121, 122, 123, 0, 0, 0, 211, 212, 213, 0, 221, 222, 223, 0, 0, 0, 0]) b = pg.subArray(a, offset=2, shape=(2,2,3), stride=(10,4,1)) From e35f59fcb77d50ab493d7c9db847c24c4fbd6006 Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Wed, 12 Oct 2016 10:26:54 -0700 Subject: [PATCH 4/5] Fix interpolateArray for order=0 --- pyqtgraph/functions.py | 65 +++++++++++++++++-------------- pyqtgraph/tests/test_functions.py | 8 ++-- 2 files changed, 40 insertions(+), 33 deletions(-) diff --git a/pyqtgraph/functions.py b/pyqtgraph/functions.py index 32d9f2bf..8593241e 100644 --- a/pyqtgraph/functions.py +++ b/pyqtgraph/functions.py @@ -602,37 +602,44 @@ def interpolateArray(data, x, default=0.0, order=1): if md > nd: raise TypeError("x.shape[-1] must be less than or equal to data.ndim") - # First we generate arrays of indexes that are needed to - # extract the data surrounding each point - fields = np.mgrid[(slice(0,order+1),) * md] - 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(md): - mask = (xmin[...,ax] >= 0) & (x[...,ax] <= data.shape[ax]-1) - # keep track of points that need to be set to default - totalMask &= mask - - # ..and keep track of indexes that are out of bounds - # (note that when x[...,ax] == data.shape[ax], then xmax[...,ax] will be out - # of bounds, but the interpolation will work anyway) - mask &= (xmax[...,ax] < data.shape[ax]) - axisIndex = indexes[...,ax][fields[ax]] - axisIndex[axisIndex < 0] = 0 - axisIndex[axisIndex >= data.shape[ax]] = 0 - fieldInds.append(axisIndex) - prof() - - # Get data values surrounding each requested point - fieldData = data[tuple(fieldInds)] - prof() - - ## Interpolate if order == 0: - result = fieldData[0,0] - else: + xinds = np.round(x).astype(int) # NOTE: for 0.5 this rounds to the nearest *even* number + for ax in range(md): + mask = (xinds[...,ax] >= 0) & (xinds[...,ax] <= data.shape[ax]-1) + xinds[...,ax][~mask] = 0 + # keep track of points that need to be set to default + totalMask &= mask + result = data[tuple([xinds[...,i] for i in range(xinds.shape[-1])])] + + elif order == 1: + # First we generate arrays of indexes that are needed to + # extract the data surrounding each point + fields = np.mgrid[(slice(0,order+1),) * md] + xmin = np.floor(x).astype(int) + xmax = xmin + 1 + indexes = np.concatenate([xmin[np.newaxis, ...], xmax[np.newaxis, ...]]) + fieldInds = [] + for ax in range(md): + mask = (xmin[...,ax] >= 0) & (x[...,ax] <= data.shape[ax]-1) + # keep track of points that need to be set to default + totalMask &= mask + + # ..and keep track of indexes that are out of bounds + # (note that when x[...,ax] == data.shape[ax], then xmax[...,ax] will be out + # of bounds, but the interpolation will work anyway) + mask &= (xmax[...,ax] < data.shape[ax]) + axisIndex = indexes[...,ax][fields[ax]] + axisIndex[axisIndex < 0] = 0 + axisIndex[axisIndex >= data.shape[ax]] = 0 + fieldInds.append(axisIndex) + prof() + + # Get data values surrounding each requested point + fieldData = data[tuple(fieldInds)] + prof() + + ## Interpolate s = np.empty((md,) + fieldData.shape, dtype=float) dx = x - xmin # reshape fields for arithmetic against dx diff --git a/pyqtgraph/tests/test_functions.py b/pyqtgraph/tests/test_functions.py index 4c9cabfe..7ad3bf91 100644 --- a/pyqtgraph/tests/test_functions.py +++ b/pyqtgraph/tests/test_functions.py @@ -58,8 +58,8 @@ def check_interpolateArray(order): x = np.array([[ 0.3, 0.6], [ 1. , 1. ], - [ 0.5, 1. ], - [ 0.5, 2.5], + [ 0.501, 1. ], # NOTE: testing at exactly 0.5 can yield different results from map_coordinates + [ 0.501, 2.501], # due to differences in rounding [ 10. , 10. ]]) result = interpolateArray(data, x) @@ -82,8 +82,8 @@ def check_interpolateArray(order): # test mapping 2D array of locations - x = np.array([[[0.5, 0.5], [0.5, 1.0], [0.5, 1.5]], - [[1.5, 0.5], [1.5, 1.0], [1.5, 1.5]]]) + x = np.array([[[0.501, 0.501], [0.501, 1.0], [0.501, 1.501]], + [[1.501, 0.501], [1.501, 1.0], [1.501, 1.501]]]) r1 = interpolateArray(data, x) r2 = scipy.ndimage.map_coordinates(data, x.transpose(2,0,1), order=order) From 4e7773fa0bb90e733976093259a7a69eb5631edb Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Wed, 14 Dec 2016 11:02:40 -0800 Subject: [PATCH 5/5] Add scipy to travis requirements--some unit tests require this --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 2c7b7769..c4a67ac3 100644 --- a/.travis.yml +++ b/.travis.yml @@ -51,7 +51,7 @@ install: - conda update conda --yes - conda create -n test_env python=${PYTHON} --yes - source activate test_env - - conda install numpy pyopengl pytest flake8 six coverage --yes + - conda install numpy scipy pyopengl pytest flake8 six coverage --yes - echo ${QT} - echo ${TEST} - echo ${PYTHON}