
mapch(f, x, y) = mapc(f, x, y)
mapch(f, x::C, y::C) where C <: Union{HSV, HSL, HSI} =
    C(f(Hue(x), Hue(y)), f(x.s, y.s), f(comp3(x), comp3(y)))
mapch(f, x::C, y::C) where C <: Union{LCHab, LCHuv} =
    C(f(x.l, y.l), f(x.c, y.c), f(Hue(x), Hue(y)))

mapca(fc, a, x::C) where C <: TransparentColorN{2} = C(fc(comp1(x)), a)
mapca(fc, a, x::C) where C <: TransparentColorN{3} = C(fc(comp1(x)), fc(comp2(x)), a)
mapca(fc, a, x::C) where C <: TransparentColorN{4} = C(fc(comp1(x)), fc(comp2(x)), fc(comp3(x)), a)
mapca(fc, a, x::C, y::C) where C <: TransparentColorN{2} =
    C(fc(comp1(x), comp1(y)), a)
mapca(fc, a, x::C, y::C) where C <: TransparentColorN{3} =
    C(fc(comp1(x), comp1(y)), fc(comp2(x), comp2(y)), a)
mapca(fc, a, x::C, y::C) where C <: TransparentColorN{4} =
    C(fc(comp1(x), comp1(y)), fc(comp2(x), comp2(y)), fc(comp3(x), comp3(y)), a)
mapca(fc, a, x::C, y::C) where C <: TransparentColorN{4, <:Union{HSV, HSL, HSI}} =
    C(fc(Hue(x), Hue(y)), fc(x.s, y.s), fc(comp3(x), comp3(y)), a)
mapca(fc, a, x::C, y::C) where C <: TransparentColorN{4, <:Union{LCHab, LCHuv}} =
    C(fc(x.l, y.l), fc(x.c, y.c), fc(Hue(x), Hue(y)))

# complement
_n(v::T) where T = oneunit(T) - v

# linear interpolation
_w(v1::T, v2::T, w) where T = _w(v1, _n(w), v2, w)

function _w(h1::Hue{T}, h2::Hue{T}, w) where T
    h1a = ifelse(h1.isgray, h2.angle, h1.angle)
    h2a = ifelse(h2.isgray, h1.angle, h2.angle)
    d0 = h2a - h1a
    d = ifelse(abs(d0) > T(180), d0 - copysign(T(360), d0), d0)
    a = muladd(d, w, h1a)
    convert(T, ifelse(a > T(360), a - T(360), ifelse(a < T(0), a + T(360), a)))
end

_w(v1::T, w1, v2::T, w2) where T = convert(T, muladd(w1, v1, w2 * v2))

function _w(v1::T, w1::T, v2::T, w2::T) where T <: Union{N0f8, N0f16}
    r1, rw1 = reinterpret(v1), reinterpret(w1)
    r2, rw2 = reinterpret(v2), reinterpret(w2)
    m = widemul(r1, rw1) + widemul(r2, rw2)
    if T === N0f8
        z = (m + ((m + 0x80) >> 0x8) + 0x80) >> 0x8
    else
        z = (m + ((m + 0x8000) >> 0x10) + 0x8000) >> 0x10
    end
    reinterpret(T, unsafe_trunc(typeof(r1), z))
end

function _w(h1::Hue{T}, w1, h2::Hue{T}, w2) where T
    w = w1 + w2
    w == oneunit(w) && return _w(h1, h2, w2)
    w == zero(w) && return zero(T)
    _w(h1, h2, w2 / w) * w
end

_w_safe(v1::T, w1, v2::T, w2) where T = min(oneunit(T), _w(v1, w1, v2, w2))
_w_safe(v1::T, w1, v2::T, w2) where T <: FixedPoint =
    convert(T, min(float(typemax(T)), muladd(w1, float(v1), w2 * float(v2))))

_comp(op::CompositeOperation{:clear}, c1, c2) = mapc(v1 -> zero(v1), c1)

