using ColorSchemes, Colors, Interpolations

# code from peterkovesi/PerceptualColourMaps.jl
# https://github.com/peterkovesi/PerceptualColourMaps.jl/
# Peter is no longer updating his package
# this code of his provides equalizecolormap

#-------------------------------------------------------------------------------
"""
equalisecolourmap/equalizecolormap - Equalise colour contrast over a colour map
```
Usage: newrgbmap = equalisecolourmap(rgblab, map, formula, W, sigma, diagnostics)
                   equalizecolormap(....
Arguments:     rgblab - String "RGB" or "LAB" indicating the type of data
                        in map.
                  map - A Nx3 RGB or CIELAB colour map
                        or an array of ColorTypes.RGBA{Float64} values
              formula - String "CIE76" or "CIEDE2000"
                    W - A 3-vector of weights to be applied to the
                        lightness, chroma and hue components of the
                        difference equation. It is recommended that you
                        use [1, 0, 0] to only take into account lightness.
                        If desired use  [1, 1, 1] for the full formula.
                        See note below.
                sigma - Optional Gaussian smoothing parameter, see
                        explanation below.
               cyclic - Boolean flag indicating whether the colour map is
                        cyclic. This affects how smoothing is applied at
                        the end points.
          diagnostics - Optional boolean flag indicating whether diagnostic
                        plots should be displayed.  Defaults to false.
Returns:    newrgbmap - RGB colour map adjusted so that the perceptual
                        contrast of colours along the colour map is constant.
                        This is a Nx3 Array of Float64 values.
```
Suggested parameters:
The CIE76 and CIEDE2000 colour difference formulas were developed for
much lower spatial frequencies than we are typically interested in.
Neither is ideal for our application.  The main thing to note is that
at *fine* spatial frequencies perceptual contrast is dominated by
*lightness* difference, chroma and hue are relatively unimportant.
For colour maps with a significant range of lightness use:
```
                       formula = "CIE76" or "CIEDE2000"
                             W = [1, 0, 0]  (Only correct for lightness)
                         sigma = 5 - 7
```
For isoluminant or low lightness gradient colour maps use:
```
                       formula = "CIE76"
                             W = [1, 1, 1]  (Correct for colour and lightness)
                         sigma = 5 - 7
```
Ideally, for a colour map to be effective the perceptual contrast along the
colour map should be constant.  Many colour maps are very poor in this regard.
Try testing your favourite colour map on the sineramp() test image.  The
perceptual contrast is very much dominated by the contrast in colour lightness
values along the map.  This function attempts to equalise the chosen
perceptual contrast measure along a colour map by stretching and/or
compressing sections of the colour map.
This function's primary use is for the correction of colour maps generated by
cmap() however it can be applied to any colour map.  There are limitations to
what this function can correct.  When applied to some of MATLAB's colour maps
such as 'jet', 'hsv' and 'cool' you get colour discontinuity artifacts because
these colour maps have segments that are nearly constant in lightness.
However, it does a nice job of fixing up MATLAB's 'hot', 'winter', 'spring'
and 'autumn' colour maps.  If you do see colour discontinuities in the
resulting colour map try changing W from [1, 0, 0] to [1, 1, 1], or some
intermediate weighting of [1, 0.5, 0.5], say.
Difference formula: Neither CIE76 or CIEDE2000 difference measures are ideal
for the high spatial frequencies that we are interested in.  Empirically I
find that CIEDE2000 seems to give slightly better results on colour maps where
there is a significant lightness gradient (this applies to most colour maps).
In this case you would be using a weighting vector W = [1, 0, 0].  For
isoluminant, or low lightness gradient colour maps where one is using a
weighting vector W = [1, 1, 1] CIE76 should be used as the CIEDE2000 chroma
correction is inapropriate for the spatial frequencies we are interested in.
Weighting vetor W: The CIEDE2000 colour difference formula incorporates the
scaling parameters kL, kC, kH in the demonimator of the lightness, chroma, and
hue difference components respectively.  The 3 components of W correspond to
the reciprocal of these 3 parameters.  (I do not know why they chose to put
kL, kC, kH in the denominator. If you wanted to ignore, say, the chroma
component you would have to set kC to Inf, rather than setting W[2] to 0 which
seems more sensible to me).  If you are using CIE76 then W[2] amd W[3] are
applied to the differences in a and b.  In this case you should ensure W[2] =
W[3].  In general, for the spatial frequencies of interest to us, lightness
differences are overwhelmingly more important than chroma or hue and W shoud
be set to [1, 0, 0]
Smoothing parameter sigma:
The output colour map will have lightness values of constant slope magnitude.
However, it is possible that the sign of the slope may change, for example at
the mid point of a bilateral colour map.  This slope discontinuity of lightness
can induce a false apparent feature in the colour map.  A smaller effect is
also occurs for slope discontinuities in a and b.  For such colour maps it can
be useful to introduce a small amount of smoothing of the Lab values to soften
the transition of sign in the slope to remove this apparent feature.  However
in doing this one creates a small region of suppressed luminance contrast in
the colour map which induces a 'blind spot' that compromises the visibility of
features should they fall in that data range.  Accordingly the smoothing
should be kept to a minimum.  A value of sigma in the range 5 to 7 in a 256
element colour map seems about right.  As a guideline sigma should not be more
than about 1/25 of the number of entries in the colour map, preferably less.
Reference: Peter Kovesi. Good Colour Maps: How to Design
Them. [arXiv:1509.03700 [cs.GR] 2015](https://arXiv:1509.03700)
See also: cmap, applycycliccolourmap, applydivergingcolourmap,
sineramp, circlesineramp
"""
function equalisecolourmap(rgblab::AbstractString, cmap::AbstractMatrix{Float64},
                           formula::AbstractString="CIE76", W::Array=[1.0, 0.0, 0.0],
                           sigma::Real = 0.0, cyclic::Bool = false, diagnostics::Bool = false)
   # October  2015  - Ported from MATLAB to Julia

    N = size(cmap,1)   # No of colour map entries

    if N/sigma < 25
        @warn "It is not recommended that sigma be larger than 1/25 of the colour map length"
    end

    rgblab = uppercase(rgblab)
    formula = uppercase(formula)

    if rgblab == "RGB" && (maximum(cmap) > 1.01 || minimum(cmap) < -0.01)
        error("If map is RGB values should be in the range 0-1")
    elseif rgblab == "LAB" && maximum(abs.(cmap)) < 10
        error("If map is LAB magnitude of values are expected to be > 10")
    end

    # If input is RGB convert colour map to Lab. Also, ensure that we
    # have both RGB and Lab representations of the colour map.  I am
    # assuming that the Colors.convert() function uses a default white
    # point of D65
    if rgblab == "RGB"
        rgbmap = copy(cmap)
        labmap = srgb2lab(cmap)
        L = labmap[:,1]
        a = labmap[:,2]
        b = labmap[:,3]
    elseif rgblab == "LAB"
        labmap = copy(cmap)
        rgbmap = lab2srgb(cmap)
        L = cmap[:,1]
        a = cmap[:,2]
        b = cmap[:,3]
    else
        error("Input must be RGB or LAB")
    end

    # The following section of code computes the locations to interpolate into
    # the colour map in order to achieve equal steps of perceptual contrast.
    # The process is repeated recursively on its own output. This helps overcome
    # the approximations induced by using linear interpolation to estimate the
    # locations of equal perceptual contrast. This is mainly an issue for
    # colour maps with only a few entries.

    initialdeltaE = 0  # Define these variables in main scope
    initialcumdE = 0
    initialequicumdE = 0
    initialnewN = 0

    for iter = 1:3
        # Compute perceptual colour difference values along the colour map using
        # the chosen formula and weighting vector.
        if formula ==  "CIE76"
            deltaE = cie76(L, a, b, W)
        elseif formula == "CIEDE2000"
            deltaE = ciede2000(L, a, b, W)
        else
            error("Unknown colour difference formula")
        end

        # Form cumulative sum of of delta E values.  However, first ensure all
        # values are larger than 0.001 to ensure the cumulative sum always
        # increases.
        deltaE[deltaE .< 0.001] .= 0.001
        cumdE = cumsum(deltaE, dims=1)

        # Form an array of equal steps in cumulative contrast change.
        equicumdE =  collect(0:(N-1))./(N-1) .* (cumdE[end]-cumdE[1]) .+ cumdE[1]

        # Solve for the locations that would give equal Delta E values.
        newN = interp1(cumdE, 1:N, equicumdE)

        # newN now represents the locations where we want to interpolate into the
        # colour map to obtain constant perceptual contrast
        Li = interpolate(L, BSpline(Linear()))
        L = [Li(v) for v in newN]

        ai = interpolate(a, BSpline(Linear()))
        a = [ai(v) for v in newN]

        bi = interpolate(b, BSpline(Linear()))
        b = [bi(v) for v in newN]

        # Record initial colour differences for evaluation at the end
        if iter == 1
            initialdeltaE = deltaE
            initialcumdE = cumdE
            initialequicumdE = equicumdE
            initialnewN = newN
        end
    end

    # Apply smoothing of the path in CIELAB space if requested.  The aim is to
    # smooth out sharp lightness/colour changes that might induce the perception
    # of false features.  In doing this there will be some cost to the
    # perceptual contrast at these points.
    if sigma > 0.0
        L = smooth(L, sigma, cyclic)
        a = smooth(a, sigma, cyclic)
        b = smooth(b, sigma, cyclic)
    end

    # Convert map back to RGB
    newlabmap = [L a b]
    newrgbmap = lab2srgb(newlabmap)

    return newrgbmap
