Clean up "src/differences.jl" (#476)

This exports the re-implemented version of `mean_hue`.

docs/src/colordifferences.md

@@ -4,7 +4,7 @@ The [`colordiff`](@ref) function gives an approximate value for the difference b
 
 ```jldoctest example; setup = :(using Colors)
 julia> colordiff(colorant"red", colorant"darkred")
-23.754149863643036
+23.75414986364304
 
 julia> colordiff(colorant"red", colorant"blue")
 52.88136782250768

docs/src/constructionandconversion.md

@@ -146,4 +146,5 @@ Depending on the source and destination colorspace, this may not be perfectly lo
 parse
 hex
 normalize_hue
+mean_hue
 ```

src/Colors.jl

@@ -18,7 +18,7 @@ export weighted_color_mean,
        distinguishable_colors, whitebalance,
        colordiff, DE_2000, DE_94, DE_JPC79, DE_CMC, DE_BFD, DE_AB, DE_DIN99, DE_DIN99d, DE_DIN99o,
        MSC, sequential_palette, diverging_palette, colormap,
-       normalize_hue,
+       normalize_hue, mean_hue,
        colormatch, CIE1931_CMF, CIE1964_CMF, CIE1931J_CMF, CIE1931JV_CMF, CIE2006_2_CMF, CIE2006_10_CMF
 
 # Early utilities

src/differences.jl

@@ -151,22 +151,6 @@ Construct a metric using Euclidean color difference equation applied in the
 """
 DE_DIN99o()
 
-
-# Compute the mean of two hue angles
-function mean_hue(h1, h2)
-    if abs(h2 - h1) > 180
-        if h1 + h2 < 360
-            mh = (h1 + h2 + 360) / 2
-        else
-            mh = (h1 + h2 - 360) / 2
-        end
-    else
-        mh = (h1 + h2) / 2
-    end
-
-    mh
-end
-
 # Color difference metrics
 # ------------------------
 
@@ -181,12 +165,11 @@ end
 #   The CIEDE2000 color difference metric evaluated between a and b.
 #
 
-pow7(x) = (y = x*x*x; y*y*x)
-pow7(x::Integer) = pow7(Float64(x))
-const twentyfive7 = pow7(25)
 
 # Delta E 2000
 function _colordiff(ai::Color, bi::Color, m::DE_2000)
+    twentyfive7 = pow7(25)
+
     # Ensure that the input values are in L*a*b* space
     a_Lab = convert(Lab, ai)
     b_Lab = convert(Lab, bi)
@@ -195,31 +178,18 @@ function _colordiff(ai::Color, bi::Color, m::DE_2000)
     mc = (chroma(a_Lab) + chroma(b_Lab))/2
     g = (1 - sqrt(pow7(mc) / (pow7(mc) + twentyfive7))) / 2
 
+    a_Lab_r = Lab(a_Lab.l, a_Lab.a * (1 + g), a_Lab.b)
+    b_Lab_r = Lab(b_Lab.l, b_Lab.a * (1 + g), b_Lab.b)
+
     # Convert to L*C*h, where the remainder of the calculations are performed
-    a = convert(LCHab, Lab(a_Lab.l, a_Lab.a * (1 + g), a_Lab.b))
-    b = convert(LCHab, Lab(b_Lab.l, b_Lab.a * (1 + g), b_Lab.b))
+    a = convert(LCHab, a_Lab_r)
+    b = convert(LCHab, b_Lab_r)
 
     # Calculate the delta values for each channel
-    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
-    if a.c * b.c == 0
-        dh = zero(dh)
-    elseif dh > 180
-        dh -= 360
-    elseif dh < -180
-        dh += 360
-    end
-    # Calculate H*
-    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)
-
-    # Calculate mean L* and C* values
-    ml, mc = (a.l + b.l) / 2, (a.c + b.c) / 2
+    dl, dc, dh = b.l - a.l, b.c - a.c, delta_h(b_Lab_r, a_Lab_r)
 
-    # Calculate mean hue value
-    if a.c * b.c == 0
-        mh = a.h + b.h
-    else
-        mh = mean_hue(a.h, b.h)
-    end
+    # Calculate mean L*, C* and hue values
+    ml, mc, mh = (a.l + b.l) / 2, (a.c + b.c) / 2, mean_hue(a, b)
 
     # lightness weight
     mls = (ml - 50)^2
@@ -236,9 +206,9 @@ function _colordiff(ai::Color, bi::Color, m::DE_2000)
     sh = 1 + 0.015 * mc * t
 
     # rotation term
-    dtheta = 30 * exp(-((mh - 275)/25)^2)
+    dtheta = 60 * exp(-((mh - 275)/25)^2)
     cr = 2 * sqrt(pow7(mc) / (pow7(mc) + twentyfive7))
-    tr = -sind(2*dtheta) * cr
+    tr = -sind(dtheta) * cr
 
     # Final calculation
     sqrt((dl/(m.kl*sl))^2 + (dc/(m.kc*sc))^2 + (dh/(m.kh*sh))^2 +
@@ -247,59 +217,31 @@ end
 
 # Delta E94
 function _colordiff(ai::Color, bi::Color, m::DE_94)
-
     a = convert(LCHab, ai)
     b = convert(LCHab, bi)
 
     # Calculate the delta values for each channel
-    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
-    if a.c * b.c == 0
-        dh = zero(dh)
-    elseif dh > 180
-        dh -= 360
-    elseif dh < -180
-        dh += 360
-    end
-    # Calculate H*
-    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)
+    dl, dc, dh = b.l - a.l, b.c - a.c, delta_h(b, a)
 
     # Lightness, hue, chroma correction terms
     sl = 1
     sc = 1 + m.k1 * a.c
     sh = 1 + m.k2 * a.c
 
     sqrt((dl/(m.kl*sl))^2 + (dc/(m.kc*sc))^2 + (dh/(m.kh*sh))^2)
-
 end
 
 # Metric form of jpc79 color difference equation (mostly obsolete)
 function _colordiff(ai::Color, bi::Color, m::DE_JPC79)
-
     # Convert directly into LCh
     a = convert(LCHab, ai)
     b = convert(LCHab, bi)
 
     # Calculate deltas in each direction
-    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
-    if a.c * b.c == 0
-        dh = zero(dh)
-    elseif dh > 180
-        dh -= 360
-    elseif dh < -180
-        dh += 360
-    end
-    # Calculate H* from C*'s and h
-    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)
-
-    #Calculate mean lightness
-    ml = (a.l + b.l)/2
-    mc = (a.c + b.c)/2
-    # Calculate mean hue value
-    if a.c * b.c == 0
-        mh = a.h + b.h
-    else
-        mh = mean_hue(a.h, b.h)
-    end
+    dl, dc, dh = b.l - a.l, b.c - a.c, delta_h(b, a)
+
+    # Calculate mean lightness, chroma and hue
+    ml, mc, mh = (a.l + b.l) / 2, (a.c + b.c) / 2, mean_hue(a, b)
 
     # L* adjustment term
     sl = 0.08195*ml/(1+0.01765*ml)
@@ -318,28 +260,17 @@ function _colordiff(ai::Color, bi::Color, m::DE_JPC79)
 
     # Calculate the final difference
     sqrt((dl/sl)^2 + (dc/sc)^2 + (dh/sh)^2)
-
 end
 
 
 # Metric form of the cmc color difference
 function _colordiff(ai::Color, bi::Color, m::DE_CMC)
-
     # Convert directly into LCh
     a = convert(LCHab, ai)
     b = convert(LCHab, bi)
 
     # Calculate deltas in each direction
-    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
-    if a.c * b.c == 0
-        dh = zero(dh)
-    elseif dh > 180
-        dh -= 360
-    elseif dh < -180
-        dh += 360
-    end
-    # Calculate H* from C*'s and h
-    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)
+    dl, dc, dh = b.l - a.l, b.c - a.c, delta_h(b, a)
 
     # L* adjustment term
     if (a.l <= 16)
@@ -364,7 +295,6 @@ function _colordiff(ai::Color, bi::Color, m::DE_CMC)
 
     # Calculate the final difference
     sqrt((dl/(m.kl*sl))^2 + (dc/(m.kc*sc))^2 + (dh/sh)^2)
-
 end
 
 # The BFD color difference equation
@@ -394,36 +324,27 @@ function _colordiff(ai::Color, bi::Color, m::DE_BFD)
     b = LCHab(lb, chroma(b_Lab), hue(b_Lab))
 
     # Calculate deltas in each direction
-    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
-    if a.c * b.c == 0
-        dh = zero(dh)
-    elseif dh > 180
-        dh -= 360
-    elseif dh < -180
-        dh += 360
-    end
-
-    # Calculate H* from C*'s and h
-    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)
+    dl, dc, dh = b.l - a.l, b.c - a.c, delta_h(b, a)
 
     # Find the mean value of the inputs to use as the "standard"
-    ml, mc = (a.l + b.l)/2, (a.c + b.c)/2
-
-    # Calculate mean hue value
-    if a.c * b.c == 0
-        mh = a.h + b.h
-    else
-        mh = mean_hue(a.h, b.h)
-    end
+    ml, mc, mh = (a.l + b.l) / 2, (a.c + b.c) / 2, mean_hue(a, b)
 
     # Correction terms for a variety of nonlinearities in CIELAB.
     g = sqrt(mc^4/(mc^4 + 14000))
 
-    t = 0.627 + 0.055*cosd(mh - 245) - 0.04*cosd(2*mh - 136) + 0.07*cosd(3*mh - 32) + 0.049*cosd(4*mh + 114) - 0.015*cosd(5*mh + 103)
+    t = 0.627 + 0.055 * cosd( mh - 245) -
+                0.040 * cosd(2mh - 136) +
+                0.070 * cosd(3mh -  32) +
+                0.049 * cosd(4mh + 114) -
+                0.015 * cosd(5mh + 103)
 
     rc = sqrt(mc^6/(mc^6 + 7e7))
 
-    rh = -0.26cosd(mh-308) - 0.379cosd(2*mh-160) - 0.636*cosd(3*mh - 254) + 0.226cosd(4*mh + 140) - 0.194*cosd(5*mh + 280)
+    rh = -0.260 * cosd( mh - 308) -
+          0.379 * cosd(2mh - 160) -
+          0.636 * cosd(3mh - 254) +
+          0.226 * cosd(4mh + 140) -
+          0.194 * cosd(5mh + 280)
 
     dcc = 0.035*mc/(1+0.00365*mc) + 0.521
     dhh = dcc*(g*t+1-g)

src/utilities.jl

@@ -26,6 +26,8 @@ pow5_12(x) = pow3_4(x) / cbrt(x) # 5/12 == 1/2 + 1/4 - 1/3 == 3/4 - 1/3
 # x^y ≈ exp(y*log(x)) ≈ exp2(y*log2(y)); the middle form is faster
 @noinline pow12_5(x) = x^2 * exp(0.4 * log(x)) # 12/5 == 2.4 == 2 + 0.4
 
+pow7(x) = (y = x*x*x; y*y*x)
+pow7(x::Integer) = pow7(Float64(x)) # avoid overflow
 
 # Linear interpolation in [a, b] where x is in [0,1],
 # or coerced to be if not.
@@ -195,6 +197,51 @@ normalize_hue(c::C) where C <: Union{LCHab, LCHuv} = C(c.l, c.c, normalize_hue(c
 normalize_hue(c::C) where {Cb <: Union{LCHab, LCHuv}, C <: Union{AlphaColor{Cb}, ColorAlpha{Cb}}} =
     C(c.l, c.c, normalize_hue(c.h), c.alpha)
 
+"""
+    mean_hue(h1::Real, h2::Real)
+    mean_hue(a::C, b::C) where {C <: Colorant}
+
+Compute the mean of two hue angles in degrees.
+
+If the inputs are HSV-like or Lab-like color objects, this will also return a
+hue, not a color. If one of the colors is achromatic, i.e. has zero saturation
+or chroma, the hue of the other color is returned instead of the mean.
+"""
+function mean_hue(h1::T, h2::T) where {T <: Real}
+    @fastmath hmin, hmax = minmax(h1, h2)
+    d = 180 - normalize_hue(hmin - hmax - 180)
+    F = typeof(zero(T) / 2)
+    mh = muladd(F(0.5), d, hmin)
+    return mh < 0 ? mh + 360 : mh
+end
+
+function mean_hue(a::C, b::C) where {Cb <: Union{Lab, LCHab, Luv, LCHuv},
+                                     C <: Union{Cb, AlphaColor{Cb}, ColorAlpha{Cb}}}
+    h1 = chroma(a) == 0 ? hue(b) : hue(a)
+    h2 = chroma(b) == 0 ? hue(a) : hue(b)
+    mean_hue(h1, h2)
+end
+function mean_hue(a::C, b::C) where {Cb <: Union{HSV, HSL, HSI},
+                                     C <: Union{Cb, AlphaColor{Cb}, ColorAlpha{Cb}}}
+    mean_hue(a.s == 0 ? b.h : a.h, b.s == 0 ? a.h : b.h)
+end
+mean_hue(a, b) = mean_hue(promote(a, b)...)
+
+function delta_h(a::C, b::C) where {Cb <: Union{Lab, Luv},
+                                    C <: Union{Cb, AlphaColor{Cb}, ColorAlpha{Cb}}}
+    da, db, dc = comp2(a) - comp2(b), comp3(a) - comp3(b), chroma(a) - chroma(b)
+    d = comp3(a) * comp2(b) - comp2(a) * comp3(b)
+    dh = @fastmath sqrt(max(muladd(dc, -dc, muladd(da, da, db^2)), 0))
+    return copysign(dh, d)
+end
+function delta_h(a::C, b::C) where {Cb <: Union{LCHab, LCHuv},
+                                    C <: Union{Cb, AlphaColor{Cb}, ColorAlpha{Cb}}}
+    sh = @fastmath (hue(a) - hue(b) + 180) / 360
+    d =  @fastmath sh - floor(sh) - oftype(sh, 0.5)
+    2 * sqrt(chroma(a) * chroma(b)) * sinpi(d)
+end
+delta_h(a, b) = delta_h(promote(a, b)...)
+
 """
     weighted_color_mean(w1, c1, c2)
 

test/utilities.jl

@@ -160,6 +160,29 @@ using InteractiveUtils # for `subtypes`
         @test normalize_hue(LCHuvA{Float64}(30, 40, -0.5, 0.6)) === LCHuvA{Float64}(30, 40, 359.5, 0.6)
     end
 
+    @testset "mean_hue" begin
+        @test @inferred(mean_hue(0.0, 90.0)) === 45.0
+        @test @inferred(mean_hue(90.0f0, 270.0f0)) === 180.0f0
+        @test @inferred(mean_hue(90, 271)) === 0.5
+
+        @test @inferred(mean_hue(HSV(50, 1, 1), HSV(100, 1, 0.5))) === 75.0
+        @test @inferred(mean_hue(HSL(50, 1, 1), HSL(100, 0, 0.5))) === 50.0
+        @test @inferred(mean_hue(Lab(50, 0, 10), Lab(50, 10, -10))) === 22.5f0
+        @test @inferred(mean_hue(ALuv(50, 0, 0), ALuv(50, 0, 10))) === 90.0f0
+
+        @test_throws Exception mean_hue(RGB(0.1, 0.2, 0.3), RGB(0.4, 0.5, 0.6))
+        @test_throws Exception mean_hue(LChab(10, 20, 30), LCHuv(10, 20, 30))
+    end
+
+    @testset "delta_h" begin
+        @test @inferred(Colors.delta_h(LCHab(50, 60, 70), LCHab(90, 80, 70))) === 0.0f0
+        @test @inferred(Colors.delta_h(LCHuv(50, 60, 80), LCHuv(90,  0, 70))) === 0.0f0
+        @test @inferred(Colors.delta_h(LCHab(50, 60, 70), LCHab(90, 80, 190))) ≈ -120.0f0
+        @test @inferred(Colors.delta_h(LCHab(90, 80, 70), LCHab(50, 60, 310.0))) ≈ 120.0
+        @test @inferred(Colors.delta_h(Lab(50, -10, 10), Lab(50, 10, -10))) ≈ 28.284271f0
+        @test @inferred(Colors.delta_h(ALuv(50, 0, 0), ALuv(50, 0, 10))) === 0.0f0
+    end
+
     # test utility function weighted_color_mean
     parametric2 = [GrayA,AGray32,AGray]
     parametric3 = ColorTypes.parametric3