_comp(op::CompositeOperation{:copy}, c1, c2) = c2

_comp(op::CompositeOperation{:destination}, c1, c2) = c1

function _comp(op::typeof(CompositeSourceOver), c1, c2)
    k1 = mul(alpha(c1), _n(alpha(c2)))
    k2 = alpha(c2)
    a = k1 + k2
    k1a = a == zero(a) ? a : k1 / a
    k2a = a == zero(a) ? a : oneunit(a) - k1a
    mapca((v1, v2) -> _w(v1, k1a, v2, k2a), a, c1, c2)
end

_comp(op::typeof(CompositeDestinationOver), c1, c2) = _comp(CompositeSourceOver, c2, c1)

function _comp(op::typeof(CompositeSourceIn), c1, c2)
    a = mul(alpha(c1), alpha(c2))
    ifelse(a == zero(a), mapca(v -> a, a, c2), mapca(v -> v, a, c2))
end

_comp(op::typeof(CompositeDestinationIn), c1, c2) = _comp(CompositeSourceIn, c2, c1)

function _comp(op::typeof(CompositeSourceOut), c1, c2)
    a = mul(_n(alpha(c1)), alpha(c2))
    ifelse(a == zero(a), mapca(v -> a, a, c2), mapca(v -> v, a, c2))
end

_comp(op::typeof(CompositeDestinationOut), c1, c2) = _comp(CompositeSourceOut, c2, c1)

function _comp(op::typeof(CompositeSourceAtop), c1, c2)
    a1, a2 = alpha(c1), alpha(c2)
    k1a = a1 == zero(a1) ? a1 : _n(a2)
    k2a = a1 == zero(a1) ? a1 : a2
    mapca((v1, v2) -> _w(v1, k1a, v2, k2a), a1, c1, c2)
end

_comp(op::typeof(CompositeDestinationAtop), c1, c2) = _comp(CompositeSourceAtop, c2, c1)

function _comp(op::CompositeOperation{:xor}, c1, c2)
    m = mul(alpha(c1), alpha(c2))
    k1, k2 = alpha(c1) - m, alpha(c2) - m
    a = k1 + k2
    k1a = a == zero(a) ? a : k1 / a
    k2a = a == zero(a) ? a : oneunit(a) - k1a
    mapca((v1, v2) -> _w(v1, k1a, v2, k2a), a, c1, c2)
end

function _comp(op::CompositeOperation{:lighter}, c1, c2)
    k1, k2 = float(alpha(c1)), float(alpha(c2))
    a = k1 + k2
    ac = min(oneunit(a), a)
    k1a = a == zero(a) ? a : k1 / ac
    k2a = a == zero(a) ? a : k2 / ac
    mapca((v1, v2) -> _w_safe(v1, k1a, v2, k2a), ac, c1, c2)
end

mix_alpha(opacity::Nothing, a) = a
mix_alpha(opacity::Real, a) = opacity * a

@inline (mode::BlendMode)(c1, c2; opacity=nothing, op=CompositeSourceOver) =
    _blend_c(mode, c1, c2, opacity, op)

@inline (op::CompositeOperation)(c1, c2; opacity=nothing, mode=BlendNormal) =
    _blend_c(mode, c1, c2, opacity, op)

"""
    blend(c1, c2; mode=BlendNormal, opacity=1, op=CompositeSourceOver)

Create the mixed color of two colors `c1` and `c2`. The `c1` means the backdrop
color and the `c2` means the source color.

`mode` specifies the blend mode, e.g. [`BlendMultiply`](@ref).

`opacity` modifies the source (i.e. `c2`-side) alpha by multiplication.

`op` specifies the composite operations, e.g. [`CompositeSourceAtop`](@ref).

The return type is the same as `c1`.

# Examples
```jldoctest
julia> blend(RGB(1, 0.5, 0), RGB(0, 0.5, 1), mode=BlendLighten)
RGB{Float64}(1.0,0.5,1.0)
```
"""
@inline blend(c1, c2; mode::BlendMode=BlendNormal, opacity=nothing, op=CompositeSourceOver) =
    _blend_c(mode, c1, c2, opacity, op)


