def test_sg_coeffs_exact():
    polyorder = 4
    window_length = 9
    halflen = window_length // 2

    x = np.linspace(0, 21, 43)
    delta = x[1] - x[0]

    # The data is a cubic polynomial.  We'll use an order 4
    # SG filter, so the filtered values should equal the input data
    # (except within half window_length of the edges).
    y = 0.5 * x ** 3 - x
    h = savgol_coeffs(window_length, polyorder)
    y0 = convolve1d(y, h)
    assert_allclose(y0[halflen:-halflen], y[halflen:-halflen])

    # Check the same input, but use deriv=1.  dy is the exact result.
    dy = 1.5 * x ** 2 - 1
    h = savgol_coeffs(window_length, polyorder, deriv=1, delta=delta)
    y1 = convolve1d(y, h)
    assert_allclose(y1[halflen:-halflen], dy[halflen:-halflen])

    # Check the same input, but use deriv=2. d2y is the exact result.
    d2y = 3.0 * x
    h = savgol_coeffs(window_length, polyorder, deriv=2, delta=delta)
    y2 = convolve1d(y, h)
    assert_allclose(y2[halflen:-halflen], d2y[halflen:-halflen])

def move_nanvar_filter(arr, window, axis=-1):
    "Moving window variance ignoring NaNs, implemented with a filter."
    arr = np.array(arr, copy=False)
    global convolve1d
    if convolve1d is None:
        try:
            from scipy.ndimage import convolve1d
        except ImportError:
            raise ValueError("'filter' method requires SciPy.")
    if axis is None:
        raise ValueError("An `axis` value of None is not supported.")
    if window < 1:
        raise ValueError("`window` must be at least 1.")
    if window > arr.shape[axis]:
        raise ValueError("`window` is too long.")
    arr = arr.astype(float)
    nrr = np.isnan(arr)
    arr[nrr] = 0
    nrr = nrr.astype(int)
    w = np.ones(window, dtype=int)
    x0 = (1 - window) // 2
    convolve1d(nrr, w, axis=axis, mode="constant", cval=0, origin=x0, output=nrr)
    y = convolve1d(arr, w, axis=axis, mode="constant", cval=np.nan, origin=x0)
    y /= window - nrr
    y *= y
    arr *= arr
    convolve1d(arr, w, axis=axis, mode="constant", cval=np.nan, origin=x0, output=arr)
    arr /= window - nrr
    arr -= y
    arr[nrr == window] = np.nan
    return arr

 def test_03_02_adaptive_threshold_different(self):
     r = np.random.RandomState()
     r.seed(31)
     block = r.uniform(size=(10,10))
     i,j = np.mgrid[0:10:2,0:10:2]
     block[i,j] *= .5
     i,j = np.mgrid[0:50,0:50]
     img = block[i%10, j%10] * .5
     #
     # Make the middle higher in intensity
     #
     img[20:30, 20:30] *= 2
     global_threshold = T.get_global_threshold(T.TM_OTSU, block)
     adaptive_threshold = T.get_adaptive_threshold(
         T.TM_OTSU, img, global_threshold,
         adaptive_window_size = 10)
     #
     # Check that the gradients are positive for i,j<15 and negative
     # for i,j>=15
     #
     gradient = convolve1d(adaptive_threshold, [-1, 0, 1], 0)
     self.assertTrue(np.all(gradient[20:25, 20:30] < 0))
     self.assertTrue(np.all(gradient[25:30, 20:30] > 0))
     gradient = convolve1d(adaptive_threshold, [-1, 0, 1], 1)
     self.assertTrue(np.all(gradient[20:30, 20:25] < 0))
     self.assertTrue(np.all(gradient[20:30, 25:30] > 0))

	def conv_2_sep_Y_circ(self,image,xkernel,thetakernel):#,upsample=10):
		# up_image = np.repeat(image,upsample)
		# x_convolved = ndimage.convolve1d(image,xkernel,mode='mirror',axis=1)
		# both_convolved = ndimage.convolve1d(x_convolved,thetakernel,mode='mirror',axis=0)
		x_convolved = ndimage.convolve1d(image,xkernel,mode='constant',axis=1)
		both_convolved = ndimage.convolve1d(x_convolved,thetakernel,mode='wrap',axis=0)
		return both_convolved

