readme.txt

Compilation of convolution and correlation functions for Python

SciPy has functions built-in:
http://www.scipy.org/doc/api_docs/SciPy.signal.signaltools.html#convolve
but this is just a call to SciPy.signal.sigtools._convolve2d or sigtools._correlateND, and:

    "Your large convolutions are usually done using the Fourier Transform (as 
    the direct method implemented by convolveND will be slow for large data 
    -- though it currently could use some optimizations)."

scipy.signal.fftconvolve is probably the simplest, best way:

http://stackoverflow.com/questions/1100100/fft-based-convolution-and-correlation-in-python

Explicit method here:
http://bradmontgomery.blogspot.com/2007/12/computing-correlation-coefficients-in.html
though FFT-based is much faster, and this can be simplified:

>>> I = numpy.asarray(Image.open('test.jpg').convert('L'))
>>> Image.fromarray(numpy.uint8(I)).show()

numarray had a version with FFT possibilities, but this is now part of NumPy/SciPy, or maybe isn't?
http://structure.usc.edu/numarray/node59.html
http://docs.scipy.org/doc/scipy-0.7.x/reference/generated/scipy.stsci.convolve.correlate2d.html
http://www.scipy.org/doc/api_docs/SciPy.stsci.html

This was moved into scipy.stsci.convolve, but doesn't work on my machine.  This is probably what I'm looking for, though.

fromthread.py

# From http://mail.scipy.org/pipermail/scipy-user/2005-January/003967.html

def convolvefft(arr1,arr2):
       s1 = array(arr1.shape())
       s2 = array(arr2.shape())
       fftsize = s1 + s2 - 1
       # finds the closest power of 2 in each dimension (you may comment this out and compare speeds)
       fftsize= pow(2,ceil(log2(fftsize))) 

       ARR1 = fftn(arr1,fftsize)
       ARR2 = fftn(arr2,fftsize)
      
       RES = ifftn(ARR1*ARR2) 
       #RES=RES[validpart]  # I'm not sure how to get the correct part --- first try would be to just truncate to the shape you wanted
       return RES  

haraszti.py

# From http://www.rzuser.uni-heidelberg.de/~ge6/Programing/convolution.html

from numpy.fft import fft2, ifft2