end


# Case when colour map is an array of ColorTypes.RGB{Float64}

function equalisecolourmap(rgblab::AbstractString, cmap::AbstractVector{ColorTypes.RGB{Float64}},
                           formula::AbstractString="CIE76", W::Array=[1.0, 0.0, 0.0],
                           sigma::Real = 0.0, cyclic::Bool = false, diagnostics::Bool = false)

    return equalisecolourmap(rgblab, RGB2FloatArray(cmap), formula, W,
                             sigma, cyclic, diagnostics)
end


# Case when colour map is an array of ColorTypes.RGBA{Float64}

function equalisecolourmap(rgblab::AbstractString, cmap::AbstractVector{ColorTypes.RGBA{Float64}},
                           formula::AbstractString="CIE76", W::Array=[1.0, 0.0, 0.0],
                           sigma::Real = 0.0, cyclic::Bool = false, diagnostics::Bool = false)

    return equalisecolourmap(rgblab, RGBA2FloatArray(cmap), formula, W,
                             sigma, cyclic, diagnostics)
end


# Convenience functions for those who spell colour without a 'u' and equalise with a 'z' ...
"""
equalisecolourmap - Equalise colour contrast over a colourmap
equalizecolormap
```
Usage: newrgbmap = equalisecolourmap(rgblab, map, formula, W, sigma, diagnostics)
                   equalizecolormap(....
Arguments:     rgblab - String "RGB" or "LAB" indicating the type of data
                        in map.
                  map - A Nx3 RGB or CIELAB colour map
                        or an array of ColorTypes.RGB{Float64} values
              formula - String "CIE76" or "CIEDE2000"
                    W - A 3-vector of weights to be applied to the
                        lightness, chroma and hue components of the
                        difference equation. It is recommended that you
                        use [1, 0, 0] to only take into account lightness.
                        If desired used  [1, 1, 1] for the full formula.
                sigma - Optional Gaussian smoothing parameter.
               cyclic - Boolean flag indicating whether the colour map is
                        cyclic. This affects how smoothing is applied at
                        the end points.
          diagnostics - Optional boolean flag indicating whether diagnostic
                        plots should be displayed.  Defaults to false.
Returns:    newrgbmap - RGB colour map adjusted so that the perceptual
                        contrast of colours along the colour map is constant.
                        This is a Nx3 Array of Float64 values.
For full documentation see equalisecolourmap()
                                 ^     ^
```
See also: cmap, sineramp, circlesineramp
"""
function equalizecolormap(rgblab::AbstractString, cmap::Array{Float64,2},
                           formula::AbstractString="CIE76", W::Array=[1.0, 0.0, 0.0],
                           sigma::Real = 0.0, cyclic::Bool = false, diagnostics::Bool = false)

    return equalisecolourmap(rgblab, cmap, formula, W, sigma, cyclic, diagnostics)
