def width_filter(data, angles=None, FILT_SIZE=5):
    if angles == None:
        #estimate angle from intensity data
        return ndimage.generic_filter(data.astype('f'), width, FILT_SIZE, extra_arguments=genCoords(FILT_SIZE))
    else:
        d = np.concatenate([data[:,:,None], angles[:,:,None]], 2)
        return ndimage.generic_filter(d.astype('f'), width_o, [FILT_SIZE, FILT_SIZE, 2], extra_arguments=genCoords(FILT_SIZE))[:,:,0].squeeze()

    def __init__(self, label_image=None, connectivity=1, data=None, **attr):

        super(RAG, self).__init__(data, **attr)
        if self.number_of_nodes() == 0:
            self.max_id = 0
        else:
            self.max_id = max(self.nodes_iter())

        if label_image is not None:
            fp = ndi.generate_binary_structure(label_image.ndim, connectivity)
            # In the next ``ndi.generic_filter`` function, the kwarg
            # ``output`` is used to provide a strided array with a single
            # 64-bit floating point number, to which the function repeatedly
            # writes. This is done because even if we don't care about the
            # output, without this, a float array of the same shape as the
            # input image will be created and that could be expensive in
            # memory consumption.
            ndi.generic_filter(
                label_image,
                function=_add_edge_filter,
                footprint=fp,
                mode='nearest',
                output=as_strided(np.empty((1,), dtype=np.float_),
                                  shape=label_image.shape,
                                  strides=((0,) * label_image.ndim)),
                extra_arguments=(self,))

def prepareInputs(temp_folder, lastoolsPath):

	# lastoolsPath = "C:/lastools/bin/"

	# Run lasground
	# Note: Nasa laszip dapat
	os.chdir(lastoolsPath)

	print "Running LASground..."

	subprocess.call(["lasground_new", "-i", temp_folder + "/pointcloud.laz","-metro", "-compute_height","-odir", temp_folder + "/", "-o","ground.laz"], stdout=subprocess.PIPE)

	print "Running LASClassify..."

	# Prepare file_list.txt

	# subprocess.call(["lasclassify", "-i", "C:/bertud_temp/ground.laz","-odir", "C:/bertud_temp/", "-o","classified.laz"], stdout=subprocess.PIPE)
	
	# Added fine tuning parameter -planar
	subprocess.call(["lasclassify", "-i", temp_folder + "/ground.laz","-planar","0.15","-odir", temp_folder + "/", "-o","classified.laz"], stdout=subprocess.PIPE)

	print "Running LASGrid for classification raster..."

	# Prepare file_list.txt
	# subprocess.call(["lasgrid", "-i", "C:/bertud_temp/classified.laz","-step","0.5","-classification","-odir", "C:/bertud_temp/", "-o","classified.tif"], stdout=subprocess.PIPE)
	
	# Added fine tuning parameter -subsample 8
	subprocess.call(["lasgrid", "-i", temp_folder + "/classified.laz","-step","0.5","-classification","-subsample","8","-odir", temp_folder + "/", "-o","classified.tif"], stdout=subprocess.PIPE)

	print "Running LASGrid for number of returns raster..."
	subprocess.call(["lasgrid", "-i", temp_folder + "/classified.laz","-step","0.5","-number_returns","-lowest","-subsample","8","-odir", temp_folder + "/", "-o","numret.tif"], stdout=subprocess.PIPE)

	print "Running blast2DEM..."
	
	subprocess.call(["blast2dem", "-i", temp_folder + "/classified.laz", "-first_only","-step","0.5","-elevation","-odir", temp_folder + "/", "-o","dsm.tif"], stdout=subprocess.PIPE)
	subprocess.call(["blast2dem", "-i", temp_folder + "/classified.laz", "-keep_classification","2","-keep_classification","8","-step","0.5","-elevation","-odir", temp_folder + "/", "-o","dtm.tif"], stdout=subprocess.PIPE)

	dsm = io.imread(temp_folder + "/dsm.tif")
	dtm = io.imread(temp_folder + "/dtm.tif")

	# nDSM

	dtm[dtm<0] = 9999
	dsm[dsm<0] = 0

	ndsm = dsm-dtm

	# Revised nDSM generation
	# ndsm[ndsm<2] = 0
	ndsm[ndsm<0] = 0

	io.imsave(temp_folder + "/ndsm.tif",ndsm)
	
	# Slope 

	slope = ndimage.generic_filter(ndsm,slopeFilter,size=3)
	io.imsave(temp_folder + "/slope.tif",slope)

	slopeslope = ndimage.generic_filter(slope,slopeFilter,size=3)
	io.imsave(temp_folder + "/slopeslope.tif",slopeslope)

def performSoftMatting(im=None, transmission=None, *args, **kwargs):
    global neighbors

    width, height, depth = im.shape
    windowRadius = 1
    numWindowPixels = 9
    epsilon = 10 ** - 8
    _lambda = 10 ** - 4

    totalPixels = numWindowPixels ** 2
    
    windowIndicies = np.reshape(xrange(1, width * height + 1), (width, height), order='F')
    totalElements = totalPixels * (width - 2) * (height - 2)
    xIndicies = np.ones((1, totalElements))
    yIndicies = np.ones((1, totalElements))
    laplacian = np.zeros((1, totalElements))
    count = 0

    neighbors = np.empty((width * height, numWindowPixels))

    footprint = np.array([[1,1,1],
                          [1,1,1],
                          [1,1,1]])

    ndimage.generic_filter(windowIndicies, getWindow, footprint=footprint)

    U = epsilon / numWindowPixels * identity(3)

    for i in xrange(0 + windowRadius, height - windowRadius):
        for j in xrange(0 + windowRadius, width - windowRadius):
            print 'i', i
            print 'j', j
            window = im[j - windowRadius: j + windowRadius + 1, i - windowRadius : i + windowRadius + 1, :]

            reshapedWindow = np.reshape(window, (numWindowPixels, 3), order='F')

            diffFromMean = reshapedWindow.T - np.tile(np.mean(reshapedWindow, axis=0).T, (numWindowPixels, 1)).T

            window_covariance = np.dot(diffFromMean, diffFromMean.T) / numWindowPixels

            entry = identity(numWindowPixels) - (1 + np.dot(np.dot(diffFromMean.T, np.linalg.inv(window_covariance + U)), diffFromMean)) / float(numWindowPixels)

            temp = count * totalPixels
            temp2 = count * totalPixels + totalPixels

            iterationNeighbors = np.reshape(np.reshape(neighbors[height * j + i], (3, 3)), (1, numWindowPixels), order='F')

            x = np.tile(iterationNeighbors, (numWindowPixels, 1))
            y = (x.T).flatten(1)

            xIndicies[0][temp : temp2] = x.flatten(1)
            yIndicies[0][temp : temp2] = y
            laplacian[0][temp : temp2] = entry.flatten(1)
            count += 1

    L = csc_matrix((laplacian.flatten(), (xIndicies.flatten(), yIndicies.flatten())))
    tBar = np.append(np.reshape(transmission.T, (width * height, 1)), [0])

    T = spla.spsolve(L + _lambda * identity(L.shape[0]), tBar * _lambda) 
    return np.reshape(np.delete(T, len(T) - 1), transmission.shape, order='F')

    def _filter_meshes(self):
        """
        Apply a 2D median filter to the low-resolution 2D mesh,
        including only pixels inside the image at the borders.
        """

        from scipy.ndimage import generic_filter
        try:
            nanmedian_func = np.nanmedian    # numpy >= 1.9
        except AttributeError:    # pragma: no cover
            from scipy.stats import nanmedian
            nanmedian_func = nanmedian

        if self.filter_threshold is None:
            # filter the entire arrays
            self.background_mesh = generic_filter(
                self.background_mesh, nanmedian_func, size=self.filter_size,
                mode='constant', cval=np.nan)
            self.background_rms_mesh = generic_filter(
                self.background_rms_mesh, nanmedian_func,
                size=self.filter_size, mode='constant', cval=np.nan)
        else:
            # selectively filter
            indices = np.nonzero(self.background_mesh > self.filter_threshold)
            self.background_mesh = self._selective_filter(
                self.background_mesh, indices)
            self.background_rms_mesh = self._selective_filter(
                self.background_rms_mesh, indices)

        return

def checkTile(tile,title=''):
    nRow,nCol= np.shape(tile)
    freqList = np.sort(np.reshape(tile,(-1,)))

    nearestDist = neighborDistClass(tile)
    minDists = ndimage.generic_filter(tile, nearestDist.minFilter, footprint=footprint,mode='constant',cval=np.inf)

    nearestDist = neighborDistClass(tile)
    minDists2 = ndimage.generic_filter(tile, nearestDist.minFilter, footprint=secondNeighborFootprint,mode='constant',cval=np.inf)

    nearestDist = neighborDistClass(tile)
    minDistsWrap = ndimage.generic_filter(tile, nearestDist.minFilter, footprint=sideFootprint,mode='wrap',cval=np.inf)

    nearestDist = neighborDistClass(tile)
    maxDists = ndimage.generic_filter(tile, nearestDist.maxFilter, footprint=footprint,mode='reflect')

    plotArray(title=title,image=tile,normNSigma=2.,origin='upper')
    plotArray(title='{} min dists wrap'.format(title),image=minDistsWrap,normNSigma=2.,origin='upper')
    plotArray(title='{} min dists'.format(title),image=minDists,normNSigma=2.,origin='upper')
    plotArray(title='{} min dists 2nd nearest'.format(title),image=minDists2,normNSigma=2.,origin='upper')
    plotArray(title='{} max dists'.format(title),image=maxDists,normNSigma=2.,origin='upper')

#    def f(thing):
#        thing.axes.hist(minDists.ravel(),bins=100)
#        thing.axes.set_title('{} min dists'.format(title))
#        
#    pop = PopUp(plotFunc=f)

    def f(thing):
        thing.axes.hist(minDists2.ravel(),bins=100)
        thing.axes.set_title('{} second neighbor min dists'.format(title))
        
    pop = PopUp(plotFunc=f)

def test_ticket_701():
    # Test generic filter sizes
    arr = np.arange(4).reshape((2,2))
    func = lambda x: np.min(x)
    res = sndi.generic_filter(arr, func, size=(1,1))
    # The following raises an error unless ticket 701 is fixed
    res2 = sndi.generic_filter(arr, func, size=1)
    assert_equal(res, res2)

def get_windowed_vals(x, k=3):
    nrows, ncols = x.shape
    targets = np.concatenate([np.zeros(k - 1, dtype=int), x["target"], np.zeros(k - 1, dtype=int)])
    preds = np.concatenate([np.zeros(k - 1, dtype=int), x["pred"], np.zeros(k - 1, dtype=int)])
    wtargets = generic_filter(targets, np.max, size=(k,), mode="constant")
    wpreds = generic_filter(preds, np.max, size=(k,), mode="constant")
    starts = range(-k, nrows + 1)
    ends = range(nrows + k + 1)
    wvals = pd.DataFrame({"start": starts, "target": wtargets, "pred": wpreds}, index=[ends])
    wvals["conv"] = x["conv"].iloc[0]

    return wvals

def currents_function(ax, data_file, bmap, key_ax, time_index, downsample_ratio):
    def compute_average(array):
        avg = numpy.average(array)
        return numpy.nan if avg > 10**3 else avg

    print "Currents Downsample Ratio:", downsample_ratio

    currents_u = data_file.variables['u'][time_index][39]
    currents_v = data_file.variables['v'][time_index][39]
    rho_mask = get_rho_mask(data_file)

    # average nearby points to align grid, and add the edge column/row so it's the right size.
    #-------------------------------------------------------------------------
    right_column = currents_u[:, -1:]
    currents_u_adjusted = ndimage.generic_filter(scipy.hstack((currents_u, right_column)),
                                                 compute_average, footprint=[[1], [1]], mode='reflect')
    bottom_row = currents_v[-1:, :]
    currents_v_adjusted = ndimage.generic_filter(scipy.vstack((currents_v, bottom_row)),
                                                 compute_average, footprint=[[1], [1]], mode='reflect')

    # zoom
    #-------------------------------------------------------------------------
    u_zoomed = crop_and_downsample(currents_u_adjusted, downsample_ratio)
    v_zoomed = crop_and_downsample(currents_v_adjusted, downsample_ratio)
    rho_mask[rho_mask == 1] = numpy.nan
    rho_mask_zoomed = crop_and_downsample(rho_mask, downsample_ratio)
    longs = data_file.variables['lon_rho'][:]
    lats = data_file.variables['lat_rho'][:]

    longs_zoomed = crop_and_downsample(longs, downsample_ratio, False)
    lats_zoomed = crop_and_downsample(lats, downsample_ratio, False)

    u_zoomed[rho_mask_zoomed == 1] = numpy.nan
    v_zoomed[rho_mask_zoomed == 1] = numpy.nan

    x, y = bmap(longs_zoomed, lats_zoomed)

    bmap.drawmapboundary(linewidth=0.0, ax=ax)

    overlay = bmap.quiver(x, y, u_zoomed, v_zoomed, ax=ax, color='black', units='inches',
                          scale=10.0, headwidth=2, headlength=3,
                          headaxislength=2.5, minlength=0.5, minshaft=.9)

    # Multiplying .5, 1, and 2 by .5144 is converting from knots to m/s
    #-------------------------------------------------------------------------
    quiverkey = key_ax.quiverkey(overlay, .95, .4, 0.5*.5144, ".5 knots", labelpos='S', labelcolor='white',
                                 color='white', labelsep=.5, coordinates='axes')
    quiverkey1 = key_ax.quiverkey(overlay, 3.75, .4, 1*.5144, "1 knot", labelpos='S', labelcolor='white',
                                  color='white', labelsep=.5, coordinates='axes')
    quiverkey2 = key_ax.quiverkey(overlay, 6.5, .4, 2*.5144, "2 knots", labelpos='S', labelcolor='white',
                                  color='white', labelsep=.5, coordinates='axes')
    key_ax.set_axis_off()

def rag_solidity(labels, connectivity=2):

    graph = RAG()

    # The footprint is constructed in such a way that the first
    # element in the array being passed to _add_edge_filter is
    # the central value.
    fp = ndi.generate_binary_structure(labels.ndim, connectivity)
    for d in range(fp.ndim):
        fp = fp.swapaxes(0, d)
        fp[0, ...] = 0
        fp = fp.swapaxes(0, d)

    # For example
    # if labels.ndim = 2 and connectivity = 1
    # fp = [[0,0,0],
    #       [0,1,1],
    #       [0,1,0]]
    #
    # if labels.ndim = 2 and connectivity = 2
    # fp = [[0,0,0],
    #       [0,1,1],
    #       [0,1,1]]

    ndi.generic_filter(
        labels,
        function=_add_edge_filter,
        footprint=fp,
        mode='nearest',
        output=np.zeros(labels.shape, dtype=np.uint8),
        extra_arguments=(graph,))

    # remove bg_label
    # graph.remove_node(-1)
    graph.remove_node(0)

    for n in graph:
        mask = (labels == n)
        solidity = 1. * mask.sum() / convex_hull_image(mask).sum()
        graph.node[n].update({'labels': [n],
                              'solidity': solidity,
                              'mask': mask})

    for x, y, d in graph.edges_iter(data=True):
        new_mask = np.logical_or(graph.node[x]['mask'], graph.node[y]['mask'])
        new_solidity = 1. * new_mask.sum() / convex_hull_image(new_mask).sum()
        org_solidity = np.mean([graph.node[x]['solidity'],
                                graph.node[y]['solidity']])
        d['weight'] = org_solidity / new_solidity

    return graph

def texture(gray_img, ksize, threshold, offset=3, texture_method='dissimilarity', borders='nearest',
            max_value=255):
    """Creates a binary image from a grayscale image using skimage texture calculation for thresholding.
    This function is quite slow.

    Inputs:
    gray_img       = Grayscale image data
    ksize          = Kernel size for texture measure calculation
    threshold      = Threshold value (0-255)
    offset         = Distance offsets
    texture_method = Feature of a grey level co-occurrence matrix, either
                     'contrast', 'dissimilarity', 'homogeneity', 'ASM', 'energy',
                     or 'correlation'.For equations of different features see
                     scikit-image.
    borders        = How the array borders are handled, either 'reflect',
                     'constant', 'nearest', 'mirror', or 'wrap'
    max_value      = Value to apply above threshold (usually 255 = white)

    Returns:
    bin_img        = Thresholded, binary image

    :param gray_img: numpy.ndarray
    :param ksize: int
    :param threshold: int
    :param offset: int
    :param texture_method: str
    :param borders: str
    :param max_value: int
    :return bin_img: numpy.ndarray
    """

    # Function that calculates the texture of a kernel
    def calc_texture(inputs):
        inputs = np.reshape(a=inputs, newshape=[ksize, ksize])
        inputs = inputs.astype(np.uint8)
        # Greycomatrix takes image, distance offset, angles (in radians), symmetric, and normed
        # http://scikit-image.org/docs/dev/api/skimage.feature.html#skimage.feature.greycomatrix
        glcm = greycomatrix(inputs, [offset], [0], 256, symmetric=True, normed=True)
        diss = greycoprops(glcm, texture_method)[0, 0]
        return diss

    # Make an array the same size as the original image
    output = np.zeros(gray_img.shape, dtype=gray_img.dtype)

    # Apply the texture function over the whole image
    generic_filter(gray_img, calc_texture, size=ksize, output=output, mode=borders)

    # Threshold so higher texture measurements stand out
    bin_img = binary(gray_img=output, threshold=threshold, max_value=max_value, object_type='light')
    return bin_img

def simulate_fire(grid_size, prob_tree, prob_burning, prob_lightning,
                  prob_immune, t):
    grids = []
    grid = init_grid(grid_size, prob_tree, prob_burning)
    grids.append(grid)
    for i in range(t):
        new_grid = np.zeros_like(grid)
        ndimage.generic_filter(grids[-1], spread, size=3, mode="constant",
                               output=new_grid,
                               # these are passed to spread
                               extra_arguments=(prob_immune,
                                                prob_lightning))
        grids.append(new_grid.copy())
    return grids

    def check(j):
        func = FILTER2D_FUNCTIONS[j]

        im = np.ones((20, 20))
        im[:10,:10] = 0
        footprint = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
        footprint_size = np.count_nonzero(footprint)
        weights = np.ones(footprint_size)/footprint_size

        res = ndimage.generic_filter(im, func(weights),
                                     footprint=footprint)
        std = ndimage.generic_filter(im, filter2d, footprint=footprint,
                                     extra_arguments=(weights,))
        assert_allclose(res, std, err_msg="#{} failed".format(j))

def test_generic_filter():
    def filter2d(footprint_elements, weights):
        return (weights*footprint_elements).sum()

    im = np.ones((20, 20))
    im[:10,:10] = 0
    footprint = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
    footprint_size = np.count_nonzero(footprint)
    weights = np.ones(footprint_size)/footprint_size
    for mod in MODULES:
        res = ndimage.generic_filter(im, mod.filter2d(weights),
                                     footprint=footprint)
        std = ndimage.generic_filter(im, filter2d, footprint=footprint,
                                     extra_arguments=(weights,))
        assert_allclose(res, std, err_msg="{} failed".format(mod.__name__))

def background_variance_filter(data, bbox):
    """
    Determine the background variance for each pixel from a box with size of
    bbox.

    Parameters
    ----------
    data : `~numpy.ndarray`
        Data to measure background variance

    bbox :  int
        Box size for calculating background variance

    Raises
    ------
    ValueError
        A value error is raised if bbox is less than 1

    Returns
    -------
    background : `~numpy.ndarray` or `~numpy.ma.MaskedArray`
        An array with the measured background variance in each pixel

    """
    # Check to make sure the background box is an appropriate size
    if bbox < 1:
        raise ValueError('bbox must be greater than 1')

    return ndimage.generic_filter(data, sigma_func, size=(bbox, bbox))

def next_seq(curr, cnr=False):
  old = curr
  footprint = np.array([[1,1,1],
                        [1,0,1],
                        [1,1,1]])

  sums= ndimage.generic_filter(curr, sum, footprint=footprint,
                              mode='constant', cval=0)
  curr = curr.flatten()
  sums = sums.flatten()

  for i in range(len(curr)):
    if curr[i] == 1:
      if sums[i] != 2 and sums[i] != 3:
        curr[i] = 0
    else:
      if sums[i] == 3:
        curr[i] = 1
  curr = curr.reshape(old.shape)
  if cnr:
    curr[0,0] = 1
    curr[0,-1] = 1
    curr[-1,0] = 1
    curr[-1,-1] = 1
  return curr

    def _filter_meshes(self, data_low_res):
        """
        Apply a 2d median filter to the low-resolution background map,
        including only pixels inside the image at the borders.
        """

        from scipy.ndimage import generic_filter
        try:
            nanmedian_func = np.nanmedian    # numpy >= 1.9
        except AttributeError:
            from scipy.stats import nanmedian
            nanmedian_func = nanmedian

        if self.filter_threshold is None:
            return generic_filter(data_low_res, nanmedian_func,
                                  size=self.filter_shape, mode='constant',
                                  cval=np.nan)
        else:
            data_out = np.copy(data_low_res)
            for i, j in zip(*np.nonzero(data_low_res >
                                        self.filter_threshold)):
                yfs, xfs = self.filter_shape
                hyfs, hxfs = yfs // 2, xfs // 2
                y0, y1 = max(i - hyfs, 0), min(i - hyfs + yfs,
                                               data_low_res.shape[0])
                x0, x1 = max(j - hxfs, 0), min(j - hxfs + xfs,
                                               data_low_res.shape[1])
                data_out[i, j] = np.median(data_low_res[y0:y1, x0:x1])
            return data_out