@inline _blend_c(mode::BlendMode, c1::Color, c2::Color, opacity, op) =
    _blend_cc(mode, c1, convert(typeof(c1), c2), opacity, op)

@inline _blend_c(mode::BlendMode, c1::Color, c2::TransparentColor, opacity, op) =
    _blend_cc(mode, c1, convert(typeof(c1), color(c2)), mix_alpha(opacity, alpha(c2)), op)

@inline _blend_c(mode::BlendMode, c1::TransparentColor, c2::Color, opacity, op) =
    _blend_tc(mode, c1, convert(color_type(c1), c2), opacity, op)

@inline _blend_c(mode::BlendMode, c1::TransparentColor, c2::TransparentColor, opacity, op) =
    _blend_tc(mode, c1, convert(color_type(c1), color(c2)), mix_alpha(opacity, alpha(c2)), op)


# without opacity
@inline _blend_cc(mode::BlendMode, c1, c2, ::Nothing, ::OpaqueInSourceOutOperation) =
    _blend(mode, c1, c2)

_blend_cc(mode::BlendMode, c1, c2, ::Nothing, ::OpaqueInDestinationOutOperation) = c1

@inline _blend_cc(mode::BlendMode, c1, c2, ::Nothing, op::CompositeOperation{:lighter}) =
    mapc((v1, v2) -> _w_safe(v1, oneunit(v1), v2, oneunit(v2)), c1, _blend(mode, c1, c2))

@inline function _blend_tc(mode::BlendMode, c1, c2, ::Nothing, op)
    cm = mapch((v1, v2) -> _w(v1, v2, alpha(c1)), c2, _blend(mode, color(c1), c2))
    _comp(op, c1, typeof(c1)(cm))
end


# with opacity
@inline _blend_cc(mode::BlendMode, c1, c2, opacity::Real, ::typeof(CompositeSourceOver)) =
    mapch((v1, v2) -> _w(v1, v2, opacity), c1, _blend(mode, c1, c2))

@inline _blend_cc(mode::BlendMode, c1, c2, opacity::Real, ::typeof(CompositeSourceAtop)) =
    _blend_cc(mode, c1, c2, opacity, CompositeSourceOver)

_blend_cc(mode::BlendMode, c1, c2, opacity::Real, ::typeof(CompositeDestination)) = c1

_blend_cc(mode::BlendMode, c1, c2, opacity::Real, ::typeof(CompositeDestinationOver)) = c1

@inline _blend_cc(mode::BlendMode, c1, c2, opacity::Real, ::CompositeOperation{:lighter}) =
    mapc((v1, v2) -> _w_safe(v1, oneunit(opacity), v2, opacity), c1, _blend(mode, c1, c2))

@inline function _blend_tc(mode::BlendMode, c1, c2, opacity::Real, op)
    cm = mapch((v1, v2) -> _w(v1, v2, alpha(c1)), c2, _blend(mode, color(c1), c2))
    _comp(op, c1, typeof(c1)(cm, opacity))
end


_blend(::BlendMode{:normal}, c1::C, c2::C) where C <: Color = c2

_blend(::BlendMode{:multiply}, c1::C, c2::C) where C <: Color = mapch(mul, c1, c2)

mul(v1, v2) = v1 * v2
@fastmath function mul(v1::N, v2::N) where N <: Normed
    rv1, rv2, r1 = reinterpret(v1), reinterpret(v2), reinterpret(oneunit(v2))
    m = widemul(rv1, rv2) / convert(floattype(N), r1)
    reinterpret(N, unsafe_trunc(typeof(r1), round(m)))
