import numpy as np


def conf_test(rad, 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 ---------->
    '''

    radshape = rad.shape
    rad = rad.reshape(np.prod(radshape))

    if thr.ndim == 1:
        thr = np.full((rad.shape[0], 4), thr[:4]).T

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

    alpha, gamma = np.empty(rad.shape), np.empty(rad.shape)
    flipped = np.zeros(rad.shape)

#    if thr.ndim == 2:
    gamma[hicut > locut] = thr[0, hicut > locut]
    alpha[hicut > locut] = thr[2, hicut > locut]
    flipped[hicut > locut] = False
    gamma[hicut < locut] = thr[2, hicut < locut]
    alpha[hicut < locut] = thr[0, hicut < locut]
    flipped[hicut < locut] = True
#    else:
#        gamma[hicut > locut] = thr[0]
#        alpha[hicut > locut] = thr[2]
#        flipped[hicut > locut] = False
#        gamma[hicut < locut] = thr[2]
#        alpha[hicut < locut] = thr[0]
#        flipped[hicut < locut] = True

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

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

    # Rad outside alpha-gamma interval
    confidence[(rad > gamma) & (flipped is False)] = 1
    confidence[(rad < alpha) & (flipped is False)] = 0
    confidence[(rad > gamma) & (flipped is True)] = 0
    confidence[(rad < alpha) & (flipped is True)] = 1

    confidence[confidence > 1] = 1
    confidence[confidence < 0] = 0

    return confidence