def _activate(self, inputs, out, cval=np.nan):
     """Generate input activities for neurons (convolve and add bias)."""
     out[...] = 0.
     ndimage.correlate(inputs, self.filter, mode='constant', cval=cval, 
         output=out)
     if self.biases.size > 0:
         out[self.active_slice] += self.biases

def get_part_scores(E,log_part_blocks,log_invpart_blocks,
                    frequency_mode='valid',time_mode='same'):
    part_block_shape = log_part_blocks[0].shape
    if time_mode=='valid':
        t_start = part_block_shape[0]/2
        t_end = t_start + E.shape[0] - part_block_shape[0] + 1
    else: # feature_mode == 'same'
        # TODO: check whether I should make the modes available
        #       type-safe
        t_start = 0
        t_end = E.shape[0]
    # check the frequency parameters
    if frequency_mode=='valid':
        f_start = part_block_shape[1]/2
        f_end = f_start + E.shape[1] - part_block_shape[1] + 1
    else: # frequency_mode == 'same'
        f_start = 0
        f_end = E.shape[1]
    e_pos = part_block_shape[2]/2
    return np.array([
            (ndimage.correlate(E,log_part_block)[t_start:t_end,
                                                 f_start:f_end,
                                                 e_pos]
             +ndimage.correlate(1-E,log_invpart_block)[t_start:t_end,
                                                       f_start:f_end,
                                                       e_pos])
            for log_part_block,log_invpart_block in zip(log_part_blocks,
                                                        log_invpart_blocks)])

        def step(self):
            """
            Perform single automaton step.
            """

            gsum = sum(self.grid)

            # Compute burn touched cells
            maximum(1, self.grid, self.burn_map)
            self.burn_map -= 1

            # Correlate cells for next set of fires
            correlate(self.burn_map, self.spread, mode='constant', cval=0,
                output=self.next_burn_map)

            # And cutoff at 1 and multiply by grid to remove
            # barren cells.
            self.next_burn_map *= self.grid
            minimum(1, self.next_burn_map, self.next_burn_map)

            # Finally ignite next set of trees and top at barren
            self.grid += self.next_burn_map
            self.grid += self.burn_map
            minimum(3, self.grid, self.grid)

            if p.sleep:
                __import__('time').sleep(p.sleep)

            # No more fire?
            return gsum < sum(self.grid)

def sobel(img):
    '''
    edges = sobel(img)


    Compute edges using Sobel's algorithm

    `edges` is a binary image of edges computed according to Matlab's
    implementation of Sobel's algorithm.


    Parameters
    ----------
      img : Any 2D-ndarray
    Returns
    -------
      edges : Binary image of edges
    '''
    # This is based on Octave's implementation,
    # but with some reverse engineering to match Matlab exactly
    img = img.astype(np.float)
    img = (img-img.min())/img.ptp()
    vfiltered = ndimage.correlate(img, _vsobel_filter, mode='nearest') # This emulates Matlab's implementation
    hfiltered = ndimage.correlate(img, _hsobel_filter, mode='nearest')
    filtered = vfiltered**2 + hfiltered**2
    thresh = 2*np.sqrt(filtered.mean())
    filtered *= (np.sqrt(filtered) < thresh)

    r,c = filtered.shape
    x = (filtered > np.hstack((np.zeros((r,1)),filtered[:,:-1]))) & (filtered > np.hstack((filtered[:,1:], np.zeros((r,1)))))
    y = (filtered > np.vstack((np.zeros(c),filtered[:-1,:]))) & (filtered > np.vstack((filtered[1:,:], np.zeros(c))))
    return x | y

    def warp(self):
        if self.display:
            print 'Warp %d' % (self.counter,)

        u0, v0 = self.u.copy(), self.v.copy()

        idxx = self.idx + u0
        idyy = self.idy + v0

        idxx = np.clip(idxx, 0, self.N-1)
        idyy = np.clip(idyy, 0, self.M-1)

        I1warped = interp2linear(self.I1, idxx, idyy).astype(theano.config.floatX)
        if self.display:
            #plt.figure()
            #plt.imshow(I1warped)
            print "I1warped", I1warped
            print "I0", self.I0
            print "u0", u0
            print "v0", v0
            pass

        It = I1warped - self.I0
        if self.display:
            print "It", It
            
        Ix = ndimage.correlate(I1warped, self.mask, mode='nearest')
        Iy = ndimage.correlate(I1warped, self.mask.T, mode='nearest')

        # boundary handling
        m = (idxx > self.N - 1) | (idxx < 0) | (idyy > self.M - 1) | (idyy < 0)
        Ix[m] = 0.0
        Iy[m] = 0.0
        It[m] = 0.0

        self.Ix = Ix
        self.Iy = Iy
        self.It = It

        self.train.init(np.zeros_like(self.I0), np.zeros_like(self.I0), Ix, Iy, It)
        cnt_step = 0
        while not self.train.done():
            e = self.train.step()[0]
            if self.display and (cnt_step % 10 == 0 or self.train.done()):
                print e,
            cnt_step += 1

        if self.display:
            print # finish line


        self.u += self.train.tu.get_value()
        self.v += self.train.tv.get_value()
        self.counter += 1

        if self._intermediate_saver:
            self._intermediate_saver.save_members(self)
            self._intermediate_saver.save_locals(locals())