end
function mul(v1::N0f8, v2::N0f8)
    m = UInt16(reinterpret(v1)) * UInt16(reinterpret(v2))
    z = m + ((m + 0x80) >> 0x8) + 0x80
    reinterpret(N0f8, unsafe_trunc(UInt8, z >> 0x8))
end
function mul(v1::N0f16, v2::N0f16)
    m = UInt32(reinterpret(v1)) * UInt32(reinterpret(v2))
    z = m + ((m + 0x8000) >> 0x10) + 0x8000
    reinterpret(N0f16, unsafe_trunc(UInt16, z >> 0x10))
end

function _blend(::BlendMode{:multiply}, c1::C, c2::C) where {T, C <: Union{Lab{T}, Luv{T}}}
    l = c1.l * c2.l / T(100)
    C(l, mul_lab(comp2(c1), comp2(c2)), mul_lab(comp3(c1), comp3(c2)))
end

function mul_lab(v1::T, v2::T) where T
    t, s = T(128), T(1/128)
    v1r = min(t - v1, t + v1)
    v2r = t - flipsign(v2, v1)
    flipsign(muladd(v1r, -s * v2r, t), v1)
end

_blend(::BlendMode{:screen}, c1::C, c2::C) where C <: Color =
    mapch((v1, v2) -> _n(mul(_n(v1), _n(v2))), c1, c2)

function _blend(::BlendMode{:screen}, c1::C, c2::C) where {T, C <: Union{Lab{T}, Luv{T}}}
    l = muladd(T(100) - c1.l, (c2.l - T(100)) / T(100), T(100))
    C(l, mul_lab(comp2(c1), comp2(c2)), mul_lab(comp3(c1), comp3(c2)))
end

_blend(::BlendMode{:overlay}, c1::C, c2::C) where C <: Color =
    _blend(BlendHardLight, c2, c1)

function _blend(::BlendMode{:overlay}, c1::C, c2::C) where {T, C <: Union{Lab{T}, Luv{T}}}
    l1r = min(c1.l, T(100) - c1.l)
    mr = (l1r + l1r) * ifelse(c1.l == l1r, c2.l, T(100) - c2.l) / T(100)
    l = ifelse(c1.l == l1r, mr, T(100) - mr)
    C(l, mul_lab(comp2(c1), comp2(c2)), mul_lab(comp3(c1), comp3(c2)))
end

_blend(::BlendMode{:darken}, c1::C, c2::C) where C <: Color = mapch(min, c1, c2)

_blend(::BlendMode{:lighten}, c1::C, c2::C) where C <: Color = mapch(max, c1, c2)

_blend(::typeof(BlendColorDodge), c1::C, c2::C) where C <: Color = mapch(dodge, c1, c2)

function dodge(v1, v2)
    v1 == zero(v1) && return v1
    v2 == oneunit(v2) && return v2
    div01(v1, _n(v2))
end

div01(v1, v2) = convert(typeof(v1), min(1, v1 / v2))
@fastmath function div01(v1::N0f8, v2::N0f8)
    rv1, rv2, r1 = reinterpret(v1), reinterpret(v2), reinterpret(oneunit(v2))
    d = Float32(rv1) * r1 / Float32(rv2)
    reinterpret(N0f8, unsafe_trunc(typeof(r1), min(Float32(r1), round(d))))
end
@fastmath function div01(v1::N0f16, v2::N0f16)
    rv1, rv2, r1 = reinterpret(v1), reinterpret(v2), reinterpret(oneunit(v2))
    d = Float64(rv1) * r1 / Float64(rv2)
    reinterpret(N0f16, unsafe_trunc(typeof(r1), min(Float64(r1), round(d))))
end

_blend(::typeof(BlendColorBurn), c1::C, c2::C) where C <: Color = mapch(burn, c1, c2)

# Note that "color-burn" tends to cause the loss of significance.
function burn(v1, v2)
    v1 == oneunit(v1) && return v1
    v2 == zero(v2) && return v2
    _n(div01(_n(v1), v2))
