Fit a bidimensional polynomial to an image.

    :param data: a bidimensional array
    :param integer order: order of the polynomial
    :param bool output_fit: return the fitted image
    :returns: a tuple with an array with the coefficients of the polynomial terms

    >>> import numpy
    >>> xx, yy = numpy.mgrid[-1:1:100j,-1:1:100j]
    >>> z = 456.0 + 0.3 * xx - 0.9* yy
    >>> imsurfit(z, order=1) #doctest: +NORMALIZE_WHITESPACE
    (array([  4.56000000e+02,   3.00000000e-01,  -9.00000000e-01]),)

def imsurfit(data, order, output_fit=False):
    """Fit a bidimensional polynomial to an image.

    :param data: a bidimensional array
    :param integer order: order of the polynomial
    :param bool output_fit: return the fitted image
    :returns: a tuple with an array with the coefficients of the polynomial terms

    >>> import numpy
    >>> xx, yy = numpy.mgrid[-1:1:100j,-1:1:100j]
    >>> z = 456.0 + 0.3 * xx - 0.9* yy
    >>> imsurfit(z, order=1) #doctest: +NORMALIZE_WHITESPACE
    (array([  4.56000000e+02,   3.00000000e-01,  -9.00000000e-01]),)

    """

    # we create a grid with the same number of points
    # between -1 and 1
    c0 = complex(0, data.shape[0])
    c1 = complex(0, data.shape[1])
    xx, yy = numpy.ogrid[-1:1:c0, -1:1:c1]

    ncoeff = (order + 1) * (order + 2) // 2

    powerlist = list(_powers(order))

    # Array with ncoff x,y moments of the data
    bb = numpy.zeros(ncoeff)
    # Moments
    for idx, powers in enumerate(powerlist):
        p1, p2 = powers
        bb[idx] = (data * xx ** p1 * yy ** p2).sum() / data.size

    # Now computing aa matrix
    # it contains \sum x^a y^b
    # most of the terms are zero
    # only those with a and b even remain
    # aa is symmetric

    x = xx[:, 0]
    y = yy[0]

    # weights are stored so we compute them only once
    wg = {(0, 0): 1, (-1, 0): x**2, (0, -1): y**2}
    # wg[(2,0)] = wg[(-1,0)].sum() / len(x)
    # wg[(0,2)] = wg[(0,-1)].sum() / len(y)

    aa = numpy.zeros((ncoeff, ncoeff))
    for j, ci in enumerate(powerlist):
        for i, ri in enumerate(powerlist[j:]):
            p1 = ci[0] + ri[0]
            p2 = ci[1] + ri[1]
            if p1 % 2 == 0 and p2 % 2 == 0:
                # val = (x ** p1).sum() / len(x) * (y ** p2).sum() / len(y)
                val = _compute_weight((p1, p2), wg)
                aa[j, i+j] = val

    # Making symmetric the array
    aa += numpy.triu(aa, k=1).T

    polycoeff = numpy.linalg.solve(aa, bb)

    if output_fit:
        index = 0
        result = 0
        for o in range(order + 1):
            for b in range(o + 1):
                a = o - b
                result += polycoeff[index] * (xx ** a) * (yy ** b)
                index += 1

        return polycoeff, result

    return (polycoeff,)