end

# Case when colour map is an array of ColorTypes.RGB{Float64}

function equalizecolormap(rgblab::AbstractString, cmap::Array{ColorTypes.RGB{Float64},1},
                           formula::AbstractString="CIE76", W::Array=[1.0, 0.0, 0.0],
                           sigma::Real = 0.0, cyclic::Bool = false, diagnostics::Bool = false)

    return equalisecolourmap(rgblab, RGB2FloatArray(cmap), formula, W,
                             sigma, cyclic, diagnostics)
end


# Case when colour map is an array of ColorTypes.RGBA{Float64}

function equalizecolormap(rgblab::AbstractString, cmap::Array{ColorTypes.RGBA{Float64},1},
                           formula::AbstractString="CIE76", W::Array=[1.0, 0.0, 0.0],
                           sigma::Real = 0.0, cyclic::Bool = false, diagnostics::Bool = false)

    return equalisecolourmap(rgblab, RGBA2FloatArray(cmap), formula, W,
                             sigma, cyclic, diagnostics)
end


#----------------------------------------------------------------------------
#
# Function to smooth an array of values but also ensure end values are
# not altered or, if the map is cyclic, ensures smoothing is applied
# across the end points in a cyclic manner.  It is assumed that the
# input data is a column vector.

function smooth(L::Array{T,1}, sigma::Real, cyclic::Bool) where {T<:Real}

    if cyclic
        Le = [L; L; L] # Form a concatenation of 3 repetitions of the array.

        Ls = gaussfilt1d(Le, sigma)               # Apply smoothing filter
        Ls = Ls[length(L)+1: length(L)+length(L)] # and then return the center section

    else  # Non-cyclic colour map: Pad out input array L at both ends by 3*sigma
          # with additional values at the same slope.  The aim is to eliminate
          # edge effects in the filtering
        extension = collect(1:ceil(3*sigma))

        dL1 = L[2]-L[1]
        dL2 = L[end]-L[end-1]
        Le = [-reverse(dL1*extension,dims=1).+L[1]; L;  dL2*extension.+L[end]]

        Ls = gaussfilt1d(Le, sigma) # Apply smoothing filter

        # Trim off extensions
        Ls = Ls[length(extension)+1 : length(extension)+length(L)]
    end

    return Ls
