Fix erroneous MSC(h) and improve accuracy of MSC(h, l) (JuliaGrap…

…hics#349)

src/algorithms.jl

@@ -171,35 +171,24 @@ function MSC(h)
     #Wrap h to [0, 360] range]
     h = mod(h, 360)
 
-    p=0 #variable
-    o=0 #min
-    t=0 #max
-
+    # p=0 #variable
+    # o=0 #min
+    # t=0 #max
     #Selecting edge of RGB cube; R=1 G=2 B=3
-    if h0 <= h < h1
-        p=2; o=3; t=1
-    elseif h1 <= h < h2
-        p=1; o=3; t=2
-    elseif h2 <= h < h3
-        p=3; o=1; t=2
-    elseif h3 <= h < h4
-        p=2; o=1; t=3
-    elseif h4 <= h < h5
-        p=1; o=2; t=3
-    elseif h5 <= h || h < h0
-        p=3; o=2; t=1
-    end
-
-    col=zeros(3)
-
-    #check if we are directly on the edge of the RGB cube (within some tolerance)
-    for edge in [h0, h1, h2, h3, h4, h5]
-        if edge - 200eps() < h < edge + 200eps()
-            col[p] = edge in [h0, h2, h4] ? 0.0 : 1.0
-            col[o] = 0.0
-            col[t] = 1.0
-            return convert(LCHuv, RGB(col[1],col[2],col[3]))
-        end
+    if h < h0
+        p=3; t=1
+    elseif h < h1
+        p=2; t=1
+    elseif h < h2
+        p=1; t=2
+    elseif h < h3
+        p=3; t=2
+    elseif h < h4
+        p=2; t=3
+    elseif h < h5
+        p=1; t=3
+    else
+        p=3; t=1
     end
 
     alpha=-sind(h)
@@ -214,7 +203,6 @@ function MSC(h)
     M = [0.4124564  0.3575761  0.1804375;
          0.2126729  0.7151522  0.0721750;
          0.0193339  0.1191920  0.9503041]'
-    g = 2.4
 
     m_tx=M[t,1]
     m_ty=M[t,2]
@@ -223,20 +211,17 @@ function MSC(h)
     m_py=M[p,2]
     m_pz=M[p,3]
 
-    f1=(4alpha*m_px+9beta*m_py)
-    a1=(4alpha*m_tx+9beta*m_ty)
-    f2=(m_px+15m_py+3m_pz)
-    a2=(m_tx+15m_ty+3m_tz)
+    f1 = 4alpha*m_px+9beta*m_py
+    a1 = 4alpha*m_tx+9beta*m_ty
+    f2 = m_px+15m_py+3m_pz
+    a2 = m_tx+15m_ty+3m_tz
 
     cp=((alpha*un+beta*vn)*a2-a1)/(f1-(alpha*un+beta*vn)*f2)
 
-    #gamma inversion
-    cp = cp <= 0.003 ? 12.92cp : 1.055cp^(1.0/g)-0.05
-    #cp = 1.055cp^(1.0/g)-0.05
-
-    col[p]=clamp01(cp)
-    col[o]=0.0
-    col[t]=1.0
+    col = zeros(3)
+    col[p] = clamp01(srgb_compand(cp))
+    # col[o] = 0.0
+    col[t] = 1.0
 
     return convert(LCHuv, RGB(col[1],col[2],col[3]))
 end
@@ -247,15 +232,32 @@ end
 
 # Maximally saturated color for a specific hue and lightness
 # is found by looking for the edge of LCHuv space.
-function MSC(h,l)
-    pmid=MSC(h)
-
-    if l <= pmid.l
-        pend=LCHuv(0,0,0)
-    elseif l > pmid.l
-        pend=LCHuv(100,0,0)
+function MSC(h, l; linear::Bool=false)
+    if linear
+        pmid = MSC(h)
+        pend_l = l > pmid.l ? 100 : 0
+        return (pend_l-l)/(pend_l-pmid.l) * pmid.c
     end
+    return find_maximum_chroma(LCHuv(l, 0, h))
+end
 
-    a=(pend.l-l)/(pend.l-pmid.l)
-    a*(pmid.c-pend.c)+pend.c
+function find_maximum_chroma(c::T,
+                             low::Real=0,
+                             high::Real=180) where {T<:Union{LCHab, LCHuv}}
+    err = 1e-6
+    high-low < err && return low
+
+    mid = (low + high) / 2
+    lchm = T(c.l, mid, c.h)
+    rgbm = convert(RGB, lchm)
+    clamped = max(red(rgbm), green(rgbm), blue(rgbm)) > 1-err ||
+              min(red(rgbm), green(rgbm), blue(rgbm)) < err
+    if clamped
+        return find_maximum_chroma(c, low, mid)
+    else
+        return find_maximum_chroma(c, mid, high)
+    end
 end
+
+find_maximum_chroma(c::Lab) = find_maximum_chroma(convert(LCHab, c))
+find_maximum_chroma(c::Luv) = find_maximum_chroma(convert(LCHuv, c))

test/algorithms.jl

@@ -1,21 +1,29 @@
 using Test, Colors
 
 @testset "Algorithms" begin
-    @test isconcretetype(eltype(colormap("Grays")))
 
-    @test_throws ArgumentError colormap("Grays", N=10)
+    # issue 346
+    msc_h_diff = 0
+    for hsv_h in 0:0.1:360
+        hsv = HSV(hsv_h, 1.0, 1.0) # most saturated
+        lch = convert(LCHuv, hsv)
+        msc = MSC(lch.h)
+        msc_h_diff = max(msc_h_diff, colordiff(msc, lch))
+    end
+    @test msc_h_diff < 1
+
+    @test MSC(0, 100) ≈ 0
+    @test MSC(0, 0) ≈ 0
 
-    col = distinguishable_colors(10)
-    @test isconcretetype(eltype(col))
-    local mindiff
-    mindiff = Inf
-    for i = 1:10
-        for j = i+1:10
-            mindiff = min(mindiff, colordiff(col[i], col[j]))
-        end
+    msc_h_l_sat = 1
+    for h = 0:0.1:359.9, l = 1:1:99
+        c = MSC(h, l)
+        hsv = convert(HSV, LCHuv(l, c, h))
+        msc_h_l_sat = min(msc_h_l_sat, max(hsv.s, hsv.v))
     end
-    @test mindiff > 8
+    @test msc_h_l_sat > 1 - 1e-4
 
-    cols = distinguishable_colors(1)
-    @test colordiff(distinguishable_colors(1, cols; dropseed=true)[1], cols[1]) > 50
+    # issue 346
+    @test MSC(0, 90, linear=true) > MSC(0, 90)
+    @test MSC(280, 50, linear=true) < MSC(280, 50)
 end

test/colormaps.jl

@@ -1,5 +1,45 @@
 @testset "Colormaps" begin
     @test length(colormap("RdBu", 100)) == 100
 
+    @test isconcretetype(eltype(colormap("Grays")))
+
+    @test_throws ArgumentError colormap("Grays", N=10)
+
+    # The outputs of `colormap()` were slightly affected by the bug fix of
+    # `MSC(h)` (issue #349).
+    # The following were generated by `colormap()` in Colors.jl v0.9.6.
+    blues_old   = ["F4FDFF", "B3E3F4", "65B9E7", "2978BE", "0B2857"]
+    greens_old  = ["FAFFF7", "B6EEA0", "75C769", "308B40", "00391A"]
+    grays_old   = ["FFFFFF", "DCDCDC", "A9A9A9", "626262", "000000"]
+    oranges_old = ["FFFBF6", "FFD6B4", "FF9D5F", "E2500D", "732108"]
+    purples_old = ["FBFBFB", "DBD7F6", "AFA7E6", "7666BF", "3C0468"]
+    reds_old    = ["FFF1EE", "FFC4B9", "FF8576", "E72823", "6D0B0C"]
+    rdbu_old    = ["610102", "FF8D7B", "F9F8F9", "76B4E8", "092C58"]
+
+    to_rgb(s) = parse(RGB, "#"*s)
+    max_colordiff(a1, a2) = reduce(max, colordiff.(a1, a2))
+    @test max_colordiff(colormap("Blues", 5), to_rgb.(blues_old)) < 1
+    @test max_colordiff(colormap("Greens", 5), to_rgb.(greens_old)) < 1
+    @test max_colordiff(colormap("Grays", 5), to_rgb.(grays_old)) < 1
+    @test max_colordiff(colormap("Oranges", 5), to_rgb.(oranges_old)) < 1
+    @test max_colordiff(colormap("Purples", 5), to_rgb.(purples_old)) < 1
+    @test max_colordiff(colormap("Reds", 5), to_rgb.(reds_old)) < 1
+    @test max_colordiff(colormap("RdBu", 5), to_rgb.(rdbu_old)) < 1
+
     # TODO: add more tests
+
+    col = distinguishable_colors(10)
+    @test isconcretetype(eltype(col))
+    local mindiff
+    mindiff = Inf
+    for i = 1:10
+        for j = i+1:10
+            mindiff = min(mindiff, colordiff(col[i], col[j]))
+        end
+    end
+    @test mindiff > 8
+
+    cols = distinguishable_colors(1)
+    @test colordiff(distinguishable_colors(1, cols; dropseed=true)[1], cols[1]) > 50
+
 end