end
burn(v1::Float32, v2::Float32) = Float32(burn(Float64(v1), Float64(v2)))
@fastmath function burn(v1::N0f8, v2::N0f8)
    v1 == oneunit(v1) && return v1
    v2 == zero(v2) && return v2
    rv1, rv2, r1 = reinterpret(v1), reinterpret(v2), reinterpret(oneunit(v2))
    d = Float32(r1) - (Float32(r1) - Float32(rv1)) * r1 / Float32(rv2)
    reinterpret(N0f8, unsafe_trunc(typeof(r1), max(0.0f0, round(d))))
end
@fastmath function burn(v1::N0f16, v2::N0f16)
    v1 == oneunit(v1) && return v1
    v2 == zero(v2) && return v2
    rv1, rv2, r1 = reinterpret(v1), reinterpret(v2), reinterpret(oneunit(v2))
    d = Float64(r1) - (Float64(r1) - Float64(rv1)) * r1 / Float64(rv2)
    reinterpret(N0f16, unsafe_trunc(typeof(r1), max(0.0, round(d))))
end

_blend(::typeof(BlendHardLight), c1::C, c2::C) where C <: Color = mapch(hardlight, c1, c2)

function hardlight(v1, v2)
    v2r = min(v2, _n(v2))
    mr = mul(v2r + v2r, ifelse(v2 == v2r, v1, _n(v1)))
    ifelse(v2 == v2r, mr, _n(mr))
end

_blend(::typeof(BlendHardLight), c1::C, c2::C) where {T, C <: Union{Lab{T}, Luv{T}}} =
    _blend(BlendOverlay, c2, c1)

_blend(::typeof(BlendSoftLight), c1::C, c2::C) where C <: Color = mapch(softlight, c1, c2)

function softlight(v1, v2)
    v2r = min(v2, _n(v2))
    nv2r2 = _n(v2r + v2r)
    if v2r == v2 # if v2 <= 0.5
        return mul(v1, _n(mul(nv2r2, _n(v1))))
    else
        m = mul(v1, _n(v2))
        @fastmath d(v) = v <= 0.25 ? muladd(4v - 3, v, 1) * 4v : sqrt(v)
        return m + m + mul(nv2r2, convert(typeof(v1), d(v1)))
    end
end
softlight(v1::X, v2::X) where X <: FixedPoint = softlight(float(v1), float(v2)) % X

function _blend(::typeof(BlendSoftLight), c1::C, c2::C) where {T, C <: Union{Lab{T}, Luv{T}}}
    l = softlight(c1.l / T(100), c2.l / T(100))
    a = softlight(muladd(comp2(c1), T(1/256), T(0.5)),
                  muladd(comp2(c2), T(1/256), T(0.5)))
    b = softlight(muladd(comp3(c1), T(1/256), T(0.5)),
                  muladd(comp3(c2), T(1/256), T(0.5)))
    C(l * T(100), muladd(a, T(256), T(-128)) , muladd(b, T(256), T(-128)))
end

_blend(::BlendMode{:difference}, c1::C, c2::C) where C <: Color = mapch(difference, c1, c2)

difference(v1, v2) = max(v1, v2) - min(v1, v2)

_blend(::BlendMode{:exclusion}, c1::C, c2::C) where C <: Color = mapch(exclusion, c1, c2)

function exclusion(v1, v2)
    m = mul(v1, v2)
    (v1 - m) + (v2 - m)
end
function exclusion(v1::N0f8, v2::N0f8)
    rv1, rv2 = UInt16(reinterpret(v1)), UInt16(reinterpret(v2))
    m = rv1 * rv2
    z = muladd(rv1 + rv2, 0xFF, -(m + m))
    reinterpret(N0f8, unsafe_trunc(UInt8, (z + ((z + 0x80) >> 0x8) + 0x80) >> 0x8))