def unsharp_masking(X):
    lp = np.array(X)
    for i, ws in zip([0, 1, 2], [50, 50, 25]):
        h = hamming(ws)
        h /= h.sum()
        convolve1d(lp, h, axis=i, output=lp)
    return X - lp

def second_derivatives(array, smooth=2):
    """
    Compute the second derivatives of all dimensions pairs of the input array
    
    :Inputs:
        array:  any ndarray
        smooth: the second derivative are computed by convolving the input array
                by [1,-2,1]. The smooth parameter set how many times this basic
                filter is convoluted with [1,2,1]/2, which smooth it.
    
    :Output:
        Return a tuple of the second derivative arrays in the order (where dd_ij 
        is the the second derivative d2(array)/didj for a N-dimensional array):
        (dd_00, dd_01, ..., dd_0N, dd_11, dd_12, ..., dd_1N, ..., dd_N-1N, dd_NN)

    :Example:
       for 3d array 'volume_array'
       dv_00, dv_01, dv_02, dv_11, dv_12, dv_22 = second_derivatives(volume_array)
       
    See also:
      numpy.gradient    
    """
    # compute the derivative filter
    dd  = [1,-1]
    for i in xrange(smooth):
        dd = _np.convolve(dd,[1,2,1])/2. 
    
    # compute the second derivatives
    res = ()
    for i in xrange(array.ndim):
        tmp = _nd.convolve1d(array,dd,axis=i)
        for j in xrange(i,array.ndim):
            res += _nd.convolve1d(tmp,dd,axis=j),
    
    return res

def move_var_filter(arr, window, axis=-1):
    "Moving window variance implemented with a filter."
    global convolve1d
    if convolve1d is None:
        try:
            from scipy.ndimage import convolve1d
        except ImportError:
            raise ValueError("'filter' method requires SciPy.")
    if axis == None:
        raise ValueError, "An `axis` value of None is not supported."
    if window < 1:  
        raise ValueError, "`window` must be at least 1."
    if window > arr.shape[axis]:
        raise ValueError, "`window` is too long."  
    arr = arr.astype(float)
    w = np.empty(window)
    w.fill(1.0 / window)
    x0 = (1 - window) // 2
    y = convolve1d(arr, w, axis=axis, mode='constant', cval=np.nan, origin=x0)
    y *= y
    arr *= arr
    convolve1d(arr, w, axis=axis, mode='constant', cval=np.nan, origin=x0,
               output=arr)
    arr -= y 
    return arr

def extract_oriented_patches2D( img, r, coordinates, nb_proj=100 ):

    img = img.astype('float64')

    projection_matrix = np.zeros( (nb_proj,2),
                                  dtype='float64' )
    for i in xrange(nb_proj):
        theta = float(i) * 2.0 * math.pi / nb_proj
        projection_matrix[i,0] = math.sin(theta)
        projection_matrix[i,1] = math.cos(theta)

    print "computing gradients..."
    # grad = np.dstack( ( nd.sobel( img, mode='constant', axis=0 ),
    #                     nd.sobel( img, mode='constant', axis=1 ) ) )

    grad = np.dstack( ( nd.convolve1d( img, [-1,1], mode='constant', axis=0 ),
                        nd.convolve1d( img, [-1,1], mode='constant', axis=1 ) ) )
    
    print "projecting gradients..."
    hist = _patches.project_gradients2D( grad, projection_matrix )
    hist = integral_image2D( hist )

    print hist

    print "extracting patches..."
    Y = coordinates[:,0].copy().astype('int32')
    X = coordinates[:,1].copy().astype('int32')
    return _patches.extract_oriented_patches2D( img,
                                                hist,
                                                projection_matrix,
                                                X,
                                                Y,
                                                r )