end

#----------------------------------------------------------------------------
"""
deltaE: Compute weighted Delta E between successive entries in a
colour map using the CIE76 formula + weighting
```
Usage: deltaE = cie76(L::Array, a::Array, b::Array, W::Array)
```
"""
function cie76(L::Array, a::Array, b::Array, W::Array)

    N = length(L)

    # Compute central differences
    dL = zeros(size(L))
    da = zeros(size(a))
    db = zeros(size(b))

    dL[2:end-1] = (L[3:end] - L[1:end-2])/2
    da[2:end-1] = (a[3:end] - a[1:end-2])/2
    db[2:end-1] = (b[3:end] - b[1:end-2])/2

    # Differences at end points
    dL[1] = L[2] - L[1];  dL[end] = L[end] - L[end-1]
    da[1] = a[2] - a[1];  da[end] = a[end] - a[end-1]
    db[1] = b[2] - b[1];  db[end] = b[end] - b[end-1]

    return deltaE = sqrt.(W[1]*dL.^2 + W[2]*da.^2 + W[3]*db.^2)
end

#----------------------------------------------------------------------------
"""
ciede2000: Compute weighted Delta E between successive entries in a
colour map using the CIEDE2000 formula + weighting
```
Usage: deltaE = ciede2000(L::Array, a::Array, b::Array, W::Array)
```
"""
function ciede2000(L::Array, a::Array, b::Array, W::Array)

    N = length(L)
    deltaE = zeros(N, 1)
    kl = 1/W[1]
    kc = 1/W[2]
    kh = 1/W[3]

    # Compute deltaE using central differences
    for i = 2:N-1
        deltaE[i] = Colors.colordiff(Colors.Lab(L[i+1],a[i+1],b[i+1]), Colors.Lab(L[i-1],a[i-1],b[i-1]);
                                     metric=Colors.DE_2000(kl,kc,kh))/2
    end

    # Differences at end points
    deltaE[1] = Colors.colordiff(Colors.Lab(L[2],a[2],b[2]), Colors.Lab(L[1],a[1],b[1]);
                                 metric = Colors.DE_2000(kl,kc,kh))
    deltaE[N] = Colors.colordiff(Colors.Lab(L[N],a[N],b[N]), Colors.Lab(L[N-1],a[N-1],b[N-1]);
                                 metric=Colors.DE_2000(kl,kc,kh))

    return deltaE