end
function exclusion(v1::N0f16, v2::N0f16)
    rv1, rv2 = UInt32(reinterpret(v1)), UInt32(reinterpret(v2))
    m = rv1 * rv2
    z = muladd(rv1 + rv2, 0xFFFF, -(m + m))
    reinterpret(N0f16, unsafe_trunc(UInt16, (z + ((z + 0x8000) >> 0x10) + 0x8000) >> 0x10))
end

lum100(c::AbstractRGB) = lum100(red(c), green(c), blue(c))
lum100(r::T, g::T, b::T) where T = muladd(T(59), g, muladd(T(11), b, T(30) * r))
function lum100(r::T, g::T, b::T) where T <: Normed
    F = floattype(T)
    rr = reinterpret(r)
    rg = reinterpret(g)
    rb = reinterpret(b)
    ri = reinterpret(oneunit(T))
    muladd(F(59), rg, muladd(F(11), rb, F(30) * rr)) / F(ri)
end
function lum100(r::T, g::T, b::T) where T <: Normed{UInt8}
    F = floattype(T)
    rr = reinterpret(r)
    rg = reinterpret(g)
    rb = reinterpret(b)
    unsafe_trunc(Int32, 0x268cf3 * rr + 0x4bd0f0 * rg + 0xe229d * rb) * 4.65661322835f-8
end
function setlum(c::C, l::AbstractFloat) where {T, C <: AbstractRGB{T}}
    lc = lum100(c)
    r, g, b = red(c), green(c), blue(c)
    n = min(r, g, b)
    x = max(r, g, b)

    d = (lc - l) / 100
    if n < d
        # k = l / (lc - 100n)
        # cr, cg, cb = map(c -> c * k - n * k, (r, g, b))
        if n == r
            cr = zero(T)
            cb = (b - n) * l / muladd(T(-70), r - g, T(11) * (b - g))
            cg = max(cr, muladd(T(-11), cb, l) * T(1/59))
        elseif n == g
            cg = zero(T)
            cb = (b - n) * l / muladd(T(30), r - g, T(11) * (b - g))
            cr = max(cg, muladd(T(-11), cb, l) * T(1/30))
        else
            cb = zero(T)
            cr = (r - n) * l / muladd(T(30), r - g, T(-89) * (b - g))
            cg = max(cb, muladd(T(-30), cr, l) * T(1/59))
        end
        return C(cr, cg, cb)
    elseif d < x - oneunit(x)
        # k1 = 1 / (100x - lc)
        # kl = l / (100x - lc)
        # mapc(v -> x * kl - v * kl + (100v - lc) * k1, c)
        if x == r
            k1 = muladd(T(70), r - g, T(-11) * (b - g))
            k2 = muladd(T(-30), r - g, T(89) * (b - g))
            cr = oneunit(T)
            cb = muladd(x - b, l, k2) / k1
            cg = min(cr, muladd(T(-11), cb, l - T(30)) * T(1/59))
        elseif x == g
            k1 = muladd(T(-30), r - g, T(-11) * (b - g))
            k2 = muladd(T(-30), r - g, T(89) * (b - g))
            cg = oneunit(T)
            cb = muladd(x - b, l, k2) / k1
            cr = min(cg, muladd(T(-11), cb, l - T(59)) * T(1/30))
        else
            k1 = muladd(T(-30), r - g, T(89) * (b - g))
            k2 = muladd(T(70), r - g, T(-11) * (b - g))
            cb = oneunit(T)
            cr = muladd(x - r, l, k2) / k1
            cg = min(cb, muladd(T(-30), cr, l - T(11)) * T(1/59))
        end
        return C(cr, cg, cb)
    end
    return mapc(v -> v - d, c)
end
function setlum(c::C, l::AbstractFloat) where {T <: Normed, C <: AbstractRGB{T}}
    convert(C, setlum(convert(RGB{floattype(T)}, c), l))
end


