def xyz_to_cielab (xyz, refwhite):
    """Convert CIE XYZ color values to CIE L*a*b*.

    *xyz* should be of shape (*, 3). *refwhite* is the reference white value, of
    shape (3, ).

    Return value will have same shape as *xyz*, but be in CIE L*a*b*
    coordinates.

    """
    norm = xyz / refwhite
    pow = norm**0.333333333333333
    scale = 7.787037 * norm + 16./116
    mapped = np.where (norm > 0.008856, pow, scale)

    cielab = np.empty_like (xyz)
    cielab[...,L] = 116 * mapped[...,Y] - 16
    cielab[...,A] = 500 * (mapped[...,X] - mapped[...,Y])
    cielab[...,B] = 200 * (mapped[...,Y] - mapped[...,Z])

    return cielab