From 816069c020430c5fd1356ef4357ac9ae82f71753 Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Sat, 23 Nov 2013 14:35:09 -0500 Subject: [PATCH 1/9] All scipy imports in the library are now optional (need to test each of these changes) Several examples still require scipy. --- examples/ViewBox.py | 1 - pyqtgraph/SRTTransform3D.py | 4 +- pyqtgraph/colormap.py | 6 +- pyqtgraph/flowchart/library/Filters.py | 19 +++-- pyqtgraph/flowchart/library/functions.py | 21 ++++- pyqtgraph/functions.py | 100 +++++------------------ pyqtgraph/graphicsItems/PlotDataItem.py | 5 +- pyqtgraph/graphicsItems/ROI.py | 5 +- 8 files changed, 68 insertions(+), 93 deletions(-) diff --git a/examples/ViewBox.py b/examples/ViewBox.py index 2dcbb758..3a66afe3 100644 --- a/examples/ViewBox.py +++ b/examples/ViewBox.py @@ -14,7 +14,6 @@ import initExample ## This example uses a ViewBox to create a PlotWidget-like interface -#from scipy import random import numpy as np from pyqtgraph.Qt import QtGui, QtCore import pyqtgraph as pg diff --git a/pyqtgraph/SRTTransform3D.py b/pyqtgraph/SRTTransform3D.py index 417190e1..7da3c1bd 100644 --- a/pyqtgraph/SRTTransform3D.py +++ b/pyqtgraph/SRTTransform3D.py @@ -4,7 +4,6 @@ from .Vector import Vector from .Transform3D import Transform3D from .Vector import Vector import numpy as np -import scipy.linalg class SRTTransform3D(Transform3D): """4x4 Transform matrix that can always be represented as a combination of 3 matrices: scale * rotate * translate @@ -118,6 +117,7 @@ class SRTTransform3D(Transform3D): The input matrix must be affine AND have no shear, otherwise the conversion will most likely fail. """ + import numpy.linalg for i in range(4): self.setRow(i, m.row(i)) m = self.matrix().reshape(4,4) @@ -134,7 +134,7 @@ class SRTTransform3D(Transform3D): ## rotation axis is the eigenvector with eigenvalue=1 r = m[:3, :3] / scale[:, np.newaxis] try: - evals, evecs = scipy.linalg.eig(r) + evals, evecs = numpy.linalg.eig(r) except: print("Rotation matrix: %s" % str(r)) print("Scale: %s" % str(scale)) diff --git a/pyqtgraph/colormap.py b/pyqtgraph/colormap.py index 38c12097..559e6688 100644 --- a/pyqtgraph/colormap.py +++ b/pyqtgraph/colormap.py @@ -1,5 +1,4 @@ import numpy as np -import scipy.interpolate from .Qt import QtGui, QtCore class ColorMap(object): @@ -84,6 +83,11 @@ class ColorMap(object): qcolor Values are returned as an array of QColor objects. =========== =============================================================== """ + try: + import scipy.interpolate + except: + raise Exception("Colormap.map() requires the package scipy.interpolate, but it could not be imported.") + if isinstance(mode, basestring): mode = self.enumMap[mode.lower()] diff --git a/pyqtgraph/flowchart/library/Filters.py b/pyqtgraph/flowchart/library/Filters.py index 518c8c18..e5bb2453 100644 --- a/pyqtgraph/flowchart/library/Filters.py +++ b/pyqtgraph/flowchart/library/Filters.py @@ -1,9 +1,6 @@ # -*- coding: utf-8 -*- from ...Qt import QtCore, QtGui from ..Node import Node -from scipy.signal import detrend -from scipy.ndimage import median_filter, gaussian_filter -#from ...SignalProxy import SignalProxy from . import functions from .common import * import numpy as np @@ -119,7 +116,11 @@ class Median(CtrlNode): @metaArrayWrapper def processData(self, data): - return median_filter(data, self.ctrls['n'].value()) + try: + import scipy.ndimage + except ImportError: + raise Exception("MedianFilter node requires the package scipy.ndimage.") + return scipy.ndimage.median_filter(data, self.ctrls['n'].value()) class Mode(CtrlNode): """Filters data by taking the mode (histogram-based) of a sliding window""" @@ -156,7 +157,11 @@ class Gaussian(CtrlNode): @metaArrayWrapper def processData(self, data): - return gaussian_filter(data, self.ctrls['sigma'].value()) + try: + import scipy.ndimage + except ImportError: + raise Exception("GaussianFilter node requires the package scipy.ndimage.") + return scipy.ndimage.gaussian_filter(data, self.ctrls['sigma'].value()) class Derivative(CtrlNode): @@ -189,6 +194,10 @@ class Detrend(CtrlNode): @metaArrayWrapper def processData(self, data): + try: + from scipy.signal import detrend + except ImportError: + raise Exception("DetrendFilter node requires the package scipy.signal.") return detrend(data) diff --git a/pyqtgraph/flowchart/library/functions.py b/pyqtgraph/flowchart/library/functions.py index 9efb8f36..027e1386 100644 --- a/pyqtgraph/flowchart/library/functions.py +++ b/pyqtgraph/flowchart/library/functions.py @@ -1,4 +1,3 @@ -import scipy import numpy as np from ...metaarray import MetaArray @@ -47,6 +46,11 @@ def downsample(data, n, axis=0, xvals='subsample'): def applyFilter(data, b, a, padding=100, bidir=True): """Apply a linear filter with coefficients a, b. Optionally pad the data before filtering and/or run the filter in both directions.""" + try: + import scipy.signal + except ImportError: + raise Exception("applyFilter() requires the package scipy.signal.") + d1 = data.view(np.ndarray) if padding > 0: @@ -67,6 +71,11 @@ def applyFilter(data, b, a, padding=100, bidir=True): def besselFilter(data, cutoff, order=1, dt=None, btype='low', bidir=True): """return data passed through bessel filter""" + try: + import scipy.signal + except ImportError: + raise Exception("besselFilter() requires the package scipy.signal.") + if dt is None: try: tvals = data.xvals('Time') @@ -85,6 +94,11 @@ def besselFilter(data, cutoff, order=1, dt=None, btype='low', bidir=True): def butterworthFilter(data, wPass, wStop=None, gPass=2.0, gStop=20.0, order=1, dt=None, btype='low', bidir=True): """return data passed through bessel filter""" + try: + import scipy.signal + except ImportError: + raise Exception("butterworthFilter() requires the package scipy.signal.") + if dt is None: try: tvals = data.xvals('Time') @@ -175,6 +189,11 @@ def denoise(data, radius=2, threshold=4): def adaptiveDetrend(data, x=None, threshold=3.0): """Return the signal with baseline removed. Discards outliers from baseline measurement.""" + try: + import scipy.signal + except ImportError: + raise Exception("adaptiveDetrend() requires the package scipy.signal.") + if x is None: x = data.xvals(0) diff --git a/pyqtgraph/functions.py b/pyqtgraph/functions.py index 427fb01d..62df1ce3 100644 --- a/pyqtgraph/functions.py +++ b/pyqtgraph/functions.py @@ -34,17 +34,6 @@ import decimal, re import ctypes import sys, struct -try: - import scipy.ndimage - HAVE_SCIPY = True - if getConfigOption('useWeave'): - try: - import scipy.weave - except ImportError: - setConfigOptions(useWeave=False) -except ImportError: - HAVE_SCIPY = False - from . import debug def siScale(x, minVal=1e-25, allowUnicode=True): @@ -422,7 +411,9 @@ 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)) """ - if not HAVE_SCIPY: + try: + import scipy.ndimage + except ImportError: raise Exception("This function requires the scipy library, but it does not appear to be importable.") # sanity check @@ -579,15 +570,14 @@ def solve3DTransform(points1, points2): Find a 3D transformation matrix that maps points1 onto points2. Points must be specified as a list of 4 Vectors. """ - if not HAVE_SCIPY: - raise Exception("This function depends on the scipy library, but it does not appear to be importable.") + import numpy.linalg A = np.array([[points1[i].x(), points1[i].y(), points1[i].z(), 1] for i in range(4)]) B = np.array([[points2[i].x(), points2[i].y(), points2[i].z(), 1] for i in range(4)]) ## solve 3 sets of linear equations to determine transformation matrix elements matrix = np.zeros((4,4)) for i in range(3): - matrix[i] = scipy.linalg.solve(A, B[:,i]) ## solve Ax = B; x is one row of the desired transformation matrix + matrix[i] = numpy.linalg.solve(A, B[:,i]) ## solve Ax = B; x is one row of the desired transformation matrix return matrix @@ -600,8 +590,7 @@ def solveBilinearTransform(points1, points2): mapped = np.dot(matrix, [x*y, x, y, 1]) """ - if not HAVE_SCIPY: - raise Exception("This function depends on the scipy library, but it does not appear to be importable.") + import numpy.linalg ## A is 4 rows (points) x 4 columns (xy, x, y, 1) ## B is 4 rows (points) x 2 columns (x, y) A = np.array([[points1[i].x()*points1[i].y(), points1[i].x(), points1[i].y(), 1] for i in range(4)]) @@ -610,7 +599,7 @@ def solveBilinearTransform(points1, points2): ## solve 2 sets of linear equations to determine transformation matrix elements matrix = np.zeros((2,4)) for i in range(2): - matrix[i] = scipy.linalg.solve(A, B[:,i]) ## solve Ax = B; x is one row of the desired transformation matrix + matrix[i] = numpy.linalg.solve(A, B[:,i]) ## solve Ax = B; x is one row of the desired transformation matrix return matrix @@ -629,6 +618,10 @@ def rescaleData(data, scale, offset, dtype=None): try: if not getConfigOption('useWeave'): raise Exception('Weave is disabled; falling back to slower version.') + try: + import scipy.weave + except ImportError: + raise Exception('scipy.weave is not importable; falling back to slower version.') ## require native dtype when using weave if not data.dtype.isnative: @@ -671,68 +664,13 @@ def applyLookupTable(data, lut): Uses values in *data* as indexes to select values from *lut*. The returned data has shape data.shape + lut.shape[1:] - Uses scipy.weave to improve performance if it is available. Note: color gradient lookup tables can be generated using GradientWidget. """ if data.dtype.kind not in ('i', 'u'): data = data.astype(int) - ## using np.take appears to be faster than even the scipy.weave method and takes care of clipping as well. return np.take(lut, data, axis=0, mode='clip') - ### old methods: - #data = np.clip(data, 0, lut.shape[0]-1) - - #try: - #if not USE_WEAVE: - #raise Exception('Weave is disabled; falling back to slower version.') - - ### number of values to copy for each LUT lookup - #if lut.ndim == 1: - #ncol = 1 - #else: - #ncol = sum(lut.shape[1:]) - - ### output array - #newData = np.empty((data.size, ncol), dtype=lut.dtype) - - ### flattened input arrays - #flatData = data.flatten() - #flatLut = lut.reshape((lut.shape[0], ncol)) - - #dataSize = data.size - - ### strides for accessing each item - #newStride = newData.strides[0] / newData.dtype.itemsize - #lutStride = flatLut.strides[0] / flatLut.dtype.itemsize - #dataStride = flatData.strides[0] / flatData.dtype.itemsize - - ### strides for accessing individual values within a single LUT lookup - #newColStride = newData.strides[1] / newData.dtype.itemsize - #lutColStride = flatLut.strides[1] / flatLut.dtype.itemsize - - #code = """ - - #for( int i=0; i (abs(dx[0]) / 1000.)) if not uniform: - import scipy.interpolate as interp + try: + import scipy.interpolate as interp + except: + raise Exception('Fourier transform of irregularly-sampled data requires the package scipy.interpolate.') x2 = np.linspace(x[0], x[-1], len(x)) y = interp.griddata(x, y, x2, method='linear') x = x2 diff --git a/pyqtgraph/graphicsItems/ROI.py b/pyqtgraph/graphicsItems/ROI.py index c0f9008c..79be91f8 100644 --- a/pyqtgraph/graphicsItems/ROI.py +++ b/pyqtgraph/graphicsItems/ROI.py @@ -13,11 +13,8 @@ of how to build an ROI at the bottom of the file. """ from ..Qt import QtCore, QtGui -#if not hasattr(QtCore, 'Signal'): - #QtCore.Signal = QtCore.pyqtSignal import numpy as np -from numpy.linalg import norm -import scipy.ndimage as ndimage +#from numpy.linalg import norm from ..Point import * from ..SRTTransform import SRTTransform from math import cos, sin From b398ccd0ce9cbe5e3033ff645d67c2ae782a3aa8 Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Mon, 10 Mar 2014 13:06:39 -0400 Subject: [PATCH 2/9] corrected import --- examples/ImageView.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/ImageView.py b/examples/ImageView.py index d0bbd31b..44821f42 100644 --- a/examples/ImageView.py +++ b/examples/ImageView.py @@ -14,7 +14,7 @@ displaying and analyzing 2D and 3D data. ImageView provides: import initExample import numpy as np -import scipy +import scipy.ndimage from pyqtgraph.Qt import QtCore, QtGui import pyqtgraph as pg From 18ddff76f0614cef70b1eabbf60f0c1ccb839303 Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Mon, 10 Mar 2014 13:06:54 -0400 Subject: [PATCH 3/9] colormap no longer requires scipy.interpolate --- pyqtgraph/colormap.py | 30 +++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/pyqtgraph/colormap.py b/pyqtgraph/colormap.py index 559e6688..446044e1 100644 --- a/pyqtgraph/colormap.py +++ b/pyqtgraph/colormap.py @@ -63,8 +63,8 @@ class ColorMap(object): ignored. By default, the mode is entirely RGB. =============== ============================================================== """ - self.pos = pos - self.color = color + self.pos = np.array(pos) + self.color = np.array(color) if mode is None: mode = np.ones(len(pos)) self.mode = mode @@ -83,11 +83,6 @@ class ColorMap(object): qcolor Values are returned as an array of QColor objects. =========== =============================================================== """ - try: - import scipy.interpolate - except: - raise Exception("Colormap.map() requires the package scipy.interpolate, but it could not be imported.") - if isinstance(mode, basestring): mode = self.enumMap[mode.lower()] @@ -96,15 +91,24 @@ class ColorMap(object): else: pos, color = self.getStops(mode) - data = np.clip(data, pos.min(), pos.max()) + # don't need this--np.interp takes care of it. + #data = np.clip(data, pos.min(), pos.max()) - if not isinstance(data, np.ndarray): - interp = scipy.interpolate.griddata(pos, color, np.array([data]))[0] + # Interpolate + # TODO: is griddata faster? + # interp = scipy.interpolate.griddata(pos, color, data) + if np.isscalar(data): + interp = np.empty((color.shape[1],), dtype=color.dtype) else: - interp = scipy.interpolate.griddata(pos, color, data) - - if mode == self.QCOLOR: if not isinstance(data, np.ndarray): + data = np.array(data) + interp = np.empty(data.shape + (color.shape[1],), dtype=color.dtype) + for i in range(color.shape[1]): + interp[...,i] = np.interp(data, pos, color[:,i]) + + # Convert to QColor if requested + if mode == self.QCOLOR: + if np.isscalar(data): return QtGui.QColor(*interp) else: return [QtGui.QColor(*x) for x in interp] From 0bc923b7191076be0be9ea18c27fc10e5394a0fe Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Mon, 10 Mar 2014 15:46:55 -0400 Subject: [PATCH 4/9] Added SRTTransform3D test, corrected fromMatrix bug --- pyqtgraph/SRTTransform3D.py | 5 ++-- pyqtgraph/tests/test_srttransform3d.py | 39 ++++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 2 deletions(-) create mode 100644 pyqtgraph/tests/test_srttransform3d.py diff --git a/pyqtgraph/SRTTransform3D.py b/pyqtgraph/SRTTransform3D.py index 7da3c1bd..9b54843b 100644 --- a/pyqtgraph/SRTTransform3D.py +++ b/pyqtgraph/SRTTransform3D.py @@ -122,7 +122,8 @@ class SRTTransform3D(Transform3D): self.setRow(i, m.row(i)) m = self.matrix().reshape(4,4) ## translation is 4th column - self._state['pos'] = m[:3,3] + self._state['pos'] = m[:3,3] + ## scale is vector-length of first three columns scale = (m[:3,:3]**2).sum(axis=0)**0.5 ## see whether there is an inversion @@ -132,7 +133,7 @@ class SRTTransform3D(Transform3D): self._state['scale'] = scale ## rotation axis is the eigenvector with eigenvalue=1 - r = m[:3, :3] / scale[:, np.newaxis] + r = m[:3, :3] / scale[np.newaxis, :] try: evals, evecs = numpy.linalg.eig(r) except: diff --git a/pyqtgraph/tests/test_srttransform3d.py b/pyqtgraph/tests/test_srttransform3d.py new file mode 100644 index 00000000..88aa1581 --- /dev/null +++ b/pyqtgraph/tests/test_srttransform3d.py @@ -0,0 +1,39 @@ +import pyqtgraph as pg +from pyqtgraph.Qt import QtCore, QtGui +import numpy as np +from numpy.testing import assert_array_almost_equal, assert_almost_equal + +testPoints = np.array([ + [0, 0, 0], + [1, 0, 0], + [0, 1, 0], + [0, 0, 1], + [-1, -1, 0], + [0, -1, -1]]) + + +def testMatrix(): + """ + SRTTransform3D => Transform3D => SRTTransform3D + """ + tr = pg.SRTTransform3D() + tr.setRotate(45, (0, 0, 1)) + tr.setScale(0.2, 0.4, 1) + tr.setTranslate(10, 20, 40) + assert tr.getRotation() == (45, QtGui.QVector3D(0, 0, 1)) + assert tr.getScale() == QtGui.QVector3D(0.2, 0.4, 1) + assert tr.getTranslation() == QtGui.QVector3D(10, 20, 40) + + tr2 = pg.Transform3D(tr) + assert np.all(tr.matrix() == tr2.matrix()) + + # This is the most important test: + # The transition from Transform3D to SRTTransform3D is a tricky one. + tr3 = pg.SRTTransform3D(tr2) + assert_array_almost_equal(tr.matrix(), tr3.matrix()) + assert_almost_equal(tr3.getRotation()[0], tr.getRotation()[0]) + assert_array_almost_equal(tr3.getRotation()[1], tr.getRotation()[1]) + assert_array_almost_equal(tr3.getScale(), tr.getScale()) + assert_array_almost_equal(tr3.getTranslation(), tr.getTranslation()) + + From 4263379a90fae23e3e64a40a6b1bf13a7d4289ba Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Mon, 10 Mar 2014 23:45:27 -0400 Subject: [PATCH 5/9] added test for solve3DTransform --- pyqtgraph/functions.py | 17 +++++++++++++---- pyqtgraph/tests/test_functions.py | 21 +++++++++++++++++++++ 2 files changed, 34 insertions(+), 4 deletions(-) create mode 100644 pyqtgraph/tests/test_functions.py diff --git a/pyqtgraph/functions.py b/pyqtgraph/functions.py index 62df1ce3..74d1f8a5 100644 --- a/pyqtgraph/functions.py +++ b/pyqtgraph/functions.py @@ -568,16 +568,25 @@ def transformCoordinates(tr, coords, transpose=False): def solve3DTransform(points1, points2): """ Find a 3D transformation matrix that maps points1 onto points2. - Points must be specified as a list of 4 Vectors. + Points must be specified as either lists of 4 Vectors or + (4, 3) arrays. """ import numpy.linalg - A = np.array([[points1[i].x(), points1[i].y(), points1[i].z(), 1] for i in range(4)]) - B = np.array([[points2[i].x(), points2[i].y(), points2[i].z(), 1] for i in range(4)]) + pts = [] + for inp in (points1, points2): + if isinstance(inp, np.ndarray): + A = np.empty((4,4), dtype=float) + A[:,:3] = inp[:,:3] + A[:,3] = 1.0 + else: + A = np.array([[inp[i].x(), inp[i].y(), inp[i].z(), 1] for i in range(4)]) + pts.append(A) ## solve 3 sets of linear equations to determine transformation matrix elements matrix = np.zeros((4,4)) for i in range(3): - matrix[i] = numpy.linalg.solve(A, B[:,i]) ## solve Ax = B; x is one row of the desired transformation matrix + ## solve Ax = B; x is one row of the desired transformation matrix + matrix[i] = numpy.linalg.solve(pts[0], pts[1][:,i]) return matrix diff --git a/pyqtgraph/tests/test_functions.py b/pyqtgraph/tests/test_functions.py new file mode 100644 index 00000000..da9beca2 --- /dev/null +++ b/pyqtgraph/tests/test_functions.py @@ -0,0 +1,21 @@ +import pyqtgraph as pg +import numpy as np +from numpy.testing import assert_array_almost_equal, assert_almost_equal + +np.random.seed(12345) + +def testSolve3D(): + p1 = np.array([[0,0,0,1], + [1,0,0,1], + [0,1,0,1], + [0,0,1,1]], dtype=float) + + # transform points through random matrix + tr = np.random.normal(size=(4, 4)) + tr[3] = (0,0,0,1) + p2 = np.dot(tr, p1.T).T[:,:3] + + # solve to see if we can recover the transformation matrix. + tr2 = pg.solve3DTransform(p1, p2) + + assert_array_almost_equal(tr[:3], tr2[:3]) From 1eac666d02d7dd6752e33f1090b40b5363d523ad Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Mon, 10 Mar 2014 23:54:35 -0400 Subject: [PATCH 6/9] PlotDataItem._fourierTransform now uses np.interp --- pyqtgraph/graphicsItems/PlotDataItem.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/pyqtgraph/graphicsItems/PlotDataItem.py b/pyqtgraph/graphicsItems/PlotDataItem.py index 3b8d7061..5806475d 100644 --- a/pyqtgraph/graphicsItems/PlotDataItem.py +++ b/pyqtgraph/graphicsItems/PlotDataItem.py @@ -651,16 +651,12 @@ class PlotDataItem(GraphicsObject): def _fourierTransform(self, x, y): ## Perform fourier transform. If x values are not sampled uniformly, - ## then use interpolate.griddata to resample before taking fft. + ## then use np.interp to resample before taking fft. dx = np.diff(x) uniform = not np.any(np.abs(dx-dx[0]) > (abs(dx[0]) / 1000.)) if not uniform: - try: - import scipy.interpolate as interp - except: - raise Exception('Fourier transform of irregularly-sampled data requires the package scipy.interpolate.') x2 = np.linspace(x[0], x[-1], len(x)) - y = interp.griddata(x, y, x2, method='linear') + y = np.interp(x2, x, y) x = x2 f = np.fft.fft(y) / len(y) y = abs(f[1:len(f)/2]) From 72a3902a1865ade1df05923a748b208210dabd10 Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Tue, 11 Mar 2014 02:25:05 -0400 Subject: [PATCH 7/9] Added pure-python implementation of scipy.ndimage.map_coordinates --- pyqtgraph/functions.py | 60 +++++++++++++++++-- pyqtgraph/graphicsItems/ROI.py | 99 ------------------------------- pyqtgraph/tests/test_functions.py | 21 +++++++ 3 files changed, 76 insertions(+), 104 deletions(-) diff --git a/pyqtgraph/functions.py b/pyqtgraph/functions.py index 74d1f8a5..d63a8bbb 100644 --- a/pyqtgraph/functions.py +++ b/pyqtgraph/functions.py @@ -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 ` @@ -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. diff --git a/pyqtgraph/graphicsItems/ROI.py b/pyqtgraph/graphicsItems/ROI.py index 79be91f8..5494f3e3 100644 --- a/pyqtgraph/graphicsItems/ROI.py +++ b/pyqtgraph/graphicsItems/ROI.py @@ -1082,105 +1082,6 @@ class ROI(GraphicsObject): #mapped += translate.reshape((2,1,1)) 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)): """ diff --git a/pyqtgraph/tests/test_functions.py b/pyqtgraph/tests/test_functions.py index da9beca2..86a7ff2c 100644 --- a/pyqtgraph/tests/test_functions.py +++ b/pyqtgraph/tests/test_functions.py @@ -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) + + From ff697ce4929457097baa038381556ae96453a3ad Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Tue, 11 Mar 2014 13:13:33 -0400 Subject: [PATCH 8/9] Expanded capabilities of interpolateArray function to support broadcasting --- pyqtgraph/debug.py | 66 +++++++++--------- pyqtgraph/functions.py | 112 ++++++++++++++++++++++++------ pyqtgraph/tests/test_functions.py | 34 +++++++-- 3 files changed, 151 insertions(+), 61 deletions(-) diff --git a/pyqtgraph/debug.py b/pyqtgraph/debug.py index 6b01c339..f208f3a5 100644 --- a/pyqtgraph/debug.py +++ b/pyqtgraph/debug.py @@ -399,7 +399,9 @@ class Profiler(object): only the initial "pyqtgraph." prefix from the module. """ - _profilers = os.environ.get("PYQTGRAPHPROFILE", "") + _profilers = os.environ.get("PYQTGRAPHPROFILE", None) + _profilers = _profilers.split(",") if _profilers is not None else [] + _depth = 0 _msgs = [] @@ -415,38 +417,36 @@ class Profiler(object): _disabledProfiler = DisabledProfiler() - if _profilers: - _profilers = _profilers.split(",") - def __new__(cls, msg=None, disabled='env', delayed=True): - """Optionally create a new profiler based on caller's qualname. - """ - if disabled is True: - return cls._disabledProfiler - - # determine the qualified name of the caller function - caller_frame = sys._getframe(1) - try: - caller_object_type = type(caller_frame.f_locals["self"]) - except KeyError: # we are in a regular function - qualifier = caller_frame.f_globals["__name__"].split(".", 1)[1] - else: # we are in a method - qualifier = caller_object_type.__name__ - func_qualname = qualifier + "." + caller_frame.f_code.co_name - if func_qualname not in cls._profilers: # don't do anything - return cls._disabledProfiler - # create an actual profiling object - cls._depth += 1 - obj = super(Profiler, cls).__new__(cls) - obj._name = msg or func_qualname - obj._delayed = delayed - obj._markCount = 0 - obj._finished = False - obj._firstTime = obj._lastTime = ptime.time() - obj._newMsg("> Entering " + obj._name) - return obj - else: - def __new__(cls, delayed=True): - return lambda msg=None: None + def __new__(cls, msg=None, disabled='env', delayed=True): + """Optionally create a new profiler based on caller's qualname. + """ + if disabled is True or (disabled=='env' and len(cls._profilers) == 0): + return cls._disabledProfiler + + # determine the qualified name of the caller function + caller_frame = sys._getframe(1) + try: + caller_object_type = type(caller_frame.f_locals["self"]) + except KeyError: # we are in a regular function + qualifier = caller_frame.f_globals["__name__"].split(".", 1)[1] + else: # we are in a method + qualifier = caller_object_type.__name__ + func_qualname = qualifier + "." + caller_frame.f_code.co_name + if disabled=='env' and func_qualname not in cls._profilers: # don't do anything + return cls._disabledProfiler + # create an actual profiling object + cls._depth += 1 + obj = super(Profiler, cls).__new__(cls) + obj._name = msg or func_qualname + obj._delayed = delayed + obj._markCount = 0 + obj._finished = False + obj._firstTime = obj._lastTime = ptime.time() + obj._newMsg("> Entering " + obj._name) + return obj + #else: + #def __new__(cls, delayed=True): + #return lambda msg=None: None def __call__(self, msg=None): """Register or print a new message with timing information. diff --git a/pyqtgraph/functions.py b/pyqtgraph/functions.py index d63a8bbb..57d9a685 100644 --- a/pyqtgraph/functions.py +++ b/pyqtgraph/functions.py @@ -416,6 +416,7 @@ def affineSlice(data, shape, origin, vectors, axes, order=1, returnCoords=False, have_scipy = True except ImportError: have_scipy = False + have_scipy = False # sanity check if len(shape) != len(vectors): @@ -452,14 +453,18 @@ def affineSlice(data, shape, origin, vectors, axes, order=1, returnCoords=False, #print "X values:" #print x ## 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): - ind = (Ellipsis,) + inds - if have_scipy: + if have_scipy: + extraShape = data.shape[len(axes):] + output = np.empty(tuple(shape) + extraShape, dtype=data.dtype) + for inds in np.ndindex(*extraShape): + ind = (Ellipsis,) + inds output[ind] = scipy.ndimage.map_coordinates(data[ind], x, order=order, **kargs) - else: - output[ind] = mapCoordinates(data[ind], x.T) + 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 = list(range(output.ndim)) trb = [] @@ -476,51 +481,114 @@ def affineSlice(data, shape, origin, vectors, axes, order=1, returnCoords=False, else: return output -def mapCoordinates(data, x, default=0.0): +def interpolateArray(data, x, default=0.0): """ - Pure-python alternative to scipy.ndimage.map_coordinates + N-dimensional interpolation similar scipy.ndimage.map_coordinates. - *data* is an array of any shape. - *x* is an array with shape[-1] == data.ndim + 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. Returns array of shape (x.shape[:-1] + data.shape) + + For example, assume we have the following 2D image data:: + + >>> data = np.array([[1, 2, 4 ], + [10, 20, 40 ], + [100, 200, 400]]) + + To compute a single interpolated point from this data:: + + >>> x = np.array([(0.5, 0.5)]) + >>> interpolateArray(data, x) + array([ 8.25]) + + To compute a 1D list of interpolated locations:: + + >>> x = np.array([(0.5, 0.5), + (1.0, 1.0), + (1.0, 2.0), + (1.5, 0.0)]) + >>> interpolateArray(data, x) + array([ 8.25, 20. , 40. , 55. ]) + + To compute a 2D array of interpolated locations:: + + >>> x = np.array([[(0.5, 0.5), (1.0, 2.0)], + [(1.0, 1.0), (1.5, 0.0)]]) + >>> interpolateArray(data, x) + array([[ 8.25, 40. ], + [ 20. , 55. ]]) + + ..and so on. The *x* argument may have any shape as long as + ```x.shape[-1] <= data.ndim```. In the case that + ```x.shape[-1] < data.ndim```, then the remaining axes are simply + broadcasted as usual. For example, we can interpolate one location + from an entire row of the data:: + + >>> x = np.array([[0.5]]) + >>> interpolateArray(data, x) + array([[ 5.5, 11. , 22. ]]) + + This is useful for interpolating from arrays of colors, vertexes, etc. """ + + prof = debug.Profiler() + result = np.empty(x.shape[:-1] + data.shape, dtype=data.dtype) nd = data.ndim + md = x.shape[-1] # First we generate arrays of indexes that are needed to # extract the data surrounding each point - fields = np.mgrid[(slice(0,2),) * nd] + fields = np.mgrid[(slice(0,2),) * 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(nd): - mask = (xmin[...,ax] >= 0) & (xmax[...,ax] < data.shape[ax]) - totalMask &= mask + 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]] - axisMask = mask.astype(np.ubyte).reshape((1,)*(fields.ndim-1) + mask.shape) - axisIndex *= axisMask + #axisMask = mask.astype(np.ubyte).reshape((1,)*(fields.ndim-1) + mask.shape) + axisIndex[axisIndex < 0] = 0 + axisIndex[axisIndex >= data.shape[ax]] = 0 fieldInds.append(axisIndex) + prof() # Get data values surrounding each requested point # fieldData[..., i] contains all 2**nd values needed to interpolate x[i] fieldData = data[tuple(fieldInds)] - + prof() + ## Interpolate - s = np.empty((nd,) + fieldData.shape, dtype=float) + s = np.empty((md,) + fieldData.shape, dtype=float) dx = x - xmin # reshape fields for arithmetic against dx - for ax in range(nd): + for ax in range(md): f1 = fields[ax].reshape(fields[ax].shape + (1,)*(dx.ndim-1)) - s[ax] = f1 * dx[:,ax] + (1-f1) * (1-dx[:,ax]) + 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(nd): + for i in range(md): result = result.sum(axis=0) + prof() + totalMask.shape = totalMask.shape + (1,) * (nd - md) result[~totalMask] = default + prof() return result diff --git a/pyqtgraph/tests/test_functions.py b/pyqtgraph/tests/test_functions.py index 86a7ff2c..af0dde58 100644 --- a/pyqtgraph/tests/test_functions.py +++ b/pyqtgraph/tests/test_functions.py @@ -21,10 +21,10 @@ def testSolve3D(): assert_array_almost_equal(tr[:3], tr2[:3]) -def test_mapCoordinates(): - data = np.array([[ 0., 1., 2.], - [ 2., 3., 5.], - [ 7., 7., 4.]]) +def test_interpolateArray(): + data = np.array([[ 1., 2., 4. ], + [ 10., 20., 40. ], + [ 100., 200., 400.]]) x = np.array([[ 0.3, 0.6], [ 1. , 1. ], @@ -32,11 +32,33 @@ def test_mapCoordinates(): [ 0.5, 2.5], [ 10. , 10. ]]) - result = pg.mapCoordinates(data, x) + result = pg.interpolateArray(data, x) import scipy.ndimage spresult = scipy.ndimage.map_coordinates(data, x.T, order=1) assert_array_almost_equal(result, spresult) - + # test mapping when x.shape[-1] < data.ndim + x = np.array([[ 0.3, 0], + [ 0.3, 1], + [ 0.3, 2]]) + + r1 = pg.interpolateArray(data, x) + r2 = pg.interpolateArray(data, x[0,:1]) + assert_array_almost_equal(r1, r2) + + + # 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]]]) + + r1 = pg.interpolateArray(data, x) + r2 = scipy.ndimage.map_coordinates(data, x.transpose(2,0,1), order=1) + assert_array_almost_equal(r1, r2) + + + + +if __name__ == '__main__': + test_interpolateArray() \ No newline at end of file From 34802c8aeca50225a716ecd22880bd0758297ff0 Mon Sep 17 00:00:00 2001 From: Luke Campagnola Date: Tue, 11 Mar 2014 19:01:34 -0400 Subject: [PATCH 9/9] Added pg.gaussianFilter, removed all dependency on gaussian_filter --- examples/DataSlicing.py | 1 - examples/FlowchartCustomNode.py | 7 ++--- examples/GLImageItem.py | 5 ++-- examples/GLSurfacePlot.py | 5 ++-- examples/HistogramLUT.py | 3 +- examples/ImageView.py | 3 +- examples/VideoSpeedTest.py | 8 +++-- examples/crosshair.py | 5 ++-- examples/isocurve.py | 3 +- pyqtgraph/flowchart/library/Filters.py | 3 +- pyqtgraph/functions.py | 41 +++++++++++++++++++++++++- 11 files changed, 59 insertions(+), 25 deletions(-) diff --git a/examples/DataSlicing.py b/examples/DataSlicing.py index bd201832..d766e7e3 100644 --- a/examples/DataSlicing.py +++ b/examples/DataSlicing.py @@ -11,7 +11,6 @@ a 2D plane and interpolate data along that plane to generate a slice image import initExample import numpy as np -import scipy from pyqtgraph.Qt import QtCore, QtGui import pyqtgraph as pg diff --git a/examples/FlowchartCustomNode.py b/examples/FlowchartCustomNode.py index 25ea5c77..54c56622 100644 --- a/examples/FlowchartCustomNode.py +++ b/examples/FlowchartCustomNode.py @@ -12,7 +12,6 @@ from pyqtgraph.flowchart.library.common import CtrlNode from pyqtgraph.Qt import QtGui, QtCore import pyqtgraph as pg import numpy as np -import scipy.ndimage app = QtGui.QApplication([]) @@ -44,7 +43,7 @@ win.show() ## generate random input data data = np.random.normal(size=(100,100)) -data = 25 * scipy.ndimage.gaussian_filter(data, (5,5)) +data = 25 * pg.gaussianFilter(data, (5,5)) data += np.random.normal(size=(100,100)) data[40:60, 40:60] += 15.0 data[30:50, 30:50] += 15.0 @@ -90,7 +89,7 @@ class ImageViewNode(Node): ## CtrlNode is just a convenience class that automatically creates its ## control widget based on a simple data structure. class UnsharpMaskNode(CtrlNode): - """Return the input data passed through scipy.ndimage.gaussian_filter.""" + """Return the input data passed through pg.gaussianFilter.""" nodeName = "UnsharpMask" uiTemplate = [ ('sigma', 'spin', {'value': 1.0, 'step': 1.0, 'range': [0.0, None]}), @@ -110,7 +109,7 @@ class UnsharpMaskNode(CtrlNode): # CtrlNode has created self.ctrls, which is a dict containing {ctrlName: widget} sigma = self.ctrls['sigma'].value() strength = self.ctrls['strength'].value() - output = dataIn - (strength * scipy.ndimage.gaussian_filter(dataIn, (sigma,sigma))) + output = dataIn - (strength * pg.gaussianFilter(dataIn, (sigma,sigma))) return {'dataOut': output} diff --git a/examples/GLImageItem.py b/examples/GLImageItem.py index dfdaad0c..581474fd 100644 --- a/examples/GLImageItem.py +++ b/examples/GLImageItem.py @@ -12,7 +12,6 @@ from pyqtgraph.Qt import QtCore, QtGui import pyqtgraph.opengl as gl import pyqtgraph as pg import numpy as np -import scipy.ndimage as ndi app = QtGui.QApplication([]) w = gl.GLViewWidget() @@ -22,8 +21,8 @@ w.setWindowTitle('pyqtgraph example: GLImageItem') ## create volume data set to slice three images from shape = (100,100,70) -data = ndi.gaussian_filter(np.random.normal(size=shape), (4,4,4)) -data += ndi.gaussian_filter(np.random.normal(size=shape), (15,15,15))*15 +data = pg.gaussianFilter(np.random.normal(size=shape), (4,4,4)) +data += pg.gaussianFilter(np.random.normal(size=shape), (15,15,15))*15 ## slice out three planes, convert to RGBA for OpenGL texture levels = (-0.08, 0.08) diff --git a/examples/GLSurfacePlot.py b/examples/GLSurfacePlot.py index 963cf4cf..e9896e07 100644 --- a/examples/GLSurfacePlot.py +++ b/examples/GLSurfacePlot.py @@ -10,7 +10,6 @@ import initExample from pyqtgraph.Qt import QtCore, QtGui import pyqtgraph as pg import pyqtgraph.opengl as gl -import scipy.ndimage as ndi import numpy as np ## Create a GL View widget to display data @@ -29,7 +28,7 @@ w.addItem(g) ## Simple surface plot example ## x, y values are not specified, so assumed to be 0:50 -z = ndi.gaussian_filter(np.random.normal(size=(50,50)), (1,1)) +z = pg.gaussianFilter(np.random.normal(size=(50,50)), (1,1)) p1 = gl.GLSurfacePlotItem(z=z, shader='shaded', color=(0.5, 0.5, 1, 1)) p1.scale(16./49., 16./49., 1.0) p1.translate(-18, 2, 0) @@ -46,7 +45,7 @@ w.addItem(p2) ## Manually specified colors -z = ndi.gaussian_filter(np.random.normal(size=(50,50)), (1,1)) +z = pg.gaussianFilter(np.random.normal(size=(50,50)), (1,1)) x = np.linspace(-12, 12, 50) y = np.linspace(-12, 12, 50) colors = np.ones((50,50,4), dtype=float) diff --git a/examples/HistogramLUT.py b/examples/HistogramLUT.py index 5d66cb5d..4d89dd3f 100644 --- a/examples/HistogramLUT.py +++ b/examples/HistogramLUT.py @@ -7,7 +7,6 @@ Use a HistogramLUTWidget to control the contrast / coloration of an image. import initExample import numpy as np -import scipy.ndimage as ndi from pyqtgraph.Qt import QtGui, QtCore import pyqtgraph as pg @@ -34,7 +33,7 @@ l.addWidget(v, 0, 0) w = pg.HistogramLUTWidget() l.addWidget(w, 0, 1) -data = ndi.gaussian_filter(np.random.normal(size=(256, 256)), (20, 20)) +data = pg.gaussianFilter(np.random.normal(size=(256, 256)), (20, 20)) for i in range(32): for j in range(32): data[i*8, j*8] += .1 diff --git a/examples/ImageView.py b/examples/ImageView.py index 44821f42..22168409 100644 --- a/examples/ImageView.py +++ b/examples/ImageView.py @@ -14,7 +14,6 @@ displaying and analyzing 2D and 3D data. ImageView provides: import initExample import numpy as np -import scipy.ndimage from pyqtgraph.Qt import QtCore, QtGui import pyqtgraph as pg @@ -29,7 +28,7 @@ win.show() win.setWindowTitle('pyqtgraph example: ImageView') ## Create random 3D data set with noisy signals -img = scipy.ndimage.gaussian_filter(np.random.normal(size=(200, 200)), (5, 5)) * 20 + 100 +img = pg.gaussianFilter(np.random.normal(size=(200, 200)), (5, 5)) * 20 + 100 img = img[np.newaxis,:,:] decay = np.exp(-np.linspace(0,0.3,100))[:,np.newaxis,np.newaxis] data = np.random.normal(size=(100, 200, 200)) diff --git a/examples/VideoSpeedTest.py b/examples/VideoSpeedTest.py index 7449c2cd..19c49107 100644 --- a/examples/VideoSpeedTest.py +++ b/examples/VideoSpeedTest.py @@ -13,7 +13,6 @@ import initExample ## Add path to library (just for examples; you do not need th from pyqtgraph.Qt import QtGui, QtCore, USE_PYSIDE import numpy as np import pyqtgraph as pg -import scipy.ndimage as ndi import pyqtgraph.ptime as ptime if USE_PYSIDE: @@ -95,10 +94,13 @@ def mkData(): if ui.rgbCheck.isChecked(): data = np.random.normal(size=(frames,width,height,3), loc=loc, scale=scale) - data = ndi.gaussian_filter(data, (0, 6, 6, 0)) + data = pg.gaussianFilter(data, (0, 6, 6, 0)) else: data = np.random.normal(size=(frames,width,height), loc=loc, scale=scale) - data = ndi.gaussian_filter(data, (0, 6, 6)) + print frames, width, height, loc, scale + data = pg.gaussianFilter(data, (0, 6, 6)) + print data[0] + pg.image(data) if dtype[0] != 'float': data = np.clip(data, 0, mx) data = data.astype(dt) diff --git a/examples/crosshair.py b/examples/crosshair.py index 67d3cc5f..076fab49 100644 --- a/examples/crosshair.py +++ b/examples/crosshair.py @@ -7,7 +7,6 @@ the mouse. import initExample ## Add path to library (just for examples; you do not need this) import numpy as np -import scipy.ndimage as ndi import pyqtgraph as pg from pyqtgraph.Qt import QtGui, QtCore from pyqtgraph.Point import Point @@ -33,8 +32,8 @@ p1.setAutoVisible(y=True) #create numpy arrays #make the numbers large to show that the xrange shows data from 10000 to all the way 0 -data1 = 10000 + 15000 * ndi.gaussian_filter(np.random.random(size=10000), 10) + 3000 * np.random.random(size=10000) -data2 = 15000 + 15000 * ndi.gaussian_filter(np.random.random(size=10000), 10) + 3000 * np.random.random(size=10000) +data1 = 10000 + 15000 * pg.gaussianFilter(np.random.random(size=10000), 10) + 3000 * np.random.random(size=10000) +data2 = 15000 + 15000 * pg.gaussianFilter(np.random.random(size=10000), 10) + 3000 * np.random.random(size=10000) p1.plot(data1, pen="r") p1.plot(data2, pen="g") diff --git a/examples/isocurve.py b/examples/isocurve.py index fa451063..b401dfe1 100644 --- a/examples/isocurve.py +++ b/examples/isocurve.py @@ -10,7 +10,6 @@ import initExample ## Add path to library (just for examples; you do not need th from pyqtgraph.Qt import QtGui, QtCore import numpy as np import pyqtgraph as pg -import scipy.ndimage as ndi app = QtGui.QApplication([]) @@ -18,7 +17,7 @@ app = QtGui.QApplication([]) frames = 200 data = np.random.normal(size=(frames,30,30), loc=0, scale=100) data = np.concatenate([data, data], axis=0) -data = ndi.gaussian_filter(data, (10, 10, 10))[frames/2:frames + frames/2] +data = pg.gaussianFilter(data, (10, 10, 10))[frames/2:frames + frames/2] data[:, 15:16, 15:17] += 1 win = pg.GraphicsWindow() diff --git a/pyqtgraph/flowchart/library/Filters.py b/pyqtgraph/flowchart/library/Filters.py index e5bb2453..b72fbca5 100644 --- a/pyqtgraph/flowchart/library/Filters.py +++ b/pyqtgraph/flowchart/library/Filters.py @@ -2,6 +2,7 @@ from ...Qt import QtCore, QtGui from ..Node import Node from . import functions +from ... import functions as pgfn from .common import * import numpy as np @@ -161,7 +162,7 @@ class Gaussian(CtrlNode): import scipy.ndimage except ImportError: raise Exception("GaussianFilter node requires the package scipy.ndimage.") - return scipy.ndimage.gaussian_filter(data, self.ctrls['sigma'].value()) + return pgfn.gaussianFilter(data, self.ctrls['sigma'].value()) class Derivative(CtrlNode): diff --git a/pyqtgraph/functions.py b/pyqtgraph/functions.py index 57d9a685..f76a71c9 100644 --- a/pyqtgraph/functions.py +++ b/pyqtgraph/functions.py @@ -1120,6 +1120,45 @@ def colorToAlpha(data, color): #raise Exception() return np.clip(output, 0, 255).astype(np.ubyte) + +def gaussianFilter(data, sigma): + """ + Drop-in replacement for scipy.ndimage.gaussian_filter. + + (note: results are only approximately equal to the output of + gaussian_filter) + """ + if np.isscalar(sigma): + sigma = (sigma,) * data.ndim + + baseline = data.mean() + filtered = data - baseline + for ax in range(data.ndim): + s = sigma[ax] + if s == 0: + continue + + # generate 1D gaussian kernel + ksize = int(s * 6) + x = np.arange(-ksize, ksize) + kernel = np.exp(-x**2 / (2*s**2)) + kshape = [1,] * data.ndim + kshape[ax] = len(kernel) + kernel = kernel.reshape(kshape) + + # convolve as product of FFTs + shape = data.shape[ax] + ksize + scale = 1.0 / (abs(s) * (2*np.pi)**0.5) + filtered = scale * np.fft.irfft(np.fft.rfft(filtered, shape, axis=ax) * + np.fft.rfft(kernel, shape, axis=ax), + axis=ax) + + # clip off extra data + sl = [slice(None)] * data.ndim + sl[ax] = slice(filtered.shape[ax]-data.shape[ax],None,None) + filtered = filtered[sl] + return filtered + baseline + def downsample(data, n, axis=0, xvals='subsample'): """Downsample by averaging points together across axis. @@ -1556,7 +1595,7 @@ def traceImage(image, values, smooth=0.5): paths = [] for i in range(diff.shape[-1]): d = (labels==i).astype(float) - d = ndi.gaussian_filter(d, (smooth, smooth)) + d = gaussianFilter(d, (smooth, smooth)) lines = isocurve(d, 0.5, connected=True, extendToEdge=True) path = QtGui.QPainterPath() for line in lines: