def get_gradient_y(shape, py):
    """Calculate the gradient in the y direction to the line at py

    The y gradient operator is a block matrix, where each block is the size of the image width.
    The matrix itself is made up of (img_height x img_height) blocks, most of which are all zeros.
    """
    import scipy.sparse

    height, width = shape
    rows = []
    empty = scipy.sparse.dia_matrix((width, width))
    identity = scipy.sparse.identity(width)

    # Create the blocks by row, beginning with blocks leading up to the peak row from the top
    for n in range(py):
        row = [empty]*n
        row += [-identity, identity]
        row += [empty]*(height-n-2)
        rows.append(row)
    # Set all elements in the peak row to zero
    rows.append([empty]*height)
    # Create the blocks for the rows leading up to the peak row from the bottom
    for n in range(height-py-1):
        row = [empty]*(py+n)
        row += [identity, -identity]
        row += [empty]*(height-py-n-2)
        rows.append(row)
    return scipy.sparse.bmat(rows)
Calculate the gradient in the y direction to the line at py

    The y gradient operator is a block matrix, where each block is the size of the image width.
    The matrix itself is made up of (img_height x img_height) blocks, most of which are all zeros.