def cnnbp(net, y):
    n = len(net.layers);
    
    net.e = net.o - y; # error
    net.L = 0.5 * np.sum(np.square(net.e)) / np.shape(net.e)[1];
    
    net.od = np.multiply(net.e, np.multiply(net.o, 1 - net.o));
    net.fvd = np.mat(net.od) * np.mat(net.ffw); # (50, 192)

    if net.layers[n-1]['type'] == 'c':
        net.fvd = np.multiply(net.fvd, np.multiply(net.fv, (1 - net.fv)));
    
    tm, tr, tc = np.shape(net.layers[n - 1]['a'][0]);
    fvnum = tr * tc;
    
    net.layers[n - 1]['d'] = [];
    for j in range(len(net.layers[n - 1]['a'])):
        net.layers[n - 1]['d'].append(np.array(net.fvd[:, j * fvnum : (j + 1) * fvnum]).reshape(tm, tr, tc));
    
    for l in range(n - 2, 0, -1):
        if net.layers[l]['type'] == 'c':
            net.layers[l]['d'] = [];
            for j in range(len(net.layers[l]['a'])):
                dtmp = [];
                for dl in range(len(net.layers[l]['a'][j])):
                    dtmp.append(np.multiply(np.multiply(net.layers[l]['a'][j][dl], 1 - net.layers[l]['a'][j][dl]),
                                (kronecker(net.layers[l + 1]['d'][j][dl], np.ones((net.layers[l + 1]['scale'],net.layers[l + 1]['scale'])) / np.square(float(net.layers[l + 1]['scale']))))));
                net.layers[l]['d'].append(dtmp);
        elif net.layers[l]['type'] == 's':
            net.layers[l]['d'] = [];
            for i in range(len(net.layers[l]['a'])):
                tm, tr, tc = np.shape(net.layers[l]['a'][0]);
                z = np.zeros((tm, tr, tc));
                for j in range(len(net.layers[l + 1]['a'])):
                    ztmp = [];
                    for dl in range(len(net.layers[l + 1]['d'][j])):
                        ztmp.append(convolve2d(net.layers[l + 1]['d'][j][dl], np.rot90(net.layers[l + 1]['k'][i][j], 2), mode='full'));
                    z += ztmp;
                net.layers[l]['d'].append(z);
    
    for l in range(1, n):
        if net.layers[l]['type'] == 'c':
            dk = [];
            for i in range(len(net.layers[l - 1]['a'])):
                dkj = [];
                for j in range(len(net.layers[l]['a'])):
                    tdk = [];
                    for dl in range(len(net.layers[l - 1]['a'][i])):
                        tdk.append(convolve2d(np.rot90(net.layers[l - 1]['a'][i][dl], 2), net.layers[l]['d'][j][dl], mode='valid'));
                    dkj.append(np.sum(tdk, 0) / np.shape(tdk)[0]);
                dk.append(dkj);
            net.layers[l]['dk'] = dk;
            
            net.layers[l]['db'] = [];
            for j in range(len(net.layers[l]['a'])):
                net.layers[l]['db'].append(np.sum(net.layers[l]['d'][j]) / np.shape(net.layers[l]['d'][j])[0]);
    
    net.dffw = np.mat(net.od).T * np.mat(net.fv) / np.shape(net.od)[0];
    net.dffb = np.mean(net.od, 0).T;
    return net;

def convolucion_RGB(img, kernel):

    channelR = np.zeros((img.shape[0], img.shape[1]), "uint8")
    channelG = np.zeros((img.shape[0], img.shape[1]), "uint8")
    channelB = np.zeros((img.shape[0], img.shape[1]), "uint8")

    for x in range(img.shape[0]):
        for y in range(img.shape[1]):
            channelR[x, y] = img[x, y][0]
            channelG[x, y] = img[x, y][1]
            channelB[x, y] = img[x, y][2]

    matrixR = sg.convolve2d(channelR, kernel, "same")
    matrixG = sg.convolve2d(channelG, kernel, "same")
    matrixB = sg.convolve2d(channelB, kernel, "same")

    img_conv = np.zeros((matrixR.shape[0], matrixR.shape[1], 3), "double")

    for x in range(matrixR.shape[0]):
        for y in range(matrixR.shape[1]):
            img_conv[x, y, 0] = matrixR[x, y]
            img_conv[x, y, 1] = matrixG[x, y]
            img_conv[x, y, 2] = matrixB[x, y]

    return img_conv

