g2s.py 文件源码

def g2s_weighted(x, y, Zx, Zy, Lxx, Lxy, Lyx, Lyy, N=3):
    """
    % Purpose : Computes the Global Weighted Least Squares reconstruction of a
    %   surface from its gradient field, whereby the weighting is defined by a
    %   weighted Frobenius norm
    %
    % Use (syntax):
    %   Z = g2sWeighted( Zx, Zy, x, y, N, Lxx, Lxy, Lyx, Lyy )
    %
    % Input Parameters :
    %   Zx, Zy := Components of the discrete gradient field
    %   x, y := support vectors of nodes of the domain of the gradient
    %   N := number of points for derivative formulas (default=3)
    %   Lxx, Lxy, Lyx, Lyy := Each matrix Lij is the covariance matrix of the
    %       gradient's i-component the in j-direction.
    %
    % Return Parameters :
    %   Z := The reconstructed surface
    %
    % Description and algorithms:
    %   The algorithm solves the normal equations of the Weighted Least Squares
    %   cost function, formulated by matrix algebra:
    %   e(Z) = || Lyy^(-1/2) * (D_y * Z - Zy) * Lyx^(-1/2) ||_F^2 +
    %               || Lxy^(-1/2) * ( Z * Dx' - Zx ) * Lxx^(-1/2) ||_F^2
    %   The normal equations are a rank deficient Sylvester equation which is
    %   solved by means of Householder reflections and the Bartels-Stewart
    %   algorithm.
    """

    if Zx.shape != Zy.shape:
        raise ValueError("Gradient components must be the same size")

    if np.asmatrix(Zx).shape[1] != len(x) or np.asmatrix(Zx).shape[0] != len(y):
        raise ValueError("Support vectors must have the same size as the gradient")

    m, n = Zx.shape

    Dx = diff_local(x, N, N)
    Dy = diff_local(y, N, N)

    Wxx = la.sqrtm(Lxx)
    Wxy = la.sqrtm(Lxy)
    Wyx = la.sqrtm(Lyx)
    Wyy = la.sqrtm(Lyy)

    # Solution for Zw (written here Z)
    u = mldivide(Wxy, np.ones((m, 1)))
    v = mldivide(Wyx, np.ones((n, 1)))

    A = mldivide(Wyy, np.dot(Dy, Wxy))
    B = mldivide(Wxx, np.dot(Dx, Wyx))
    F = mldivide(Wyy, mrdivide(Zy, Wyx))
    G = mldivide(Wxy, mrdivide(Zx, Wxx))

    Z = _g2sSylvester(A, B, F, G, u, v)

    # "Unweight" the solution
    Z = Wxy * Z * Wyx
    return Z