def move_nanmax_filter(arr, window, axis=-1):
    "Moving window maximium ignoring NaNs, implemented with a filter."
    global maximum_filter1d, convolve1d
    if maximum_filter1d is None:
        try:
            from scipy.ndimage import maximum_filter1d
        except ImportError:
            raise ValueError("'filter' method requires SciPy.")
    if convolve1d is None:
        try:
            from scipy.ndimage import convolve1d
        except ImportError:
            raise ValueError("'filter' method requires SciPy.")
    if axis == None:
        raise ValueError, "An `axis` value of None is not supported."
    if window < 1:  
        raise ValueError, "`window` must be at least 1."
    if window > arr.shape[axis]:
        raise ValueError, "`window` is too long."
    arr = arr.astype(float)
    nrr = np.isnan(arr)
    arr[nrr] = -np.inf
    x0 = (window - 1) // 2
    maximum_filter1d(arr, window, axis=axis, mode='constant', cval=np.nan,
                     origin=x0, output=arr)
    w = np.ones(window, dtype=int)
    nrr = nrr.astype(int)
    x0 = (1 - window) // 2
    convolve1d(nrr, w, axis=axis, mode='constant', cval=0, origin=x0,
               output=nrr)
    arr[nrr == window] = np.nan
    return arr

def seperable_gaussian_convolution():
    F = plt.imread('../images/cameraman.png')

    # Example of Gaussian convolution.
    s = 3
    G = convolve1d(F, Gauss1(s), axis=-1, mode='nearest')
    G = convolve1d(G, Gauss1(s), axis=0, mode='nearest')
    plt.subplot(1, 2, 1)
    plt.imshow(F, cmap=plt.cm.gray)
    plt.subplot(1, 2, 2)
    plt.imshow(G, cmap=plt.cm.gray)
    plt.show()

    # Gaussian convolution function run 100 times for average execution time.
    n = 10
    s_values = np.array([1, 2, 3, 5, 7, 9, 11, 15, 19])
    times = avg_runnning(F, s_values, n)
    times2d = gauss2d.avg_runnning(F, s_values, n)
    plt.clf()
    plt.title("Average execution time versus the scale\n(run {0} times)"
              .format(n))
    plt.xlabel("s")
    plt.ylabel("time (ms)")
    plt.plot(s_values, times, label="1d convolution")
    plt.plot(s_values, times2d, label="2d convolution")
    plt.legend()
    plt.show()

def move_nansum_filter(arr, window, axis=-1):
    """
    Moving sum (ignoring NaNs) along specified axis using the filter method.
    
    Parameters
    ----------
    arr : array_like
        Input array.
    window : int
        The number of elements in the moving window.
    axis : int, optional
        The axis over which to perform the moving sum. By default the moving
        sum is taken over the last axis (-1).
    
    Returns
    -------
    y : ndarray
        The moving sum (ignoring NaNs) of the input array along the specified
        axis.(A window with all NaNs returns NaN for the window sum.) The
        output has the same shape as the input.

    Notes
    -----
    The calculation of the sums uses scipy.ndimage.convolve1d. 

    Examples
    --------
    >>> from bottlechest.slow.move import move_sum_filter
    >>> arr = np.array([1, 2, np.nan, 4, 5, 6, 7])
    >>> move_nansum_filter(arr, window=2, axis=0)
    array([ NaN,   3.,   2.,   4.,   9.,  11.,  13.])

    """
    arr = np.array(arr, copy=False)
    global convolve1d
    if convolve1d is None:
        try:
            from scipy.ndimage import convolve1d
        except ImportError:
            raise ValueError("'filter' method requires SciPy.")
    if axis == None:
        raise ValueError("An `axis` value of None is not supported.")
    if window < 1:  
        raise ValueError("`window` must be at least 1.")
    if window > arr.shape[axis]:
        raise ValueError("`window` is too long.")
    arr = arr.astype(float)
    nrr = np.isnan(arr)
    arr[nrr] = 0
    nrr = nrr.astype(int)
    w = np.ones(window, dtype=int)
    x0 = (1 - window) // 2
    convolve1d(arr, w, axis=axis, mode='constant', cval=np.nan, origin=x0,
               output=arr)
    convolve1d(nrr, w, axis=axis, mode='constant', cval=0, origin=x0,
               output=nrr)
    arr[nrr == window] = np.nan
    return arr