def solve(Z, start, goal):
    Z = 1 - Z
    G = np.zeros(Z.shape)
    G[start] = 1

    # We iterate until value at exit is > 0. This requires the maze
    # to have a solution or it will be stuck in the loop.
    def diffuse(Z, gamma=0.99):
        return max(gamma*Z[0], gamma*Z[1], Z[2], gamma*Z[3], gamma*Z[4])

    G_gamma = np.empty_like(G)
    while G[goal] == 0.0:
        G = Z * generic_filter(G, diffuse, footprint=[[0, 1, 0],
                                                      [1, 1, 1],
                                                      [0, 1, 0]])
    
    # Descent gradient to find shortest path from entrance to exit
    y, x = goal
    dirs = (0,-1), (0,+1), (-1,0), (+1,0)
    P = []
    while (x, y) != start:
        P.append((y,x))
        neighbours = [-1, -1, -1, -1]
        if x > 0:            neighbours[0] = G[y, x-1]
        if x < G.shape[1]-1: neighbours[1] = G[y, x+1]
        if y > 0:            neighbours[2] = G[y-1, x]
        if y < G.shape[0]-1: neighbours[3] = G[y+1, x]
        a = np.argmax(neighbours)
        x, y  = x + dirs[a][1], y + dirs[a][0]
    P.append((y,x))
    return P

    def _filter_meshes(self, mesh2d):
        """
        Apply a 2D median filter to the low-resolution 2D meshes,
        including only pixels inside the image at the borders.
        """

        from scipy.ndimage import generic_filter
        try:
            nanmedian_func = np.nanmedian    # numpy >= 1.9
        except AttributeError:
            from scipy.stats import nanmedian
            nanmedian_func = nanmedian

        if self.filter_threshold is None:
            return generic_filter(mesh2d, nanmedian_func,
                                  size=self.filter_size, mode='constant',
                                  cval=np.nan)
        else:
            # selectively filter only pixels above ``filter_threshold``
            data_out = np.copy(mesh2d)
            for i, j in zip(*np.nonzero(mesh2d > self.filter_threshold)):
                yfs, xfs = self.filter_size
                hyfs, hxfs = yfs // 2, xfs // 2
                y0, y1 = max(i - hyfs, 0), min(i - hyfs + yfs,
                                               mesh2d.shape[0])
                x0, x1 = max(j - hxfs, 0), min(j - hxfs + xfs,
                                               mesh2d.shape[1])
                data_out[i, j] = np.median(mesh2d[y0:y1, x0:x1])
            return data_out

def lnlike(modelpatch,data,sig_smooth,sig_L2,sig_one,w_L2):
    """
    Return the negative log-likelihood given a pixel patch across a
    given set of data patches, weighted by regularization priors.
    Uniform noise case.
    """
    
    # Likelihood given current psf model
    lnlike = 0.0
    for ii in range(data.npatches):
        patch = np.ravel(data.patches[ii])
        flux  = np.dot(modelpatch.T,patch)/np.dot(modelpatch.T,modelpatch)
        model = modelpatch*flux
        lnlike += np.sum(0.5*(patch-model) ** 2
                         / data.bkg_sigmas[ii]**2. + \
                         0.5 * np.log(data.bkg_sigmas[ii]**2.))

    # Smoothness constraint
    if sig_smooth!=0:
            filt = np.array([[False,True,False],
                             [True,True,True],
                             [False,True,False]])
            nearest = ndimage.generic_filter(np.reshape(modelpatch,data.patchshape),
                                     sq_nearest, footprint=filt)
            lnlike  += np.sum(nearest) * sig_smooth

    # L2 norm
    if sig_L2!=0:
        lnlike += np.sum((modelpatch*w_L2)**2.) * sig_L2

    # PSF total ~ 1
    if sig_one!=0:
        lnlike += (np.sum(modelpatch)-1)**2. * sig_one

    return lnlike