diff --git a/sgfilter.py b/sgfilter.py index a3c8681..fca5afb 100644 --- a/sgfilter.py +++ b/sgfilter.py @@ -3,10 +3,29 @@ def sgfilter1d(degree, window_size): if not isinstance(degree, int) or degree < 0: raise ValueError("Degree of polynomial must be positive integers") - if not isinstance(window_size, int) or window_size % 2 == 0 or window_size < 0 : + if not isinstance(window_size, int) or window_size % 2 == 0 or window_size < 0: raise ValueError("Window size must be positive odd integer") - if window_size < degree + 2: - raise ValueError("Degree to high for half width") + if window_size < (degree + 1): + raise ValueError("Degree to high for window size") A = np.array([[n ** i for i in xrange(degree + 1)] for n in xrange(-(window_size - 1)/2, (window_size - 1)/2 + 1)]) Ainv = np.linalg.pinv(A) return Ainv[0] + +def sgfilter2d(degree, window_size): + if not isinstance(degree, int) or degree < 0: + raise ValueError("Degree of polynomial must be positive integers") + if not isinstance(window_size, int) or window_size % 2 == 0 or window_size < 0 : + raise ValueError("Window size must be positive odd integer") + if window_size ** 2 < ((degree + 2) * (degree + 1)) / 2.0: + raise ValueError("Degree too high for window size") + exps = [ {"x": k - n, "y": n } for k in xrange(degree + 1) for n in xrange(k + 1)] + n = np.arange(-(window_size - 1) / 2, (window_size - 1)/2 + 1, dtype = np.float64) + dx = np.repeat(n, window_size) + dy = np.tile(n, [window_size, 1]).reshape(window_size ** 2, ) + A = np.empty((window_size ** 2, len(exps))) + for i, exp in enumerate(exps): + A[:,i] = (dx ** exp["x"]) * (dy ** exp["y"]) + Ainv = np.linalg.pinv(A) + return Ainv[0] + +print sgfilter2d(3, 5)