def Convolve(image1, image2, MinPad=True, pad=True):
	""" Not so simple convolution """

	#Just for comfort:
	FFt = fft2
	iFFt = ifft2

	#The size of the images:
	r1,c1 = image1.shape
	r2,c2 = image2.shape

	#MinPad results simpler padding,smaller images:
	if MinPad:
		r = r1+r2
		c = c1+c2
	else:
		#if the Numerical Recipies says so:
		r = 2*max(r1,r2)
		c = 2*max(c1,c2)
	
	#For nice FFT, we need the power of 2:
	if pad:
		pr2 = int(log(r)/log(2.0) + 1.0 )
		pc2 = int(log(c)/log(2.0) + 1.0 )
		rOrig = r
		cOrig = c
		r = 2**pr2
		c = 2**pc2
	#end of if pad
	
	#numpy fft has the padding built in, which can save us some steps
	#here. The thing is the s(hape) parameter:
	fftimage = FFt(image1,s=(r,c)) * FFt(image2,s=(r,c))

	if pad:
		return ((iFFt(fftimage))[:rOrig,:cOrig].real
	else:
		return (iFFt(fftimage)).real

numarray.py

def _correlate2d_fft(data0, kernel0, output=None, mode="nearest", cval=0.0):
    """_correlate2d_fft does 2d correlation of 'data' with 'kernel', storing
    the result in 'output' using the FFT to perform the correlation.

    supported 'mode's include:
        'nearest'   elements beyond boundary come from nearest edge pixel.
        'wrap'      elements beyond boundary come from the opposite array edge.
        'reflect'   elements beyond boundary come from reflection on same array edge.
        'constant'  elements beyond boundary are set to 'cval' 
    """
    shape = data0.shape
    kshape = kernel0.shape
    oversized = (num.array(shape) + num.array(kshape))

    dy = kshape[0] // 2
    dx = kshape[1] // 2

    kernel = num.zeros(oversized, typecode=num.Float64)
    kernel[:kshape[0], :kshape[1]] = kernel0[::-1,::-1]   # convolution <-> correlation
    data = iraf_frame.frame(data0, oversized, mode=mode, cval=cval)

    complex_result = (isinstance(data, _nt.ComplexType) or 
                      isinstance(kernel, _nt.ComplexType))

    Fdata = fft.fft2d(data)
    del data
    
    Fkernel = fft.fft2d(kernel)
    del kernel
    
    num.multiply(Fdata, Fkernel, Fdata)
    del Fkernel

    if complex_result:
        convolved = fft.inverse_fft2d( Fdata, s=oversized)
    else:
        convolved = fft.inverse_real_fft2d( Fdata, s=oversized)
        
    result = convolved[ kshape[0]-1:shape[0]+kshape[0]-1, kshape[1]-1:shape[1]+kshape[1]-1 ]

    if output is not None:
        output._copyFrom( result )
    else:
        return result

                     
def _correlate2d_naive(data, kernel, output=None, mode="nearest", cval=0.0):
    return _correlate.Correlate2d(kernel, data, output, pix_modes[mode], cval)

def correlate2d(data, kernel, output=None, mode="nearest", cval=0.0, fft=0):
    """correlate2d does 2d correlation of 'data' with 'kernel', storing
    the result in 'output'.

    supported 'mode's include:
        'nearest'   elements beyond boundary come from nearest edge pixel.
        'wrap'      elements beyond boundary come from the opposite array edge.
        'reflect'   elements beyond boundary come from reflection on same array edge.
        'constant'  elements beyond boundary are set to 'cval'

    If fft is True,  the correlation is performed using the FFT, else the
    correlation is performed using the naive approach.

    >>> a = num.arange(20*20., shape=(20,20))
    >>> b = num.ones((5,5), typecode='Float64')
    >>> rn = correlate2d(a, b, fft=0)
    >>> rf = correlate2d(a, b, fft=1)
    >>> num.alltrue(num.ravel(rn-rf<1e-10))
    1
    """
    data, kernel = _fix_data_kernel(data, kernel)
    if fft:
        return _correlate2d_fft(data, kernel, output, mode, cval)
    else:
        return _correlate2d_naive(data, kernel, output, mode, cval)

recipes.py

# From http://www.physics.rutgers.edu/~masud/computing/WPark_recipes_in_python.html#Sec7 (1999)

def conv(x, y):
    """
    Convolution of 2 casual signals, x(t<0) = y(t<0) = 0, using discrete
    summation.
        x*y(t) = \int_{u=0}^t x(u) y(t-u) du = y*x(t)
    where the size of x[], y[], x*y[] are P, Q, N=P+Q-1 respectively.

    """
    P, Q, N = len(x), len(y), len(x)+len(y)-1
    z = []
    for k in range(N):
        t, lower, upper = 0, max(0, k-(Q-1)), min(P-1, k)
        for i in range(lower, upper+1):
            t = t + x[i] * y[k-i]
        z.append(t)
    return z

def corr(x, y):
    """
    Correlation of 2 casual signals, x(t<0) = y(t<0) = 0, using discrete
    summation.
        Rxy(t) = \int_{u=0}^{\infty} x(u) y(t+u) du = Ryx(-t)
    where the size of x[], y[], Rxy[] are P, Q, N=P+Q-1 respectively.

    The Rxy[i] data is not shifted, so relationship with the continuous
    Rxy(t) is preserved.  For example, Rxy(0) = Rxy[0], Rxy(t) = Rxy[i],
    and Rxy(-t) = Rxy[-i].  The data are ordered as follows:
        t:	-(P-1),  -(P-2),  ..., -3,  -2,  -1,  0, 1, 2, 3, ..., Q-2, Q-1
        i:	N-(P-1), N-(P-2), ..., N-3, N-2, N-1, 0, 1, 2, 3, ..., Q-2, Q-1

    """
    P, Q, N = len(x), len(y), len(x)+len(y)-1
    z1=[]
    for k in range(Q):
        t, lower, upper = 0, 0, min(P-1, Q-1-k)
        for i in range(lower, upper+1):
            t = t + x[i] * y[i+k]
        z1.append(t)		# 0, 1, 2, ..., Q-1
    z2=[]
    for k in range(1,P):
        t, lower, upper = 0, k, min(P-1, Q-1+k)
        for i in range(lower, upper+1):
            t = t + x[i] * y[i-k]
        z2.append(t)		# N-1, N-2, ..., N-(P-2), N-(P-1)
    z2.reverse()
    return z1 + z2

def fftconv(x, y):
    """
    FFT convolution using relation
        x*y <==> XY
    where x[] and y[] have been zero-padded to length N, such that N >=
    P+Q-1 and N = 2^n.

    """
    N, X, Y, x_y = len(x), fft(x), fft(y), []
    for i in range(N):
        x_y.append(X[i] * Y[i])
    return ifft(x_y)


def fftcorr(x, y):
    """
    FFT correlation using relation
        Rxy <==> X'Y
    where x[] and y[] have been zero-padded to length N, such that N >=
    P+Q-1 and N = 2^n.

    """
    N, X, Y, Rxy = len(x), fft(x), fft(y), []
    for i in range(N):
        Rxy.append(X[i].conjugate() * Y[i])
    return ifft(Rxy)

SciPy.signal.py

# Version from SciPy.signal

def correlate(in1, in2, mode='full'):
    """Cross-correlate two N-dimensional arrays.

  Description:

     Cross-correlate in1 and in2 with the output size determined by mode.

  Inputs:

    in1 -- an N-dimensional array.
    in2 -- an array with the same number of dimensions as in1.
    mode -- a flag indicating the size of the output
            'valid'  (0): The output consists only of those elements that
                            do not rely on the zero-padding.
            'same'   (1): The output is the same size as the largest input
                            centered with respect to the 'full' output.
            'full'   (2): The output is the full discrete linear
                            cross-correlation of the inputs. (Default)

  Outputs:  (out,)

    out -- an N-dimensional array containing a subset of the discrete linear
           cross-correlation of in1 with in2.

    """
    # Code is faster if kernel is smallest array.
    volume = asarray(in1)
    kernel = asarray(in2)
    if rank(volume) == rank(kernel) == 0:
        return volume*kernel
    if (product(kernel.shape,axis=0) > product(volume.shape,axis=0)):
        temp = kernel
        kernel = volume
        volume = temp
        del temp

    val = _valfrommode(mode)

    return sigtools._correlateND(volume, kernel, val)



def correlate2d(in1, in2, mode='full', boundary='fill', fillvalue=0):
    """Cross-correlate two 2-dimensional arrays.

  Description:

     Cross correlate in1 and in2 with output size determined by mode
     and boundary conditions determined by boundary and fillvalue.

  Inputs:

    in1 -- a 2-dimensional array.
    in2 -- a 2-dimensional array.
    mode -- a flag indicating the size of the output
            'valid'  (0): The output consists only of those elements that
                            do not rely on the zero-padding.
            'same'   (1): The output is the same size as the input centered
                            with respect to the 'full' output.
            'full'   (2): The output is the full discrete linear convolution
                            of the inputs. (*Default*)
    boundary -- a flag indicating how to handle boundaries
                'fill' : pad input arrays with fillvalue. (*Default*)
                'wrap' : circular boundary conditions.
                'symm' : symmetrical boundary conditions.
    fillvalue -- value to fill pad input arrays with (*Default* = 0)

  Outputs:  (out,)

    out -- a 2-dimensional array containing a subset of the discrete linear
           cross-correlation of in1 with in2.

    """
    val = _valfrommode(mode)
    bval = _bvalfromboundary(boundary)

    return sigtools._convolve2d(in1, in2, 0,val,bval,fillvalue)