sat(c::AbstractRGB) = sat(red(c), green(c), blue(c))
sat(r::T, g::T, b::T) where T = max(r, g, b) - min(r, g, b)
function setsat(c::C, s) where C <: AbstractRGB
    r, g, b = float(red(c)), float(green(c)), float(blue(c))
    cmin_rg, cmax_rg = minmax(r, g)
    cmin = min(cmin_rg, b)
    cmax = max(cmax_rg, b)
    cmax == cmin && return C(zero(eltype(C)))

    cmid = max(cmin_rg, min(cmax_rg, b))

    cmid2 = ((cmid - cmin) * float(s)) / (cmax - cmin)
    cmax2 = float(s)

    r2 = ifelse(cmax == r, cmax2, ifelse(cmin == r, zero(r), cmid2))
    g2 = ifelse(cmax == g, cmax2, ifelse(cmin == g, zero(g), cmid2))
    b2 = ifelse(cmax == b, cmax2, ifelse(cmin == b, zero(b), cmid2))

    C(r2, g2, b2)
end

_blend(m::M, c1::C, c2::C) where {C <: Color, M <: NonSeparableBlendMode} =
    convert(C, _blend(m, convert(RGB, c1), convert(RGB, c2)))

_blend(::BlendMode{:hue}, c1::C, c2::C) where C <: AbstractRGB =
    setlum(setsat(c2, sat(c1)), lum100(c1))

function _blend(::BlendMode{:hue}, c1::C, c2::C) where C <: Union{Lab, Luv}
    if abs(comp2(c2)) >= abs(comp3(c2))
        a = copysign(max(abs(comp2(c1)), abs(comp3(c1))), comp2(c2))
        b = ifelse(a == zero(a), a, a * comp3(c2) / comp2(c2))
    else
        b = copysign(max(abs(comp2(c1)), abs(comp3(c1))), comp3(c2))
        a = ifelse(b == zero(b), b, b * comp2(c2) / comp3(c2))
    end
    C(c1.l, a, b)
end

_blend(::BlendMode{:hue}, c1::C, c2::C) where C <: Union{HSV, HSL, HSI} =
    C(c2.h, c1.s, comp3(c1))

_blend(::BlendMode{:saturation}, c1::C, c2::C) where C <: AbstractRGB =
    setlum(setsat(c1, sat(c2)), lum100(c1))

function _blend(::BlendMode{:saturation}, c1::C, c2::C) where C <: Union{Lab, Luv}
    if abs(comp2(c1)) >= abs(comp3(c1))
        a = copysign(max(abs(comp2(c2)), abs(comp3(c2))), comp2(c1))
        b = ifelse(a == zero(a), a, a * comp3(c1) / comp2(c1))
    else
        b = copysign(max(abs(comp2(c2)), abs(comp3(c2))), comp3(c1))
        a = ifelse(b == zero(b), b, b * comp2(c1) / comp3(c1))
    end
    C(c1.l, a, b)
end

_blend(::BlendMode{:saturation}, c1::C, c2::C) where C <: Union{HSV, HSL, HSI} =
    C(c1.h, c2.s, comp3(c1))

_blend(::BlendMode{:color}, c1::C, c2::C) where C <: AbstractRGB =
    setlum(c2, lum100(c1))

_blend(::BlendMode{:color}, c1::C, c2::C) where C <: Union{Lab, Luv} =
    C(c1.l, comp2(c2), comp3(c2))

_blend(::BlendMode{:color}, c1::C, c2::C) where C <: Union{HSV, HSL, HSI} =
    C(c2.h, c2.s, comp3(c1))

_blend(::BlendMode{:luminosity}, c1::C, c2::C) where C <: AbstractRGB =
    setlum(c1, lum100(c2))


_blend(::BlendMode{:luminosity}, c1::C, c2::C) where C <: Union{Lab, Luv} =
    C(c2.l, comp2(c1), comp3(c1))

_blend(::BlendMode{:luminosity}, c1::C, c2::C) where C <: Union{HSV, HSL, HSI} =
    C(c1.h, c1.s, comp3(c2))
