From 50b3acbd8d997245e62c51b7467f37a561ab2ccc Mon Sep 17 00:00:00 2001 From: cormullion Date: Thu, 18 Aug 2022 16:52:31 +0100 Subject: [PATCH] bump versions for 1.3 --- CHANGELOG.md | 2 +- Project.toml | 2 +- src/equalizecolormap,jl | 758 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 760 insertions(+), 2 deletions(-) create mode 100644 src/equalizecolormap,jl diff --git a/CHANGELOG.md b/CHANGELOG.md index e436f94..26b09dc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,6 @@ # Changelog -## [v1.3.0] - forthcoming +## [v1.3.0] -2022-08-18 ### Added diff --git a/Project.toml b/Project.toml index d2c37eb..ed4baf0 100644 --- a/Project.toml +++ b/Project.toml @@ -20,7 +20,7 @@ Colors = "0.9, 0.10, 0.11, 0.12, 0.13" FileIO = "^1" ImageMagick = "0.7, 1" Images = "0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25" -Interpolations = "0.12, 0.13" +Interpolations = "0.12, 0.13, 0.14" QuartzImageIO = "^0.7" julia = "1" diff --git a/src/equalizecolormap,jl b/src/equalizecolormap,jl new file mode 100644 index 0000000..0767ac7 --- /dev/null +++ b/src/equalizecolormap,jl @@ -0,0 +1,758 @@ +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