def main():
    s = 2.0
    m = gauss_1(s)

    img = imread('cameraman.png')
    img2 = convolve1d(img, m, axis=0, mode='nearest')
    img2 = convolve1d(img2, m, axis=1, mode='nearest')
    imshow(img2, cmap=cm.gray)
    show()

def simple_gradient(Im):
    g1=np.array([-0.5, 0, 0.5])
    
    Imx=nd.convolve1d(Im, g1, axis=0)
    Imy=nd.convolve1d(Im, g1, axis=1)
    
    ImMag=(Imx**2 +Imy**2)**0.5
    
    return ImMag, Imx, Imy

def local_standardize(X, kernelsize=(17, 17, 15)):
    local_sq = X ** 2
    local_mean = np.asarray(X)
    for axis, ks in enumerate(kernelsize):
        # w = np.ones(ks) / ks
        w = np.hamming(ks)
        w /= w.sum()
        local_sq = convolve1d(local_sq, w, axis=axis, mode="reflect")
        local_mean = convolve1d(local_mean, w, axis=axis, mode="reflect")
    return (X - local_mean) / np.sqrt(local_sq - local_mean ** 2)

def GodinTypeFilter(data, n, axis=0):
    ''' perform 3 times moving average over the specified array axis.
    suitable for time averaging
    '''
    weights = sp.ones((n),sp.float32)/n
    weights2 = sp.ones((n+1),sp.float32)/(n+1)
    data=nd.convolve1d(data, weights, axis=axis, mode='constant')
    data=nd.convolve1d(data, weights, axis=axis, mode='constant')
    data=nd.convolve1d(data, weights2, axis=axis, mode='constant')
    return data

 def getSwitchTime(self, pixel=(0,0), useKernel='step', method='convolve1d'):
     """
     getSwitchTime(pixel, useKernel='step', method="convolve1d")
     
     Return the position of a step in a sequence
     and the left and the right values of the gray level (as a tuple)
     
     Parameters:
     ---------------
     pixel : tuple
         The (x,y) pixel of the image, as (row, column).
     useKernel : string
         step = [1]*5 +[-1]*5
         zero = [1]*5 +[0] + [-1]*5
         both = step & zero, the one with the highest convolution is chosen
     method : string
         For the moment, only the 1D convolution calculation
         with scipy.ndimage.convolve1d is available
     """
     pxTimeSeq = self.pixelTimeSequence(pixel)
     if method == "convolve1d":
         if useKernel == 'step' or useKernel == 'both':
             convolution_of_stepKernel = nd.convolve1d(pxTimeSeq,self.kernel)
             minStepKernel = convolution_of_stepKernel.min()
             switchStepKernel = convolution_of_stepKernel.argmin() +1
             switch = switchStepKernel
             kernel_to_use = 'step'
         if useKernel == 'zero' or useKernel == 'both':
             convolution_of_zeroKernel = nd.convolve1d(pxTimeSeq,self.kernel0)
             minZeroKernel = convolution_of_zeroKernel.min()
             switchZeroKernel = convolution_of_zeroKernel.argmin() + 1
             switch = switchZeroKernel
             kernel_to_use = 'zero'
         if useKernel == 'both':
             if minStepKernel <= minZeroKernel:
                 switch = switchStepKernel
                 kernel_to_use = 'step'
             else:
                 switch = switchZeroKernel
                 kernel_to_use = 'zero'
                 #leftLevel = np.int(np.mean(pxTimeSeq[0:switch])+0.5)
                 #rightLevel = np.int(np.mean(pxTimeSeq[switch+1:])+0.5)
                 #middle = (leftLevel+rightLevel)/2
                 #rightLevelStep = np.int(np.mean(pxTimeSeq[switchStepKernel+1:])+0.5)
                 #if abs(pxTimeSeq[switch]-middle)>abs(pxTimeSeq[switch]-rightLevelStep):
                     #switch = switchStepKernel                    
                 #switch = (switch-1)*(pxTimeSeq[switch]<middle)+switch*(pxTimeSeq[switch]>=middle)
             #switch = switchStepKernel * (minStepKernel<=minZeroKernel/1.1) + switchZeroKernel * (minStepKernel >minZeroKernel/1.1)
     else:
         raise RuntimeError("Method not yet implemented")            
     levels = self._getLevels(pxTimeSeq, switch, kernel_to_use)
     # Now redefine the switch using the correct image number
     switch = self.imageNumbers[switch]
     return switch, levels