end


#----------------------------------------------------------------------------
"""
Convenience function for converting an Nx3 array of RGB values in a
colour map to an Nx3 array of CIELAB values.  Function can also be
used to convert a 3 channel RGB image to a 3 channel CIELAB image
Note it appears that the Colors.convert() function uses a default white
point of D65
```
 Usage:  lab = srgb2lab(rgb)
 Argument:    rgb - A N x 3 array of RGB values or a 3 channel RGB image.
 Returns:     lab - A N x 3 array of Lab values of a 3 channel CIELAB image.
```
See also: lab2srgb
"""
function srgb2lab(rgb::AbstractMatrix{T}) where {T}

    N = size(rgb,1)
    lab = zeros(N,3)

    for i = 1:N
        labval = Colors.convert(ColorTypes.Lab, ColorTypes.RGB(rgb[i,1], rgb[i,2], rgb[i,3]))
        lab[i,1] = labval.l
        lab[i,2] = labval.a
        lab[i,3] = labval.b
    end

    return lab
end

#----------------------------------------------------------------------------
#
# Convenience function for converting a 3 channel RGB image to a 3
# channel CIELAB image
#
# Usage:  lab = srgb2lab(rgb)

function srgb2lab(rgb::Array{T,3}) where {T}

    (rows, cols, chan) = size(rgb)
    lab = zeros(size(rgb))

    for r = 1:rows, c = 1:cols
        labval = Colors.convert(ColorTypes.Lab, ColorTypes.RGB(rgb[r,c,1], rgb[r,c,2], rgb[r,c,3]))
        lab[r,c,1] = labval.l
        lab[r,c,2] = labval.a
        lab[r,c,3] = labval.b
    end

    return lab
end

#----------------------------------------------------------------------------
"""
Convenience function for converting an Nx3 array of CIELAB values in a
colour map to an Nx3 array of RGB values.  Function can also be
used to convert a 3 channel CIELAB image to a 3 channel RGB image
Note it appears that the Colors.convert() function uses a default white
point of D65
```
 Usage:  rgb = srgb2lab(lab)
 Argument:   lab - A N x 3 array of CIELAB values of a 3 channel CIELAB image.
 Returns:    rgb - A N x 3 array of RGB values or a 3 channel RGB image.
```
See also: srgb2lab
"""
function lab2srgb(lab::AbstractMatrix{T}) where {T}

    N = size(lab,1)
    rgb = zeros(N,3)

    for i = 1:N
        rgbval = Colors.convert(ColorTypes.RGB, ColorTypes.Lab(lab[i,1], lab[i,2], lab[i,3]))
        rgb[i,1] = rgbval.r
        rgb[i,2] = rgbval.g
        rgb[i,3] = rgbval.b
    end

    return rgb
end

#----------------------------------------------------------------------------
#
# Convenience function for converting a 3 channel Lab image to a 3
# channel RGB image
#
# Usage:  rgb = lab2srgb(lab)