def sharpen(img, k=11, lo_pass=True, min_diff=0.05, alpha=1.0):
    """
    Sharpens the input image with unsharp masking. See Wikipedia
    (http://en.wikipedia.org/wiki/Unsharp_masking) for details. Note that
    this function only changes pixel values by a maximum of max_difference
    at each iteration.

    Optionally applies a low-pass Gaussian filter at the end, 
    to reduce the effect of high frequency noise.
    """
    
    # sharpen each color channel separately
    out = np.zeros(img.shape)
    sharp_mask = bf.gauss_mask(k)
    blur_mask = bf.gauss_mask(2 * int(1.0 + (k - 1) / 4.0) + 1) 
    for i in range(0, img.shape[2]):
        blurred = convolve2d(img[:, :, i], sharp_mask, 
                             mode="same", boundary="symm")
        diff = img[:, :, i] - blurred
        scaled = alpha * diff

        if lo_pass:
            
            # if necessary, blur each color channel separately
            diff = convolve2d(diff, blur_mask, 
                              mode="same", boundary="symm")

        diff[np.abs(diff) < min_diff] = 0.0
        out[:, :, i] = img[:, :, i] + scaled


    # truncate to [0, 1]
    out = bf.truncate(out)  
    
    return out

def depthwise_conv2d_python_nhwc(input_np, filter_np, stride, padding):
    """Depthwise convolution operator in nchw layout.

    Parameters
    ----------
    input_np : numpy.ndarray
        4-D with shape [batch, in_height, in_width, in_channel]

    filter_np : numpy.ndarray
        4-D with shape [filter_height, filter_width, in_channel, channel_multiplier]

    stride : list / tuple of 2 ints
        [stride_height, stride_width]

    padding : str
        'VALID' or 'SAME'

    Returns
    -------
    output_np : np.ndarray
        4-D with shape [batch, out_height, out_width, out_channel]
    """
    batch, in_height, in_width, in_channel = input_np.shape
    filter_height, filter_width, _, channel_multiplier = filter_np.shape
    if isinstance(stride, int):
        stride_h = stride_w = stride
    else:
        stride_h, stride_w = stride

    # calculate output shape
    if padding == 'VALID':
        out_channel = in_channel * channel_multiplier
        out_height = (in_height - filter_height) // stride_h + 1
        out_width = (in_width - filter_width) // stride_w + 1
        output_np = np.zeros((batch, out_height, out_width, out_channel))
        for i in range(batch):
            for j in range(out_channel):
                output_np[i, :, :, j] = signal.convolve2d(input_np[i, :, :, j//channel_multiplier], \
                    np.rot90(filter_np[:, :, j//channel_multiplier, j%channel_multiplier], 2), \
                    mode='valid')[0:(in_height - filter_height + 1):stride_h, 0:(in_width - filter_height + 1):stride_w]
    if padding == 'SAME':
        out_channel = in_channel * channel_multiplier
        out_height = np.int(np.ceil(float(in_height) / float(stride_h)))
        out_width = np.int(np.ceil(float(in_width) / float(stride_w)))
        output_np = np.zeros((batch, out_height, out_width, out_channel))
        pad_along_height = np.int(np.max((out_height - 1) * stride_h + filter_height - in_height, 0))
        pad_along_width = np.int(np.max((out_width - 1) * stride_w + filter_width - in_width, 0))
        pad_top_tvm = np.int(np.ceil(float(pad_along_height) / 2))
        pad_left_tvm = np.int(np.ceil(float(pad_along_width) / 2))
        pad_top_scipy = np.int(np.ceil(float(filter_height - 1) / 2))
        pad_left_scipy = np.int(np.ceil(float(filter_width - 1) / 2))
        index_h = pad_top_scipy - pad_top_tvm
        index_w = pad_left_scipy - pad_left_tvm
        for i in range(batch):
            for j in range(out_channel):
                output_np[i, :, :, j] = signal.convolve2d(input_np[i, :, :, j//channel_multiplier], \
                    np.rot90(filter_np[:, :, j//channel_multiplier, j%channel_multiplier], 2), \
                    mode='same')[index_h:in_height:stride_h, index_w:in_width:stride_w]

    return output_np

 def feedforward(self, x=None):
     n = len(self.layers);
     if x is None:
         x = self.x;
     self.layers[0].a = [np.reshape(np.array(x) / 255.0, (self.inputsize[1], self.inputsize[2]))];
     inputmaps = self.inputsize[0];
     
     for l in range(1, n):
         self.layers[l].a = [];
         if self.layers[l].type == 'c':
             for j in range(self.layers[l].outputmaps):
                 rs, cs = np.shape(self.layers[l - 1].a[0]);
                 z = np.zeros((rs - self.layers[l].kernelsize + 1, cs - self.layers[l].kernelsize + 1));
                 for i in range(inputmaps):
                     z += convolve2d(self.layers[l - 1].a[i], self.layers[l].kernels[i][j], mode='valid');
                 self.layers[l].a.append(sigmoid(z + self.layers[l].b[j]));
             inputmaps = self.layers[l].outputmaps;
         elif self.layers[l].type == 's':
             for j in range(inputmaps):
                 z = convolve2d(self.layers[l - 1].a[j], np.ones((self.layers[l].scale, self.layers[l].scale)) / np.square(self.layers[l].scale), mode='valid');
                 self.layers[l].a.append(z[::self.layers[l].scale, ::self.layers[l].scale]);
         elif self.layers[l].type == 'r':
             fv = [];
             for j in range(len(self.layers[l - 1].a)):
                 rs, cs = np.shape(self.layers[l - 1].a[j]);
                 if j == 0:
                     fv = np.reshape(self.layers[l - 1].a[j], (rs * cs));
                 else:
                     fv = np.append(fv, np.reshape(self.layers[l - 1].a[j], (rs * cs)), axis = 1);
             self.layers[l].a = fv;
         elif self.layers[l].type == 'o':
             self.layers[l].a = sigmoid(np.mat(self.layers[l - 1].a) * np.mat(self.layers[l - 1].w).T + np.mat(self.layers[l - 1].b).T);
     
     return self.layers[n - 1].a;

def log_likelihoodG(theta, x, y, data, var, size=21):
    """
    Logarithm of the likelihood function.
    """
    #unpack the parameters
    amplitude, center_x, center_y, radius, focus, width_x, width_y = theta

    #1)Generate a model Airy disc
    airy = models.AiryDisk2D(amplitude, center_x, center_y, radius)
    adata = airy.eval(x, y, amplitude, center_x, center_y, radius).reshape((size, size))

    #2)Apply Focus
    f = models.Gaussian2D(1., center_x, center_y, focus, focus, 0.)
    focusdata = f.eval(x, y, 1., center_x, center_y, focus, focus, 0.).reshape((size, size))
    model = signal.convolve2d(adata, focusdata, mode='same')

    #3)Apply CCD diffusion, approximated with a Gaussian
    CCD = models.Gaussian2D(1., size/2.-0.5, size/2.-0.5, width_x, width_y, 0.)
    CCDdata = CCD.eval(x, y, 1., size/2.-0.5, size/2.-0.5, width_x, width_y, 0.).reshape((size, size))
    model = signal.convolve2d(model, CCDdata, mode='same').flatten()

    #true for Gaussian errors, but not really true here because of mixture of Poisson and Gaussian noise
    lnL = - 0.5 * np.sum((data - model)**2 / var)

    return lnL

def log_likelihood(theta, x, y, data, var, size):
    """
    Logarithm of the likelihood function.
    """
    #unpack the parameters
    peak, center_x, center_y, radius, focus, width_x, width_y = theta

    #1)Generate a model Airy disc
    amplitude = _amplitudeFromPeak(peak, center_x, center_y, radius, x_0=int(size[0]/2.-0.5), y_0=int(size[1]/2.-0.5))
    airy = models.AiryDisk2D(amplitude, center_x, center_y, radius)
    adata = airy.eval(x, y, amplitude, center_x, center_y, radius).reshape(size)

    #2)Apply Focus
    f = models.Gaussian2D(1., center_x, center_y, focus, focus, 0.)
    focusdata = f.eval(x, y, 1., center_x, center_y, focus, focus, 0.).reshape(size)
    model = signal.convolve2d(adata, focusdata, mode='same')

    #3)Apply CCD diffusion, approximated with a Gaussian
    CCDdata = np.array([[0.0, width_y, 0.0],
                        [width_x, (1.-width_y-width_y-width_x-width_x), width_x],
                        [0.0, width_y, 0.0]])
    model = signal.convolve2d(model, CCDdata, mode='same').flatten()

    #true for Gaussian errors
    #lnL = - 0.5 * np.sum((data - model)**2 / var)
    #Gary B. said that this should be from the model not data so recompute var (now contains rn**2)
    var += model.copy()
    lnL = - (np.size(var)*np.sum(np.log(var))) - (0.5 * np.sum((data - model)**2 / var))

    return lnL

def ldp(image):
	x,y = image[:,:,0].shape
	#image = image[:,:,2]	
	image = cv2.cvtColor(image, cv2.COLOR_RGB2GRAY);
	ldp_image = np.zeros((x,y,8),dtype='int')
	mask1 = np.array([[-3,-3,5],[-3,0,5],[-3,-3,5]],dtype=int);
	mask2 = np.array([[-3,5,5],[-3,0,5],[-3,-3,-3]],dtype=int);
	for i in xrange(4):
		m1 = mask1;
		m2 = mask2;
		for j in xrange(i):
			m1 = np.rot90(m1);
			m2 = np.rot90(m2);
		ldp_image[:,:,2*i] = convolve2d(image,m1,mode='same');
		ldp_image[:,:,2*i+1] = convolve2d(image,m2,mode='same');
	ldp_hist = np.zeros(image.shape,dtype='uint8');
	max_image = np.zeros(image.shape,dtype='uint8')
	for i in xrange(x):	
		for j in xrange(y):
			a = ldp_image[i,j,:];
			a = abs(a);
			b = np.sort(a);
			thresh = b[4];
			a = a -thresh;
			a[a>0] = 1;
			a[a<=0] = 0;
			c = np.packbits(a,axis=None);
			ldp_hist[i,j] = int(c);
			max_image[i,j] = b[0];
	return(ldp_image,ldp_hist,max_image);

def simple_sift(patch, grid_number):

    I_x = convolve2d(patch, sobel_vertical_kernel, mode='same', boundary='symm')
    I_y = convolve2d(patch, sobel_horizontal_kernel, mode='same', boundary='symm')

    sift_orientation_histogram_bins = np.linspace(-np.pi, np.pi, 9)

    magnitude_weights = np.sqrt(I_x**2 + I_y**2)

    orientation = np.arctan2(I_y, I_x)

    magnitude_weights_blocks = split_matrix_into_blocks(magnitude_weights, grid_number)
    orientation_blocks = split_matrix_into_blocks(orientation, grid_number)

    descriptor = []

    for magnitude_weight_block, orientation_block in zip(magnitude_weights_blocks, orientation_blocks):

        block_descriptor, _ = np.histogram(orientation_block, bins=sift_orientation_histogram_bins, weights=magnitude_weight_block)

        descriptor.extend(block_descriptor)

    descriptor = np.asarray(descriptor)
    descriptor = descriptor / norm(descriptor)
    descriptor[descriptor > 0.2] = 0.2
    descriptor = descriptor / norm(descriptor)

    return descriptor

def fresh():
    URL = 'http://img2.timeinc.net/people/i/2014/sandbox/news/140630/1994/will-smith-600x450.jpg'
    file = cStringIO.StringIO(urllib2.urlopen(URL).read())
    firstIm = Image.open(file)
    im = firstIm.convert('L')
    #im.show()
    im_array = np.asarray(im)
    im = im.convert('RGB')
    fresh_array = im_array.copy()
    #using 2d gaussian filter
    prince_convolved = gaussconvolve2d(fresh_array, 3)
    prince = Image.fromarray(prince_convolved)
    prince = prince.convert('RGB')
    #prince.show()
    #using box filter
    princebox = signal.convolve2d(fresh_array,boxfilter(19),'same')
    princebox = Image.fromarray(princebox)
    princebox = princebox.convert('RGB')
    #sharpening box filter
    princesharp = signal.convolve2d(fresh_array,sharpen(19),'same')
    princesharp = Image.fromarray(princesharp)
    princesharp = princesharp.convert('RGB')
    #display images side by side comparing gaussian and box filter at the bottom
    dims = (im.size[0]*2, im.size[1]*2)
    compare = Image.new('RGB', dims)
    compare.paste(im, (0,0))
    compare.paste(princesharp, (im.size[0],0))
    compare.paste(prince, (0, im.size[1]))
    compare.paste(princebox, (im.size[0], im.size[1]))
    return compare

    def __init__(self, phasemap, thre):
        assert thre > 0 and thre < 1.0

        super(BrayPSMap, self).__init__(*phasemap.data.shape)
        
        nablaX = np.array([ [-0.5,  0, 0.5], [ -1, 0, 1], [-0.5,  0, 0.5] ], dtype = np.float32)
        nablaY = np.array([ [ 0.5,  1, 0.5], [  0, 0, 0], [-0.5, -1,-0.5] ], dtype = np.float32)
        
        def phaseComplement(phase_diff):
            out = np.copy(phase_diff)
            mask = (phase_diff > np.pi )*1
            out -= mask*(2*np.pi)
            mask = (phase_diff < -np.pi )*1
            out += mask*(2*np.pi)
            return out
        
        diffX = np.zeros_like(phasemap.data)
        diffY = np.zeros_like(phasemap.data)
        diffX[:,:,1:] = phaseComplement( phasemap.data[:,:,1:] - phasemap.data[:,:,0:-1])
        diffY[:,1:,:] = phaseComplement( phasemap.data[:,1:,:] - phasemap.data[:,0:-1,:])
        
        for frame in range(self.data.shape[0]):
            img_dx = diffX[frame, :,:]
            img_dy = diffY[frame, :,:]
            img_x = convolve2d(img_dx, nablaY, boundary='symm', mode='same')
            img_y = convolve2d(img_dy, nablaX, boundary='symm', mode='same')
            self.data[frame, :, :] = ((np.abs(img_y+img_x)/(2*np.pi))>thre)*1.0

        self.roi = np.copy(phasemap.roi)
        self.roi = scipy.ndimage.binary_erosion(self.roi, structure=np.ones((2,2))).astype(phasemap.roi.dtype)
        self.data *= self.roi

        self.vmin = 0.0
        self.vmax = 1.0
        self.cmap = 'gray'

    def __calcHarrisMat__(self):
        self.__findGradients__()
        # M = SUM( 3X3 window )
        # H = det(M) - k*(trace(M))*(trace(M))
        # k = 0.04 <=> 0.06 , we will assume it is 0.05
        # det(M) = (Ix*Ix)*(Iy*Iy) - (Ix*Iy)*(Ix*Iy) ,  trace(M) = (Ix*Ix) + (Iy*Iy)
        gaussianKernel1D = cv2.getGaussianKernel(3, self.sigma)
        window = gaussianKernel1D * gaussianKernel1D.transpose()
        # window = np.array([ [ 1 , 1 , 1 ] ,
        #                     [ 1 , 1 , 1 ] ,
        #                     [ 1 , 1 , 1 ] ])

        gradXSquared = self.gradientX * self.gradientX
        gradYSquared = self.gradientY * self.gradientY
        gradXgradY = self.gradientX * self.gradientY

        # Calculate the summation of the window's value ( IxIx, IyIy, IxIy)

        mIxIx = convolve2d(gradXSquared, window, mode='same')
        mIyIy = convolve2d(gradYSquared, window, mode='same')
        mIxIy = convolve2d(gradXgradY, window, mode='same')

        # Calculate the |M|
        detOfMatrix = (mIxIx * mIyIy) - (mIxIy * mIxIy)
        # Calculate the trace()
        traceOfMatrix = mIxIx + mIyIy

        self.responseMat = detOfMatrix - 0.05 * traceOfMatrix * traceOfMatrix

def deformed_patch_based(lim_patch_list, rim_patch_list, image_patch, patch_list):
    deformed_list = []
    #print len(image_patch)
    #print len(patch_list)
    #print np.mean((np.asarray(image_patch)-np.asarray(patch_list))**2)**0.5

    fx = np.asarray([[1.0, -1.0]], dtype='float')
    fy = np.asarray([[1.0], [-1.0]], dtype='float')

    for i in range(len(lim_patch_list)):
        grad_x = convolve2d(patch_list[i], fx, mode='same')*0.1
        grad_y = convolve2d(patch_list[i], fy, mode='same')*0.1

        grad_x[:, 0] = grad_x[:, 0]*0+1
        grad_y[0, :] = grad_y[0, :]*0+1
        # ====================== 这里手动构造梯度 ===================
        Def = Deformed.deformed_patch()
        cp, c = Def.deform(patch_list[i], image_patch[i],grad_x,grad_y)
        deformed_list.append(cp)
        #print np.mean((patch_list[i]-image_patch[i])**2)**0.5
        #print np.mean((cp-image_patch[i])**2)**0.5

        print "变形完成：%0.2f %%" % (i*1.0/len(lim_patch_list)*100)
    print np.mean((np.asarray(deformed_list)-np.asarray(patch_list))**2)**0.5
    return deformed_list

	def prewitt(self, gray_image):
		# Construct the mask
		kx = np.matrix('\
			-1 -1 -1;\
			 0  0  0;\
			 1  1  1 \
			') / 6.0

		ky = np.matrix('\
			-1  0  1;\
			-1  0  1;\
			-1  0  1 \
			') / 6.0

		# Convolute
		gx = convolve2d(gray_image, kx, 'same')
		gy = convolve2d(gray_image, ky, 'same')

		height,width = gray_image.shape[:2]
		result = gray_image.copy()
		max_p = 0
		for y in range(0, height):
			for x in range(0, width):
				result[y,x] = math.sqrt(gx[y,x]**2 + gy[y,x]**2)
				max_p = max(max_p, result[y,x])

		return np.uint8(result > 0.08995 * max_p) * 255

def gradient_magnitude(image):
    """
    Returns the L1 gradient magnitude of the given image.
    The result is an n x m numpy array of floating point values,
    where n is the number of rows in the image and m is the number of columns.
    """
    # First, convert the input image to a 32-bit floating point grayscale.
    # Be sure to scale it such that the intensity varies between 0 and 1.
    Is = np.mean(image,axis=2,dtype=np.float64)
    Is = normalize2D(Is)
    
    # Next, compute the graient in the x direction using the sobel operator with a kernel size of 5
    sobelX = [[2.0,1.0,0.0,-1.0,-2.0],
              [3.0,2.0,0.0,-2.0,-3.0],
              [4.0,3.0,0,-3.0,-4.0],
              [3.0,2.0,0.0,-2.0,-3.0],
              [2.0,1.0,0.0,-1.0,-2.0]]
    #gradientX = applyConvolution(Is,sobelX) #apparently there's a np.convolute2d - and it is so much faster!
    gradientX = signal.convolve2d(Is,sobelX,mode='same')
    # Compute the graient in the y direction using the sobel operator with a kernel size of 5
    gradientY = signal.convolve2d(Is,np.transpose(sobelX),mode='same')
    #gradientY = applyConvolution(Is,np.transpose(sobelX))    

    # Finally, compute the l1 norm of the x and y gradients at each pixel value.
    # The l1 norm is the sum of their absolute values.
    # convert the final image from a double-precision floating point to single.
    gradientX = np.absolute(gradientX)
    gradientY = np.absolute(gradientY)
    energy = gradientX + gradientY
    
    print 'found energy ',energy.shape

    # and return the result
    return energy

def twoDsmooth(mat,ker):
    try:
        len(ker)
        kmat = ker
        
    except:
        kmat = np.ones((ker,ker))
        kmat = kmat / pow(ker, 2)

    [kr,kc] = list(kmat.shape)
    if (kr%2 == 0):
        conmat = np.ones((2,1))
        kmat = signal.convolve2d(kmat,conmat,'symm','same')
        kr = kr + 1

    if (kc%2 == 0):
        conmat = np.ones((2,1))
        kmat = signal.convolve2d(kmat,conmat,'symm','same')
        kc = kc + 1

    [mr,mc] = list(mat.shape)
    fkr = math.floor(kr/2)
    fkc = math.floor(kc/2)
    rota = np.rot90(kmat,2)
    mat=signal.convolve2d(mat,rota,'same','symm')
    return mat

def sobelEdgeFeaturesDigit(datum):
  from scipy import signal
  import numpy

  features = util.Counter()
  opx = numpy.array([[1, 0, -1],
                     [2, 0, -2],
                     [1, 0, -1]])
  opy = numpy.array([[1, 2, 1],
                     [0, 0, 0],
                     [-1, -2, -1]])
  pixels = numpy.empty([DIGIT_DATUM_WIDTH, DIGIT_DATUM_HEIGHT])

  for x in range(DIGIT_DATUM_WIDTH):
    for y in range(DIGIT_DATUM_HEIGHT):
      pixels[x, y] = datum.getPixel(x, y)

  gx = signal.convolve2d(pixels, opx)
  gy = signal.convolve2d(pixels, opy)
  gx = numpy.absolute(gx)
  gy = numpy.absolute(gy)

  for x in range(DIGIT_DATUM_WIDTH):
    for y in range(DIGIT_DATUM_HEIGHT):
      mag = gx[x, y] + gy[x, y]
      features[(x, y, "edge")] = 1 if mag >= 3 else 0

  return features

 def _get_riesz_triple(self, img):
     """:param img: 2D input image of shape (h, w)
     :return: 3D image of shape (h, w, 3), where the 3 channels correspond to
     amplitude, phase * cos(orient), phase * sin(orient)
     """
     assert img.ndim == 2
     img = img.astype(floatX)
     r1 = signal.convolve2d(img, self.riesz_kernel, mode='valid')
     r2 = signal.convolve2d(img, self.riesz_kernel.T, mode='valid')
     kh, kw = self.riesz_kernel.shape
     assert kh % 2 == 1 and kw % 2 == 1 and kh == kw
     img = img[kh/2:-(kh/2), kw/2:-(kw/2)]
     amp = np.sqrt(np.square(img) + np.square(r1) + np.square(r2))
     phase = np.arccos(img / (amp + self.EPS))
     t = np.sqrt(np.square(r1) + np.square(r2)) + self.EPS
     v0 = phase * r1 / t
     v1 = phase * r2 / t
     if self.spatial_blur:
         amp_blur = cv2.GaussianBlur(
             amp, self.spatial_ksize, self.spatial_blur)
         def blur(v):
             a = cv2.GaussianBlur(
                 amp * v, self.spatial_ksize, self.spatial_blur)
             return a / amp_blur
         v0 = blur(v0)
         v1 = blur(v1)
     rst = np.concatenate(map(lambda a: np.expand_dims(a, axis=2),
                              (amp, v0, v1)), axis=2)
     assert np.all(np.isfinite(rst))
     return rst

def log_likelihoodC(theta, x, y, data, var, size=21):
    """
    Logarithm of the likelihood function.
    """
    #unpack the parameters
    amplitude, center_x, center_y, radius, focus, width_x, width_y, width_d = theta

    #1)Generate a model Airy disc
    airy = models.AiryDisk2D(amplitude, center_x, center_y, radius)
    adata = airy.eval(x, y, amplitude, center_x, center_y, radius).reshape((size, size))

    #2)Apply Focus
    f = models.Gaussian2D(1., center_x, center_y, focus, focus, 0.)
    focusdata = f.eval(x, y, 1., center_x, center_y, focus, focus, 0.).reshape((size, size))
    focusmodel = signal.convolve2d(adata, focusdata, mode='same')

    #3)Apply CCD diffusion kernel
    kernel = np.array([[width_d, width_y, width_d],
                       [width_x, 1., width_x],
                       [width_d, width_y, width_d]])
    kernel /= kernel.sum()
    #model = ndimage.convolve(focusmodel, kernel)
    model = signal.convolve2d(focusmodel, kernel, mode='same').flatten()

    #true for Gaussian errors, but not really true here because of mixture of Poisson and Gaussian noise
    lnL = - 0.5 * np.sum((data - model)**2 / var)

    return lnL

    def convolve(self, npImage, filters):
        """
        Return the 2D convolution of each colorband by its respective filter.

        Parameters
        ----------
        npImage : 3D numpy array where the dimensions represent respectively
        the height, the width and the colorbands. There must be 3 colorbands.
            The image to filter
        filters : a triplet of filters of the same size. The filters are
        2D array like structure (usually numpy array).
            The respective filters. They are applied in the same order to the
            colorbands.

        Return
        ------
        filtered :3D numpy array where the dimensions represent respectively
        the height, the width and the colorbands
            The result of the convolution of each colorband by its respective
            filter
        """
        redFilter, greenFilter, blueFilter = filters
        red, green, blue = npImage[:,:,0], npImage[:,:,1], npImage[:,:,2]

        newRed = sg.convolve2d(red, redFilter, "same")
        newGreen = sg.convolve2d(green, greenFilter, "same")
        newBlue = sg.convolve2d(blue, blueFilter, "same")

        return np.dstack((newRed, newGreen, newBlue))