def _mean_std(image, w):
    """Return local mean and standard deviation of each pixel using a
    neighborhood defined by a rectangular window size ``w``.
    The algorithm uses integral images to speedup computation. This is
    used by :func:`threshold_niblack` and :func:`threshold_sauvola`.

    Parameters
    ----------
    image : ndarray
        Input image.
    w : int, or iterable of int
        Window size specified as a single odd integer (3, 5, 7, …),
        or an iterable of length ``image.ndim`` containing only odd
        integers (e.g. ``(1, 5, 5)``).

    Returns
    -------
    m : ndarray of float, same shape as ``image``
        Local mean of the image.
    s : ndarray of float, same shape as ``image``
        Local standard deviation of the image.

    References
    ----------
    .. [1] F. Shafait, D. Keysers, and T. M. Breuel, "Efficient
           implementation of local adaptive thresholding techniques
           using integral images." in Document Recognition and
           Retrieval XV, (San Jose, USA), Jan. 2008.
           :DOI:`10.1117/12.767755`
    """
    if not isinstance(w, Iterable):
        w = (w,) * image.ndim
    _validate_window_size(w)

    pad_width = tuple((k // 2 + 1, k // 2) for k in w)
    padded = np.pad(image.astype('float'), pad_width,
                    mode='reflect')
    padded_sq = padded * padded

    integral = integral_image(padded)
    integral_sq = integral_image(padded_sq)

    kern = np.zeros(tuple(k + 1 for k in w))
    for indices in itertools.product(*([[0, -1]] * image.ndim)):
        kern[indices] = (-1) ** (image.ndim % 2 != np.sum(indices) % 2)

    total_window_size = np.prod(w)
    sum_full = ndi.correlate(integral, kern, mode='constant')
    m = crop(sum_full, pad_width) / total_window_size
    sum_sq_full = ndi.correlate(integral_sq, kern, mode='constant')
    g2 = crop(sum_sq_full, pad_width) / total_window_size
    # Note: we use np.clip because g2 is not guaranteed to be greater than
    # m*m when floating point error is considered
    s = np.sqrt(np.clip(g2 - m * m, 0, None))
    return m, s

def bwmorph_thin(img,n_iter=0):
    # http://www.mathworks.com/help/images/ref/bwmorph.html#f1-500491
    # code adapted from skimage.morphology.skeletonize
    # FIXME support infinite iteration
    LUT=[0,0,0,0,0,1,0,1,0,0,0,0,0,1,3,1,0,0,0,0,2,0,2,0,0,0,
         0,0,3,1,3,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
         2,0,2,0,2,0,0,0,2,0,2,0,0,1,0,1,0,1,0,1,0,0,0,0,0,1,
         0,1,2,0,0,0,2,0,2,0,2,0,0,0,2,0,2,0,0,1,0,1,0,1,0,1,
         0,0,0,0,0,0,0,0,2,0,0,0,2,0,2,0,2,0,0,0,2,0,2,0,0,0,
         0,1,0,1,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
         0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
         0,0,0,0,0,0,0,0,0,0,0,3,0,1,0,1,0,1,0,0,0,0,0,1,0,1,
         2,2,0,0,2,0,0,0,2,2,0,0,2,0,0,0,3,3,0,1,0,1,0,1,0,0,
         0,0,0,0,0,0,2,2,0,0,2,0,0,0,2,2,0,0,2,0,0,0]
    mask = np.array([[ 1, 2, 4],
                     [128, 0, 8],
                     [ 64, 32, 16]], np.uint8)

    skeleton = np.array(img).astype(np.uint8)
    
    done = False
    while(not(done)):
        # assign each pixel a unique value based on its foreground neighbours
        neighbours = ndimage.correlate(skeleton, mask, mode='constant')
        # ignore background
        neighbours *= skeleton
        # use LUT to extract code for deletion stages of thinning algorithm
        codes = np.take(LUT, neighbours)

        done = True
        # pass 1 - remove the 1's
        code_mask = (codes == 1)
        if np.any(code_mask):
            done = False
            skeleton[code_mask] = 0

        neighbours = ndimage.correlate(skeleton, mask, mode='constant')
        # ignore background
        neighbours *= skeleton
        # use LUT to extract code for deletion stages of thinning algorithm
        codes = np.take(LUT, neighbours)

        # pass 2, delete the 2's
        code_mask = (codes > 1)
        if np.any(code_mask):
            done = False
            skeleton[code_mask] = 0

        n_iter -= 1
        if n_iter == 0:
            done = True

    return skeleton.astype(np.bool)

def _mean_std(image, w):
    """Return local mean and standard deviation of each pixel using a
    neighborhood defined by a rectangular window with size w times w.
    The algorithm uses integral images to speedup computation. This is
    used by threshold_niblack and threshold_sauvola.

    Parameters
    ----------
    image : ndarray
        Input image.
    w : int
        Odd window size (e.g. 3, 5, 7, ..., 21, ...).

    Returns
    -------
    m : 2-D array of same size of image with local mean values.
    s : 2-D array of same size of image with local standard
        deviation values.

    References
    ----------
    .. [1] F. Shafait, D. Keysers, and T. M. Breuel, "Efficient
           implementation of local adaptive thresholding techniques
           using integral images." in Document Recognition and
           Retrieval XV, (San Jose, USA), Jan. 2008.
           DOI:10.1117/12.767755
    """
    if w == 1 or w % 2 == 0:
        raise ValueError(
            "Window size w = %s must be odd and greater than 1." % w)

    left_pad = w // 2 + 1
    right_pad = w // 2
    padded = np.pad(image.astype('float'), (left_pad, right_pad),
                    mode='reflect')
    padded_sq = padded * padded

    integral = integral_image(padded)
    integral_sq = integral_image(padded_sq)

    kern = np.zeros((w + 1,) * image.ndim)
    for indices in itertools.product(*([[0, -1]] * image.ndim)):
        kern[indices] = (-1) ** (image.ndim % 2 != np.sum(indices) % 2)

    sum_full = ndi.correlate(integral, kern, mode='constant')
    m = crop(sum_full, (left_pad, right_pad)) / (w ** image.ndim)
    sum_sq_full = ndi.correlate(integral_sq, kern, mode='constant')
    g2 = crop(sum_sq_full, (left_pad, right_pad)) / (w ** image.ndim)
    # Note: we use np.clip because g2 is not guaranteed to be greater than
    # m*m when floating point error is considered
    s = np.sqrt(np.clip(g2 - m * m, 0, None))
    return m, s

 def fbcorr(im, f, mode, stride, plugin=None, plugin_kwargs=None):
     assert mode=='same'
     if f.ndim == 3 and f.shape[0] == 1:
         f = f[0, :, :]
         output = correlate(im, f, mode='constant')[::stride, ::stride]
         return output
     elif f.ndim == 4:
         if f.shape[3] == 1:
             f = f[:, :, :, 0]
             output = [correlate(im, fi, mode='nearest')
                     [::stride, ::stride]
                 for fi in f]
             return np.asarray(output).transpose(1, 2, 0)
         else:
             raise NotImplementedError()

def fonts_template(fn=None,ttf=None):
    ttf = ttf or 'c:/windows/fonts/ariali.ttf'
    fn = fn or 'd:/temp/fonts.arrs'

    m = dict()
    font = ImageFont.truetype(ttf, 18)
    for e in  '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ':
        arr = font.getmask(e, mode='L')
        arr = Image.Image()._new(arr)
        arr = misc.fromimage(arr)
        h,w = arr.shape
        
        print '%s:<%s,%s> %s'%(e, h,w, arr[0,0])
        if w<10:
            tmp = numpy.ndarray( (h,10), dtype=arr.dtype )
            tmp.fill(0)
            i = (10-w)/2
            tmp[:,i:i+w] = arr
            arr = tmp
        #arr = arr[3:18,:]
        arr = arr[2:19,:]
        arr = im_norm(arr)
        rs = ndimage.correlate(arr,arr, mode='constant', cval=0.0)
        limit = numpy.max(rs)
        m[e] = (limit,arr)

    cPickle.dump(m, open(fn,'wb'))

  def scipy_convolution( self ):
    """Performs the convolution without OpenCL, as a comparison"""

    im = self.im
    if self.larger_buffer:
      if self.dim == 1:
        im = self.im[self.offset:-self.offset]
      elif self.dim == 2:
        im = self.im[self.offset:-self.offset,self.offset:-self.offset]
      elif self.dim == 3:
        im = self.im[self.offset:-self.offset,self.offset:-self.offset,self.offset:-self.offset]
    
    tstart = time.time()
    
    if self.sep:
      if self.dim == 1:
        out = ndimage.correlate1d( input=im, weights=self.fil.tolist(), axis=0, mode='wrap', origin=0 )
      elif self.dim == 2:
        out = ndimage.correlate1d( input=im, weights=self.fil.tolist(), axis=0, mode='wrap', origin=0 )
        out = ndimage.correlate1d( input=out, weights=self.fil.tolist(), axis=1, mode='wrap', origin=0 )
      elif self.dim == 3:
        out = ndimage.correlate1d( input=im, weights=self.fil.tolist(), axis=0, mode='wrap', origin=0 )
        out = ndimage.correlate1d( input=out, weights=self.fil.tolist(), axis=1, mode='wrap', origin=0 )
        out = ndimage.correlate1d( input=out, weights=self.fil.tolist(), axis=2, mode='wrap', origin=0 )
    else:
      out = ndimage.correlate( im, self.fil, mode='wrap' )
    print "filtered scipy image", time.time() - tstart, "\n", out
    
    assert numpy.array_equal( out, self.out ), "The PyOpenCL result does not match with Scipy"

def get_part_scores(E,log_part_blocks,log_invpart_blocks):
    part_block_shape = log_part_blocks[0].shape
    t_start = part_block_shape[0]/2
    t_end = t_start + E.shape[0] - part_block_shape[0] + 1
    f_start = part_block_shape[1]/2
    f_end = f_start + E.shape[1] - part_block_shape[1] + 1
    e_pos = part_block_shape[2]/2
    return np.array([
            (ndimage.correlate(E,log_part_block)[t_start:t_end,
                                            f_start:f_end,
                                            e_pos]
             +ndimage.correlate(1-E,log_invpart_block)[t_start:t_end,
                                            f_start:f_end,
                                            e_pos])
            for log_part_block,log_invpart_block in zip(log_part_blocks,
                                                        log_invpart_blocks)])

 def compute(self, dataset_pool):
     constants = dataset_pool.get_dataset('constants')
     footprint = constants["FOOTPRINT"]
     lct = ma.filled(self.get_dataset().get_2d_attribute(self.land_cover_types), 0)
     temp = equal(lct, constants[self.lct_type.upper()])
     values = correlate(temp.astype(int32), footprint, mode="reflect")
     return self.get_dataset().flatten_by_id(ma.filled(values, 0))

	def run(self, im, skin_thresh=[-1,1], n_peaks=3):
		'''
		im : color image
		'''
		im_skin = rgb2lab(im.astype(np.int16))[:,:,2]
		self.im_skin = im_skin
		# im_skin = skimage.exposure.equalize_hist(im_skin)
		# im_skin = skimage.exposure.rescale_intensity(im_skin, out_range=[0,1])
		im_skin *= im_skin > skin_thresh[0]
		im_skin *= im_skin < skin_thresh[1]

		skin_match_c = nd.correlate(-im_skin, self.hand_template)
		self.skin_match = skin_match_c

		# Display Predictions - Color Based matching
		optima = peak_local_max(skin_match_c, min_distance=20, num_peaks=n_peaks, exclude_border=False)
		# Visualize
		if len(optima) > 0:
			optima_values = skin_match_c[optima[:,0], optima[:,1]]
			optima_thresh = np.max(optima_values) / 2
			optima = optima.tolist()

			for i,o in enumerate(optima):
				if optima_values[i] < optima_thresh:
					optima.pop(i)
					break
		self.markers = optima

		return self.markers

    def test_combo(self):
        """Test subsampling followed by convolution.
        """
        # Forward.
        var = Variable((2,3))
        kernel = np.array([[1,2,3]]) # 2x3
        fn = vstack([conv(kernel, subsample(var, (2,1)))])
        fn = CompGraph(fn)
        x = np.arange(6)*1.0
        x = np.reshape(x, (2, 3))
        out = np.zeros(fn.output_size)
        fn.forward(x.flatten(), out)
        y = np.zeros((1,3))

        xsub = x[::2,::1]
        y = ndimage.convolve(xsub, kernel, mode='wrap')

        self.assertItemsAlmostEqual(np.reshape(out, y.shape), y)

        # Adjoint.
        x = np.arange(3)*1.0
        x = np.reshape(x, (1, 3))
        out = np.zeros(var.size)
        fn.adjoint(x.flatten(), out)

        y = ndimage.correlate(x, kernel, mode='wrap')
        y2 = np.zeros((2,3))
        y2[::2,:] = y

        self.assertItemsAlmostEqual(np.reshape(out, y2.shape), y2)
        out = np.zeros(var.size)
        fn.adjoint(x.flatten(), out)
        self.assertItemsAlmostEqual(np.reshape(out, y2.shape), y2)

def main():
    print 'loading data ...'
    ds = Util.loadData(r'E:\yalongli\Projects\RecoverOrientationData\ds_nx32ny32nz6.dat')
    qs = Util.loadData(r'E:\yalongli\Projects\RecoverOrientationData\qs_h12s4t3.dat')
    f = Util.loadData(r'E:\yalongli\Projects\RecoverOrientationData\volumeBin_04.dat')
    
    # small test volume
    f = f[:,:,75:85]
    print f.shape
    
    # Calculate 'J' on these 'q's one by one and store the max_J and max_q in
    # non-empty cells
    print 'initialzing J_max, d_max ...'
    inf = -999999.0
    J_max = np.ones_like(f) * inf
    d_max = np.zeros((f.shape[0], f.shape[1], f.shape[2], 3))
    
    print 'calculating J, d ...'
    for i in range(qs.shape[0]):
        J = ndimage.correlate(f, qs[i], mode='constant', cval=1.0) # background = 1.0!
        update = J > J_max
        J_max[update] = J[update]
        d_max[update] = ds[i]
        if i%100 == 0:
            print i, '...'
    
    # save non-empty cells' result, set empyt cells' direction to zero
    background = (f > 0.9)     
    d_max[background] = 0.0
    
    print 'saving ...'
    ori_p = r'E:\yalongli\Projects\RecoverOrientationData\orientationXYZ7585.dat'
    J_p = r'E:\yalongli\Projects\RecoverOrientationData\JXYZ7585.dat'
    Util.dumpData(d_max, ori_p)
    Util.dumpData(J_max, J_p)

    def test_skeletonize_num_neighbours(self):
        # an empty image
        image = np.zeros((300, 300))

        # foreground object 1
        image[10:-10, 10:100] = 1
        image[-100:-10, 10:-10] = 1
        image[10:-10, -100:-10] = 1

        # foreground object 2
        rs, cs = draw.bresenham(250, 150, 10, 280)
        for i in range(10):
            image[rs + i, cs] = 1
        rs, cs = draw.bresenham(10, 150, 250, 280)
        for i in range(20):
            image[rs + i, cs] = 1

        # foreground object 3
        ir, ic = np.indices(image.shape)
        circle1 = (ic - 135)**2 + (ir - 150)**2 < 30**2
        circle2 = (ic - 135)**2 + (ir - 150)**2 < 20**2
        image[circle1] = 1
        image[circle2] = 0
        result = skeletonize(image)

        # there should never be a 2x2 block of foreground pixels in a skeleton
        mask = np.array([[1,  1],
                         [1,  1]], np.uint8)
        blocks = correlate(result, mask, mode='constant')
        assert not numpy.any(blocks == 4)

    def Ahfunc(self, f):
        '''Conjugate transform - convolve with conj. PSF'''
        fs = reshape(f, (self.height, self.width, self.depth))

        d = ndimage.correlate(fs, self.g)
        
        return ravel(d)

def calculate_vector(data1, data2):
    """
    Return translation vector based on correlation.
    """
    if data1.shape != data2.shape:
        raise ValueError("arrays must have equal shape.")
    data3 = ndimage.correlate(data1, data2, mode="constant")
    return [i - int(s / 2) for s, i in zip(data3.shape, np.unravel_index(data3.argmax(), data3.shape))]

    def _correlate(self,arr, k, alpha=0.6):
        limit,V = self.fonts[k]
        rs = ndimage.correlate(arr,V, mode='constant', cval=0.0)
        h,w = rs.shape

        for ij in ( (i,j) for i in xrange(h) for j in xrange(w)):
            rs[ij] = rs[ij] if rs[ij]>limit*alpha else 0
        return ndimage.gaussian_filter(rs, 2)