def xyz2lab(X, Y, Z, th=0.008856):
    """transform XYZ tristimulus arrays to Lab values

    Parameters
    ----------
    X : np.array, size=(m,n)
    Y : np.array, size=(m,n)
    Z : np.array, size=(m,n)

    Returns
    -------
    L : np.array, size=(m,n)
    a : np.array, size=(m,n)
    b : np.array, size=(m,n)

    See also
    --------
    rgb2xyz, xyz2lms, lms2lch

    Notes
    -----
    .. [1] Ford & Roberts. "Color space conversion", pp. 1--31, 1998.
    .. [2] Silva et al. "Near real-time shadow detection and removal in aerial
       motion imagery application" ISPRS journal of photogrammetry and remote
       sensing, vol.140 pp.104--121, 2018.
    """
    Xn,Yn,Zn = 95.047, 100.00, 108.883 # D65 illuminant

    YYn = Y/Yn

    L_1 = 116* YYn**(1/3.)
    L_2 = 903.3 * YYn

    L = L_1
    L[YYn<=th] = L_2[YYn<=th]

    def f(tau, th):
        fx = X**(1/3.)
        fx[X<=th] = 7.787*X[X<th] + 16/116
        return fx

    a = 500*( f(X/Xn, th) - f(Z/Zn, th) )
    b = 200*( f(Y/Yn, th) - f(Z/Zn, th) )
    return L, a, b