Backports for 0.12.11

docs/src/colordifferences.md

@@ -2,7 +2,7 @@
 
 The [`colordiff`](@ref) function gives an approximate value for the difference between two colors.
 
 ```jldoctest example; setup = :(using Colors)
 julia> colordiff(colorant"red", colorant"darkred")
 23.754149863643036
 
@@ -10,7 +10,7 @@
 52.88136782250768
 
 julia> colordiff(HSV(0, 0.75, 0.5), HSL(0, 0.75, 0.5))
-19.48591066257135
+19.48591066257134
 ```
 
 ```julia

docs/src/colormapsandcolorscales.md

@@ -27,15 +27,15 @@ julia> range(c1, stop=c2, length=15)
  RGB{N0f8}(1.0,0.0,0.0)
  RGB{N0f8}(0.929,0.035,0.0)
  RGB{N0f8}(0.859,0.071,0.0)
- RGB{N0f8}(0.784,0.11,0.0)
+ RGB{N0f8}(0.784,0.106,0.0)
  RGB{N0f8}(0.714,0.145,0.0)
  RGB{N0f8}(0.643,0.18,0.0)
  RGB{N0f8}(0.573,0.216,0.0)
  RGB{N0f8}(0.502,0.251,0.0)
  RGB{N0f8}(0.427,0.286,0.0)
  RGB{N0f8}(0.357,0.322,0.0)
  RGB{N0f8}(0.286,0.357,0.0)
- RGB{N0f8}(0.216,0.392,0.0)
+ RGB{N0f8}(0.216,0.396,0.0)
  RGB{N0f8}(0.141,0.431,0.0)
  RGB{N0f8}(0.071,0.467,0.0)
  RGB{N0f8}(0.0,0.502,0.0)

src/conversions.jl

@@ -87,9 +87,12 @@ correct_gamut(c::CV) where {T<:Union{N0f8,N0f16,N0f32,N0f64},
 correct_gamut(c::CV) where {CV<:TransparentRGB} =
     CV(clamp01(red(c)), clamp01(green(c)), clamp01(blue(c)), clamp01(alpha(c))) # for `hex`
 
-function srgb_compand(v::Fractional)
+@inline function srgb_compand(v)
+    F = typeof(0.5 * v)
+    vf = F(v)
+    vc = @fastmath max(vf, F(0.0031308))
     # `pow5_12` is an optimized function to get `v^(1/2.4)`
-    v <= 0.0031308 ? 12.92v : 1.055 * pow5_12(v) - 0.055
+    vf > F(0.0031308) ? muladd(F(1.055), F(pow5_12(vc)), F(-0.055)) : F(12.92) * vf
 end
 
 function _hsx_to_rgb(im::UInt8, v, n, m)
@@ -278,9 +281,11 @@ cnvt(::Type{HSI{T}}, c::Color) where {T} = cnvt(HSI{T}, convert(RGB{T}, c))
 # Everything to XYZ
 # -----------------
 
-function invert_srgb_compand(v::Fractional)
+@inline function invert_srgb_compand(v)
+    F = typeof(0.5 * v)
+    vf = F(v)
     # `pow12_5` is an optimized function to get `x^2.4`
-    v <= 0.04045 ? 1/12.92 * v : pow12_5(1/1.055 * (v + 0.055))
+    vf > F(0.04045) ? pow12_5(muladd(F(1000/1055), vf, F(55/1055))) : F(100/1292) * vf
 end
 
 # lookup table for `N0f8` (the extra two elements are for `Float32` splines)
@@ -306,7 +311,9 @@ function invert_srgb_compand(v::Float32)
 end
 
 function cnvt(::Type{XYZ{T}}, c::AbstractRGB) where T
-    r, g, b = invert_srgb_compand(red(c)), invert_srgb_compand(green(c)), invert_srgb_compand(blue(c))
+    r = invert_srgb_compand(red(c))
+    g = invert_srgb_compand(green(c))
+    b = invert_srgb_compand(blue(c))
     XYZ{T}(0.4124564*r + 0.3575761*g + 0.1804375*b,
            0.2126729*r + 0.7151522*g + 0.0721750*b,
            0.0193339*r + 0.1191920*g + 0.9503041*b)
@@ -320,20 +327,25 @@ function cnvt(::Type{XYZ{T}}, c::xyY) where T
 end
 
 
-const xyz_epsilon = 216 / 24389
-const xyz_kappa   = 24389 / 27
+const xyz_epsilon   = 216 / 24389 # (6/29)^3
+const xyz_kappa     = 24389 / 27  # (29/6)^3*8
+const xyz_kappa_inv = 27 / 24389
 
 function cnvt(::Type{XYZ{T}}, c::Lab, wp::XYZ = WP_DEFAULT) where T
     fy = (c.l + 16) / 116
-    fx = c.a / 500 + fy
+    fx = fy + c.a / 500
     fz = fy - c.b / 200
 
     fx3 = fx^3
+    fy3 = fy^3
     fz3 = fz^3
 
-    x = fx3 > xyz_epsilon ? fx3 : (116fx - 16) / xyz_kappa
-    y = c.l > xyz_kappa * xyz_epsilon ? ((c. l+ 16) / 116)^3 : c.l / xyz_kappa
-    z = fz3 > xyz_epsilon ? fz3 : (116fz - 16) / xyz_kappa
+    epsilon = oftype(fx3, xyz_epsilon)
+    kappa_inv = oftype(fx3, xyz_kappa_inv)
+
+    x = fx3 > epsilon ? fx3 : muladd(116, fx, -16) * kappa_inv
+    y = fy3 > epsilon ? fy3 : c.l * kappa_inv
+    z = fz3 > epsilon ? fz3 : muladd(116, fz, -16) * kappa_inv
 
     XYZ{T}(x*wp.x, y*wp.y, z*wp.z)
 end
@@ -349,15 +361,18 @@ end
 
 
 function cnvt(::Type{XYZ{T}}, c::Luv, wp::XYZ = WP_DEFAULT) where T
-    (u_wp, v_wp) = xyz_to_uv(wp)
+    c.l == zero(c.l) && return XYZ{T}(zero(T), zero(T), zero(T))
 
-    a = (52 * (c.l==0 ? zero(T) : c.l / (c.u + 13 * c.l * u_wp)) - 1) / 3
-    y = c.l > xyz_kappa * xyz_epsilon ? wp.y * ((c.l + 16) / 116)^3 : wp.y * c.l / xyz_kappa
-    b = -5y
-    d = y * (39 * (c.l==0 ? zero(T) : c.l / (c.v + 13 * c.l * v_wp)) - 5)
-    x = d==b ? zero(T) : (d - b) / (a + 1/3)
-    z = a * x + b + zero(T)
+    u_wp, v_wp = xyz_to_uv(wp)
 
+    ls = c.l * oftype(u_wp, xyz_kappa_inv)
+    y = c.l > 8 ? wp.y * ((c.l + 16) / 116)^3 : wp.y * ls
+
+    u = c.u / (13c.l) + u_wp
+    v = c.v / (13c.l) + v_wp
+    v4 = 4v
+    x = y * 9u / v4
+    z = y * (12 - 3u - 20v) / v4
     XYZ{T}(x, y, z)
 end
 
@@ -411,12 +426,16 @@ cnvt(::Type{xyY{T}}, c::Color) where {T} = cnvt(xyY{T}, convert(XYZ{T}, c))
 # Everything to Lab
 # -----------------
 
-function fxyz2lab(v)
-    v > xyz_epsilon ? cbrt(v) : (xyz_kappa * v + 16) / 116
+@inline function fxyz2lab(v)
+    ka = oftype(v, 841 / 108) # (29/6)^2 / 3 = xyz_kappa / 116
+    kb = oftype(v, 16 / 116) # 4/29
+    vc = @fastmath max(v, oftype(v, xyz_epsilon))
+    @fastmath min(cbrt01(vc), muladd(ka, v, kb))
 end
+
 function cnvt(::Type{Lab{T}}, c::XYZ, wp::XYZ = WP_DEFAULT) where T
-    fx, fy, fz = fxyz2lab(c.x / wp.x), fxyz2lab(c.y / wp.y), fxyz2lab(c.z / wp.z)
-    Lab{T}(116fy - 16, 500(fx - fy), 200(fy - fz))
+    f = mapc(fxyz2lab, mapc((x, y) -> x / y, c, wp))
+    Lab{T}(muladd(116, f.y, -16), 500(f.x - f.y), 200(f.y - f.z))
 end
 
 
@@ -501,10 +520,12 @@ function cnvt(::Type{Luv{T}}, c::XYZ, wp::XYZ = WP_DEFAULT) where T
     (u_, v_) = xyz_to_uv(c)
 
     y = c.y / wp.y
+    epsilon = oftype(y, xyz_epsilon)
+    kappa = oftype(y, xyz_kappa)
 
-    l = y > xyz_epsilon ? 116 * cbrt(y) - 16 : xyz_kappa * y
-    u = 13 * l * (u_ - u_wp) + zero(T)
-    v = 13 * l * (v_ - v_wp) + zero(T)
+    l = y > epsilon ? 116 * cbrt(y) - 16 : kappa * y
+    u = 13 * l * (u_ - u_wp)
+    v = 13 * l * (v_ - v_wp)
 
     Luv{T}(l, u, v)
 end
@@ -566,8 +587,8 @@ function cnvt(::Type{DIN99{T}}, c::Lab) where T
     cc = log(1+0.045*g)/(0.045*kch*ke)
 
     # DIN99 chromaticities
-    a99 = g > 0 ? cc * ee / g : zero(T)
-    b99 = g > 0 ? cc * f / g : zero(T)
+    a99 = g > 0 ? convert(T, cc / g * ee) : zero(T)
+    b99 = g > 0 ? convert(T, cc / g * f ) : zero(T)
 
     DIN99{T}(l99, a99, b99)
 

src/utilities.jl

@@ -14,17 +14,109 @@
 # mod6 supports the input `x` in [-2^28, 2^29]
 mod6(::Type{T}, x::Int32) where T = unsafe_trunc(T, x - 6 * ((widemul(x, 0x2aaaaaaa) + Int64(0x20000000)) >> 0x20))
 
-# TODO: move `pow7` from "src/differences.jl" to here
+# Approximation of the reciprocal of the cube root, x^(-1/3).
+# assuming that x > 0.003, the conditional branches are omitted.
+@inline function rcbrt(x::Float64)
+    ix = reinterpret(UInt64, x)
+    e0 = (ix >> 0x34) % UInt32
+    ed = e0 ÷ 0x3
+    er = e0 - ed * 0x3
+    a = 0x000b_f2d7 - 0x0005_5718 * er
+    e = (UInt32(1363) - ed) << 0x14 | a
+    t1 = reinterpret(Float64, UInt64(e) << 0x20)
+    h1 = muladd(t1^2, -x * t1, 1.0)
+    t2 = muladd(@evalpoly(h1, 1/3, 2/9, 14/81), h1 * t1, t1)
+    h2 = muladd(t2^2, -x * t2, 1.0)
+    t3 = muladd(muladd(2/9, h2, 1/3), h2 * t2, t2)
+    reinterpret(Float64, reinterpret(UInt64, t3) & 0xffff_ffff_8000_0000)
+end
+@inline function rcbrt(x::Float32)
+    ix = reinterpret(UInt32, x)
+    e0 = ix >> 0x17 + 0x2
+    ed = e0 ÷ 0x3
+    er = e0 - ed * 0x3
+    a = 0x005f_9cbe - 0x002a_bd7d * er
+    t1 = reinterpret(Float32, (UInt32(169) - ed) << 0x17 | a)
+    h1 = muladd(t1^2, -x * t1, 1.0f0)
+    t2 = muladd(muladd(2/9f0, h1, 1/3f0), h1 * t1, t1)
+    h2 = muladd(t2^2, -x * t2, 1.0f0)
+    t3 = muladd(1/3f0, h2 * t2, t2)
+    reinterpret(Float32, reinterpret(UInt32, t3) & 0xffff_f000)
+end
+
+cbrt01(x) = cbrt(x)
+@inline function cbrt01(x::Float64)
+    r = rcbrt(x) # x^(-1/3)
+    h = muladd(r^2, -x * r, 1.0)
+    e = muladd(2/9, h, 1/3) * h * r
+    muladd(r, x * r, x * e * (r + r + e)) # x * x^(-2/3)
+end
+@inline function cbrt01(x::Float32)
+    r = Float64(rcbrt(x)) # x^(-1/3)
+    h = muladd(r^2, -Float64(x) * r, 1.0)
+    e = muladd(muladd(14/81, h, 2/9), h, 1/3) * h
+    Float32(1 / muladd(r, e, r))
+end
 
 pow3_4(x) = (y = @fastmath(sqrt(x)); y*@fastmath(sqrt(y))) # x^(3/4)
 
 # `pow5_12` is called from `srgb_compand`.
-# `cbrt` generates a function call, so there is little benefit of `@fastmath`.
 pow5_12(x) = pow3_4(x) / cbrt(x) # 5/12 == 1/2 + 1/4 - 1/3 == 3/4 - 1/3
+@inline function pow5_12(x::Float64)
+    @noinline _cbrt(x) = cbrt01(x)
+    p3_4 = pow3_4(x)
+    # x^(-1/6)
+    if x < 0.02
+        t0 = @evalpoly(x, 3.1366722556806232,
+            -221.51395962221136, 19788.889459114234, -905934.6541469148, 1.5928561711645417e7)
+    elseif x < 0.12
+        t0 = @evalpoly(x, 2.3135905865468644,
+            -26.43664640894651, 385.0146581045545, -2890.0920682466267, 8366.343115590817)
+    elseif x < 1.2
+        t0 = @evalpoly(x, 1.7047813285940905, -3.1261253501167308,
+            7.498744828350077, -10.100319516746419, 6.820601476522508, -1.7978894213531524)
+    else
+        return p3_4 / _cbrt(x)
+    end
+    # x^(-1/3)
+    t1 = t0 * t0
+    h1 = muladd(t1^2, -x * t1, 1.0)
+    t2 = muladd(h1, 1/3 * t1, t1)
+    h2 = muladd(t2^2, -x * t2, 1.0)
+    t2h = @evalpoly(h2, 1/3, 2/9, 14/81) * h2 * t2 # Taylor series of (1-h)^(-1/3)
+    # x^(3/4) * x^(-1/3)
+    muladd(p3_4, t2, p3_4 * t2h)
+end
+@inline function pow5_12(x::Float32)
+    # x^(-1/3)
+    rc = rcbrt(x)
+    rcx = -rc * x
+    rch = muladd(muladd(rc, x, rcx), -rc^2, muladd(rc^2, rcx, 1.0f0)) # 1 - x * rc^3
+    rce = muladd(2/9f0, rch, 1/3f0) * rch * rc
+    # x^(3/4)
+    p3_4_f64 = pow3_4(Float64(x))
+    p3_4r = reinterpret(Float64, reinterpret(UInt64, p3_4_f64) & 0xffffffff_e0000000)
+    p3_4 = Float32(p3_4r)
+    p3_4e = Float32(p3_4_f64 - p3_4r)
+    # x^(3/4) * x^(-1/3)
+    muladd(p3_4, rc, muladd(p3_4, rce, p3_4e * rc))
+end
 
 # `pow12_5` is called from `invert_srgb_compand`.
-# 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
+pow12_5(x) = pow12_5(Float64(x))
+pow12_5(x::BigFloat) = x^big"2.4"
+@inline function pow12_5(x::Float64)
+    # x^0.4
+    t1 = @evalpoly(@fastmath(min(x, 1.75)), 0.24295462640373672,
+        1.7489099720303518, -1.9919942887850166, 1.3197188815160004, -0.3257258790067756)
+    t2 = muladd(2/5, muladd(x / t1^2, @fastmath(sqrt(t1)), -t1), t1) # Newton's method
+    t3 = muladd(2/5, muladd(x / t2^2, @fastmath(sqrt(t2)), -t2), t2)
+    t4 = muladd(2/5, muladd(x / t3^2, @fastmath(sqrt(t3)), -t3), t3)
+    # x^0.4 * x^2
+    rx = reinterpret(Float64, reinterpret(UInt64, x) & 0xffffffff_f8000000) # hi
+    e = x - rx # lo
+    muladd(t4, rx^2, t4 * (rx + rx + e) * e)
+end
 
 
 # Linear interpolation in [a, b] where x is in [0,1],
@@ -208,9 +300,8 @@
 end
 function _weighted_color_mean(w1::Real, c1::Colorant{T1}, c2::Colorant{T2}) where {T1,T2}
     @fastmath min(w1, oneunit(w1) - w1) >= zero(w1) || throw(DomainError(w1, "`w1` must be in [0, 1]"))
-    weight1 = convert(promote_type(T1, T2), w1) # TODO: Consider the need for this
-    weight2 = oneunit(weight1) - weight1
-    mapc((x, y) -> convert(promote_type(T1, T2), muladd(weight1, x, weight2 * y)), c1, c2)
+    w2 = oneunit(w1) - w1
+    mapc((x, y) -> convert(promote_type(T1, T2), muladd(w1, x, w2 * y)), c1, c2)
 end
 function _weighted_color_mean(w1::Integer, c1::C, c2::C) where C <: Colorant
     (w1 & 0b1) === w1 || throw(DomainError(w1, "`w1` must be in [0, 1]"))

test/colordiff.jl

@@ -33,8 +33,8 @@ using Colors, Test
         ((50.0000,   2.5000,   0.0000), (50.0000,   3.2972,   0.0000),  1.0000),
         ((50.0000,   2.5000,   0.0000), (50.0000,   1.8634,   0.5757),  1.0000),
         ((50.0000,   2.5000,   0.0000), (50.0000,   3.2592,   0.3350),  1.0000),
-        ((60.2574, -34.0099,  36.2577), (60.4626, -34.1751,  39.4387),  1.2644),
-        ((63.0109, -31.0961,  -5.8663), (62.8187, -29.7946,  -4.0864),  1.2530),
+        ((60.2574, -34.0099,  36.2677), (60.4626, -34.1751,  39.4387),  1.2644),
+        ((63.0109, -31.0961,  -5.8663), (62.8187, -29.7946,  -4.0864),  1.2630),
         ((61.2901,   3.7196,  -5.3901), (61.4292,   2.2480,  -4.9620),  1.8731),
         ((35.0831, -44.1164,   3.7933), (35.0232, -40.0716,   1.5901),  1.8645),
         ((22.7233,  20.0904, -46.6940), (23.0331,  14.9730, -42.5619),  2.0373),
@@ -45,13 +45,14 @@ using Colors, Test
         ((2.0776,    0.0795,  -1.1350), ( 0.9033,  -0.0636,  -0.5514),  0.9082)]
 
 
-    eps_cdiff = 0.01
+    eps_cdiff = 0.00005
 
     @testset "CIEDE2000" begin
         metric = DE_2000()
         for (a, b, dexpect) in abds
-            @test abs(dexpect - colordiff(Lab(a...), Lab(b...); metric=metric)) < eps_cdiff
-            @test abs(dexpect - colordiff(Lab(b...), Lab(a...); metric=metric)) < eps_cdiff
+            a64, b64 = Lab(a...), Lab(b...)
+            @test colordiff(a64, b64; metric=metric) ≈ dexpect atol=eps_cdiff
+            @test colordiff(b64, a64; metric=metric) ≈ dexpect atol=eps_cdiff
         end
     end