function lab2srgb(lab::Array{T,3}) where {T}

    (rows, cols, chan) = size(lab)
    rgb = zeros(size(lab))

    for r = 1:rows, c = 1:cols
        rgbval = Colors.convert(ColorTypes.RGB, ColorTypes.Lab(lab[r,c,1], lab[r,c,2], lab[r,c,3]))
        rgb[r,c,1] = rgbval.r
        rgb[r,c,2] = rgbval.g
        rgb[r,c,3] = rgbval.b
    end

    return rgb
end

#----------------------------------------------------------------------------
"""
Convert an array of ColorTypes RGB values to an array of UInt32 values
for use as a colour map in Winston
```
 Usage:  uint32rgb = RGBA2UInt32(rgbmap)
 Argument:   rgbmap - Vector of ColorTypes.RGBA values as
                      returned by cmap().
 Returns: uint32rgb - An array of UInt32 values packed with the 8 bit RGB values.
```
See also: cmap
"""
function RGBA2UInt32(rgb::Vector{ColorTypes.RGBA{Float64}})

    N = length(rgb)
    uint32rgb = zeros(UInt32, N)

    for i = 1:N
        r = round(UInt32, rgb[i].r*255)
        g = round(UInt32, rgb[i].g*255)
        b = round(UInt32, rgb[i].b*255)
        uint32rgb[i] = r << 16 + g << 8 + b
    end

    return uint32rgb
end

#----------------------------------------------------------------------------
"""
linearrgbmap: Linear rgb colourmap from black to a specified colour
```
Usage: cmap = linearrgbmap(C, N)
Arguments:  C - 3-vector specifying RGB colour
            N - Number of colourmap elements, defaults to 256
Returns: cmap - N element ColorTypes.RGBA colourmap ranging from [0 0 0]
                to RGB colour C
```
It is suggested that you pass the resulting colour map to equalisecolourmap()
to obtain a map with uniform steps in perceptual lightness
```
> cmap = equalisecolourmap("rgb", linearrgbmap(C, N))
```
See also: equalisecolourmap, ternarymaps
"""
function linearrgbmap(C::Array, N::Int = 256)

    @assert length(C) == 3 "Colour must be a 3 element array"

    rgbmap = zeros(N,3)
    ramp = (0:(N-1))/(N-1)

    for n = 1:3
        rgbmap[:,n] = C[n] * ramp
    end

    return FloatArray2RGBA(rgbmap)
end

#-------------------------------------------------------------------
# Convert Nx3 Float64 array to  N array of ColorTypes.RGB{Float64}

function FloatArray2RGB(cmap::Array{Float64,2})

    (N,cols) = size(cmap)
    @assert cols == 3  "Color map data must be N x 3"

    rgbmap = Array{ColorTypes.RGB{Float64}}(N)
    for i = 1:N
        rgbmap[i] = ColorTypes.RGB(cmap[i,1], cmap[i,2], cmap[i,3])
    end

    return rgbmap
end

#-------------------------------------------------------------------
# Convert Nx3 Float64 array to  N array of ColorTypes.RGBA{Float64}

function FloatArray2RGBA(cmap::Array{Float64,2})

    (N,cols) = size(cmap)
    @assert cols == 3  "Color map data must be N x 3"

    rgbmap = Array{ColorTypes.RGBA{Float64}}(undef, N)
    for i = 1:N
        rgbmap[i] = ColorTypes.RGBA(cmap[i,1], cmap[i,2], cmap[i,3],1.0)
    end

    return rgbmap
end


#-------------------------------------------------------------------
# Convert N array of ColorTypes.RGB{Float64} to Nx3 Float64 array

function RGB2FloatArray(rgbmap::Array{ColorTypes.RGB{Float64},1})

    N = length(rgbmap)

    cmap = Array{Float64}(undef,N,3)
    for i = 1:N
        cmap[i,:] = [rgbmap[i].r rgbmap[i].g rgbmap[i].b]
    end

    return cmap
end

#-------------------------------------------------------------------------