def linearconv(C0, wavelet, sliceind, queue, rorc="row"):
    if rorc == "row":
        for i in range(C0.shape[0]):
            C0[i,:] = ndimage.convolve1d(C0[i,:],wavelet)
                # linearconvker(C0[i,:], wavelet, scale)

    elif rorc == "column":
        for i in range(C0.shape[1]):
            C0[:,i] = ndimage.convolve1d(C0[:,i],wavelet)
                # linearconvker(C0[:,i], wavelet, scale)

    queue.put([sliceind,C0])

def move_sum_filter(arr, window, axis=-1):
    """
    Moving window sum along the specified axis using the filter method.
    
    Parameters
    ----------
    arr : ndarray
        Input array.
    window : int
        The number of elements in the moving window.
    axis : int, optional
        The axis over which to perform the moving sum. By default the moving
        sum is taken over the last axis (-1).
    
    Returns
    -------
    y : ndarray
        The moving sum of the input array along the specified axis. The output
        has the same shape as the input.

    Notes
    -----
    The calculation of the sums uses scipy.ndimage.convolve1d. 

    Examples
    --------
    >>> from bottleneck.slow.move import move_sum_filter
    >>> arr = np.array([1, 2, 3, 4])
    >>> move_sum_filter(arr, window=2, axis=0)
    array([ NaN,   3.,   5.,   7.])

    """
    global convolve1d
    if convolve1d is None:
        try:
            from scipy.ndimage import convolve1d
        except ImportError:
            raise ValueError("'filter' method requires SciPy.")
    if axis == None:
        raise ValueError, "An `axis` value of None is not supported."
    if window < 1:  
        raise ValueError, "`window` must be at least 1."
    if window > arr.shape[axis]:
        raise ValueError, "`window` is too long."  
    arr = arr.astype(float)
    w = np.ones(window, dtype=int)
    x0 = (1 - window) // 2
    convolve1d(arr, w, axis=axis, mode='constant', cval=np.nan, origin=x0,
               output=arr)
    return arr

def run_gaussian_convolutions(F, s_values):
    """
    This function calculates the Gaussian convolution for a given image F
    and given values of s in an array and returns an array with the execution
    times.
    """
    times = np.array([])
    for s in s_values:
        start = datetime.datetime.now()
        G = convolve1d(F, Gauss1(s), axis=-1, mode='nearest')
        G = convolve1d(G, Gauss1(s), axis=0, mode='nearest')
        end = datetime.datetime.now()
        diff = (end - start).total_seconds()
        times = np.append(times, diff)
    return times

def tv_norm(x, beta=2):
    """Computes the total variation norm and its gradient. From jcjohnson/cnn-vis."""
    x_diff = ndimage.convolve1d(x, [-1, 1], axis=2, mode='wrap')
    y_diff = ndimage.convolve1d(x, [-1, 1], axis=1, mode='wrap')
    grad_norm2 = x_diff**2 + y_diff**2 + EPS
    grad_norm_beta = grad_norm2**(beta/2)
    loss = np.sum(grad_norm_beta)
    dgrad_norm2 = (beta/2) * grad_norm2**(beta/2 - 1)
    dx_diff = 2 * x_diff * dgrad_norm2
    dy_diff = 2 * y_diff * dgrad_norm2
    dxy_diff = dx_diff + dy_diff
    dx_diff = roll2(dx_diff, (1, 0))
    dy_diff = roll2(dy_diff, (0, 1))
    grad = dxy_diff - dx_diff - dy_diff
    return loss, grad