import numpy as np


_test_rad = np.random.randint(25, size=[6, 8])
#_test_thr = [15, 10, 5, 1, 1]
_test_thr = [5, 10, 15, 1, 1]


def test_11um(rad, threshold):

    #if (~np.isnan(rad) or threshold[4] == 1.0):
    confidence = conf_test(rad, threshold)

    return confidence


def test_11_4diff(rad1, rad2, threshold):

    raddiff = rad1 - rad2;
    confidence = conf_test(raddiff, threshold)

    return confidence


def simple_threshold_test(rad, threshold):
    return conf_test(rad, threshold)


class CloudMaskTests(object):

    def __init__(self, scene, radiance, coefficients):
        self.scene = scene
        self.coefficients = coefficients

    def select_coefficients(self):
        pass

    def test_G1(self):
        pass

    def test_G2(self):
        pass

    def test_G3(self):
        pass

    def test_G4(self):
        pass

    def overall_confidence(self):
        pass


def conf_test(rad=_test_rad, thr=_test_thr):
    '''
    Assuming a linear function between min and max confidence level, the plot below shows
    how the confidence (y axis) is computed as function of radiance (x axis).
    This case illustrates alpha < gamma, obviously in case alpha > gamma, the plot would be
    flipped.
                       gamma
    c  1                 ________
    o  |                /
    n  |               /
    f  |              /
    i  |     beta    /
    d 1/2    |....../
    e  |           /
    n  |          /
    c  |         /
    e  0________/
       |      alpha
    --------- radiance ---------->
    '''

    coeff = np.power(2, (thr[3] - 1))
    hicut = thr[0]
    beta = thr[1]
    locut = thr[2]
    power = thr[3]
    radshape = rad.shape
    rad = rad.reshape((rad.shape[0]*rad.shape[1]))
    c = np.zeros(rad.shape)

    if hicut > locut:
        gamma = thr[0]
        alpha = thr[2]
        flipped = False
    else:
        gamma = thr[2]
        alpha = thr[0]
        flipped = True

    # Rad between alpha and beta
    range_ = 2. * (beta - alpha)
    s1 = (rad[rad <= beta] - alpha) / range_
    if flipped is False:
        c[rad <= beta] = coeff * np.power(s1, power)
    if flipped is True:
        c[rad <= beta] = 1. - (coeff * np.power(s1, power))

    # Rad between beta and gamma
    range_ = 2. * (beta - gamma)
    s1 = (rad[rad > beta] - gamma) / range_
    if flipped is False:
        c[rad > beta] = 1. - (coeff * np.power(s1, power))
    if flipped is True:
        c[rad > beta] = coeff * np.power(s1, power)

    # Rad outside alpha-gamma interval
    if flipped is False:
        c[rad > gamma] = 1
        c[rad < alpha] = 0
    if flipped is True:
        c[rad > gamma] = 0
        c[rad < alpha] = 1

    c[c > 1] = 1
    c[c < 0] = 0

    confidence = c.reshape(radshape)

    return confidence


def conf_test_dble(rad, coeffs):
    '''

            gamma1                         gamma2
    c  1_______                               ________
    o  |       \                             /
    n  |        \                           /
    f  |         \                         /
    i  |          \   beta1       beta2   /
    d 1/2          \....|          |...../
    e  |            \                   /
    n  |             \                 /
    c  |              \               /
    e  0               \_____________/
       |             alpha1       alpha2
    --------------------- radiance ------------------------->
    '''

    hicut = [coeffs[0], coeffs[1]]
    locut = [coeffs[2], coeffs[3]]
    midpt = [coeffs[4], coeffs[5]]
    power = coeffs[6]

    gamma1 = hicut[0]
    gamma2 = hicut[1]
    alpha1 = locut[0]
    alpha2 = locut[1]
    beta1 = midpt[0]
    beta2 = midpt[1]

    coeff = np.power(2, (power - 1))
    radshape = rad.shape
    rad = rad.reshape((rad.shape[0]*rad.shape[1]))
    c = np.zeros(rad.shape)

    # Find if interval between inner cutoffs passes or fails test
    if (alpha1 - gamma1 > 0):

       # Value is within range of lower set of limits
        range_ = 2 * (beta1 - alpha1)
        s1 = (rad[(rad <= alpha1) & (rad >= beta1)] - alpha1) / range_
        c[(rad <= alpha1) & (rad >= beta1)] = coeff * np.power(s1, power)

        range_ = 2 * (beta1 - gamma1)
        s1 = (rad[(rad <= alpha1) & (rad < beta1)] - gamma1) / range_
        c[(rad <= alpha1) & (rad < beta1)] = coeff * np.power(s1, power)

        # Value is within range of upper set of limits
        range_ = 2 * (beta2 - alpha2)
        s1 = (rad[(rad > alpha1) & (rad <= beta2)] - alpha2) / range_
        c[(rad > alpha1) & (rad <= beta2)] = coeff * np.power(s1, power)

        range_ = 2 * (beta2 - gamma2)
        s1 = (rad[(rad > alpha1) & (rad > beta2)] - gamma2) / range_
        c[(rad > alpha1) & (rad > beta2)] = coeff * np.power(s1, power)

        # Inner region fails test
        # Check for value beyond function range
        c[(rad > alpha1) & (rad < alpha2)] = 0
        c[(rad < gamma1) | (rad > gamma2)] = 1

    else:

        # Value is withing range of lower set of limits
        range_ = 2 * (beta1 - alpha1)
        s1 = (rad[(rad <= gamma1) & (rad <= beta1)] - alpha1) / range_
        c[(rad <= gamma1) & (rad <= beta1)] = coeff * np.power(s1, power)

        range_ = 2 * (beta1 - gamma1)
        s1 = (rad[(rad <= gamma1) & (rad > beta1)] - gamma1) / range_
        c[(rad <= gamma1) & (rad > beta1)] = coeff * np.power(s1, power)

        # Value is within range of upper set of limits
        range_ = 2 * (beta2 - alpha2)
        s1 = (rad[(rad > gamma1) & (rad >= beta2)] - alpha2) / range_
        c[(rad > gamma1) & (rad >= beta2)] = coeff * np.power(s1, power)

        range_ = 2 * (beta2 - gamma2)
        s1 = (rad[(rad > gamma1) & (rad < beta2)] - gamma2) / range_
        c[(rad > gamma1) & (rad < beta2)] = coeff * np.power(s1, power)

        # Inner region passes test
        # Check for value beyond function range
        c[(rad > gamma1) & (rad < gamma2)] = 1
        c[(rad < alpha1) | (rad > alpha2)] = 0

    c[c>1] = 1
    c[c<0] = 0

    confidence = c.reshape(radshape)

    return confidence


def test():
    rad = np.random.randint(50, size=[4, 8])
    #coeffs = [5, 42, 20, 28, 15, 35, 1]
    #coeffs = [20, 28, 5, 42, 15, 35, 1]
    coeffs = [35, 15, 20, 1, 1]
    #confidence = conf_test_dble(rad, coeffs)
    confidence = test_11um(rad, coeffs)
    print(rad)
    print('\n')
    print(confidence)


if __name__ == "__main__":
    test()