# Convert N array of ColorTypes.RGBA{Float64} to Nx3 Float64 array

function RGBA2FloatArray(rgbmap::Array{ColorTypes.RGBA{Float64},1})

    N = length(rgbmap)

    cmap = Array{Float64}(undef,N,3)
    for i = 1:N
        cmap[i,:] = [rgbmap[i].r rgbmap[i].g rgbmap[i].b]
    end

    return cmap
end




#----------------------------------------------------------------------------
# Function for applying a 1D Gaussian filter to a signal. Filtering at
# ends are done using zero padding
#
# Usage: sm = gaussfilt1d(s::Array, sigma::Real)

function gaussfilt1d(s::Array, sigma::Real)

    N = length(s)

    r = ceil(Int, 3*sigma)    # Determine filter size
    fw = 2*r + 1

    # Construct filter
    f = [exp(-x.^2/(2*sigma^2)) for x = -r:r]
    f = f/sum(f)

    sm = zeros(size(s))

    # Filter centre section
    for i = r+1:N-r, k = 1:fw
        sm[i] += f[k] * s[i+k-r-1]
    end

    # Filter start section of array using 0 padding
    for i = 1:r, k = 1:fw
        ind = i+k-r-1
        if ind >= 1 && ind <= N
            sm[i] += f[k] * s[ind]
        end
    end

    # Filter end section of array using 0 padding
    for i = N-r+1:N, k = 1:fw
        ind = i+k-r-1
        if ind >= 1 && ind <= N
            sm[i] += f[k] * s[ind]
        end
    end

    return sm
end


#----------------------------------------------------------------------------
"""
Simple 1D linear interpolation of an array of data
```
 Usage:  yi = interp1(x, y, xi)
Arguments:  x - Array of coordinates at which y is defined.
            y - Array of values at coordinates x.
           xi - Coordinate locations at which you wish to interpolate y values.
Returns:   yi - Values linearly interpolated from y at xi.
```
Interpolates y, defined at values x, at locations xi and returns the
corresponding values as yi
x is assumed increasing but not necessarily equi-spaced.
xi values do not need to be sorted.
If any xi are outside the range of x then the corresponding value of
yi is set to the appropriate end value of y.
"""
function interp1(x, y, xi)

    N = length(xi)
    yi = zeros(size(xi))

    minx = minimum(x)
    maxx = maximum(x)

    for i = 1:N
        # Find interval in x that each xi lies within and interpolate its value

        if xi[i] <= minx
            yi[i] = y[1]

        elseif xi[i] >= maxx
            yi[i] = y[end]

        else
            left = maximum(findall(x .<= xi[i]))
            right = minimum(findall(x .> xi[i]))

            yi[i] = y[left] +  (xi[i]-x[left])/(x[right]-x[left]) * (y[right] - y[left])
        end
    end

    return yi
end


#### quick test


cs = ColorScheme([colorant"yellow", colorant"red"])

origcolors = get(cs, 0:0.01:1) # linear interpolation, not perceptually uniform

# generate corrected colormap:

rgbdata = equalisecolourmap("RGB", RGBA{Float64}.(origcolors), "CIEDE2000", [1,0,0])

newcolors = [RGB(rgb...) for rgb in eachrow(rgbdata)]



using Test

# equalisecolourmap: Generate a colour map in Lab space with an uneven
# ramp in lightness and check that this is corrected
rgblab = "LAB"

labmap = zeros(256,3)

labmap[1:127,1] = range(0, stop=20, length=127)

labmap[128:256,1] = range(20, stop=100, length=129)

formula = "CIE76"

W = [1,0,0]

sigma = 1

rgbmap = equalisecolourmap(rgblab, labmap, formula, W, sigma)

# Convert to Nx3 array and then back to lab space. Then check that dL
# is roughly constant

labmap2 = srgb2lab(rgbmap)

#dL = gradient(labmap2[:,1])  # gradient dramas

dL = labmap2[2:end,1] - labmap2[1:end-1,1]

@test maximum(dL[2:end-1]) - minimum(dL[2:end-1]) < 1e-1
