def ssim_global(x, y, k1=0.01, k2=0.03, l=1.0):
    """
    mmsim : https://ece.uwaterloo.ca/~z70wang/publications/msssim.pdf
    wiki : https://en.wikipedia.org/wiki/Structural_similarity
    https://ece.uwaterloo.ca/~z70wang/publications/ssim.pdf
    cf. other metrics
    https://pdfs.semanticscholar.org/401e/c2881b20ce7d6dc07c274a1d8f821d4c3841.pdf

    SSIM(x, y) = ((2 µx µy + c1) / (µx^2 + µy^2 + c1)) * ((2 σxy + c2) / (σx^2 + σy^2 + c2))
    todo :
    similarity between images
    considering illumination, constrast and structure of image
    if c1 = c2 = 0, then SSIM == UQI (universal quality index)

    c1 = (k1 * L)^2, c2 = (k2 * L)^2
    k1 = 0.01 and k2 = 0.03 by default
    L = dynamic range of the pixel values (typically 2^(#bits per per pixel - 1)
    c3 - c2 / 2
    :param x: image [batch, h, w, c]
    :param y: image [batch, h, w, c]
    :param c1: todo : check value in original paper
    :param c2: todo : check value in original paper
    :return: similarity
    """
    reduction_indices = 1, 2, 3
    c1 = (k1 * l) ** 2
    c2 = (k2 * l) ** 2
    mx = x.mean(axis=reduction_indices, keepdims=True)
    my = y.mean(axis=reduction_indices, keepdims=True)
    sx = x.std(axis=reduction_indices, keepdims=True)
    sy = y.std(axis=reduction_indices, keepdims=True)
    xy = (x - mx) * (y - my)
    sxy = xy.mean(axis=reduction_indices, keepdims=True)
    similarity = (2.0 * mx * my + c1) / (mx.square() + my.square() + c1)
    similarity *= (2 * sxy + c2) / (sx.square() + sy.square() + c2)
    return similarity
