Problem with conversion of YIQ (or problem with non-sRGB gamut)

[kimikage]

I have a proposal to add abstract gamut types in ColorTypes.jl.(JuliaGraphics/ColorTypes.jl#140)
I think it will help reduce the misunderstandings about the color gamut and the misfortunes from them.

However, it also clarifies the issues that have been obscured so far. One of them is a problem with conversion of YIQ.

The color primaries of NTSC and sRGB are different, and the NTSC color gamut is actually wider than sRGB's one.

_Originally posted by @kimikage in JuliaGraphics/ColorTypes.jl#140 (comment)

In fact, the current RGB is implicitly considered as "sRGB". On the other hand, the current conversion methods from/to YIQ consider that the RGB is defined under the "NTSC" gamut.

Colors.jl/src/conversions.jl

Lines 712 to 723 in 1bc194e

	correct_gamut(c::YIQ{T}) where {T} = YIQ{T}(clamp(c.y, zero(T), one(T)),
	clamp(c.i, convert(T,-0.5957), convert(T,0.5957)),
	clamp(c.q, convert(T,-0.5226), convert(T,0.5226)))
	function cnvt(::Type{YIQ{T}}, c::AbstractRGB) where T
	rgb = correct_gamut(c)
	YIQ{T}(0.299*red(rgb)+0.587*green(rgb)+0.114*blue(rgb),
	0.595716*red(rgb)-0.274453*green(rgb)-0.321263*blue(rgb),
	0.211456*red(rgb)-0.522591*green(rgb)+0.311135*blue(rgb))
	end
	cnvt(::Type{YIQ{T}}, c::Color3) where {T} = cnvt(YIQ{T}, convert(RGB{T}, c))

The trouble is that the gamuts are often ignored when we use YIQ or YCbCr. In fact, some image processing tools do such conversions which are inaccurate from a color management perspective. However, in the testing, the inaccuracy or inconsistency is a big problem.

[kimikage]

I don't think YIQ is used so often. Therefore, it may be a good idea to add the gamut conversion to the conversion methods, and add low-level conversion methods like rgb_to_yiq_ntsc and yiq_to_rgb_ntsc.

Again, this is because I think that YIQ is not used much. The color gamut issues potentially have a similar problem to the white points.
#278 (comment)

Edit:
What makes it even more complicated is that "NTSC" sometimes means "SMPTE-C". https://en.wikipedia.org/wiki/NTSC#Colorimetry
As shown below, "SMPTE-C" is similar to "sRGB" in numerical terms, but their context is completely different.

primary_red(  ::Type{Gamut_sRGB}) = xyY(0.64, 0.33, 1)
primary_green(::Type{Gamut_sRGB}) = xyY(0.30, 0.60, 1)
primary_blue( ::Type{Gamut_sRGB}) = xyY(0.15, 0.06, 1)
whitepoint(::Type{Gamut_sRGB}) = WP_D65

# original NTSC(1953) gamut
abstract type Gamut_NTSC <: AbstractRGBGamut end
primary_red(  ::Type{Gamut_NTSC}) = xyY(0.67, 0.33, 1)
primary_green(::Type{Gamut_NTSC}) = xyY(0.21, 0.71, 1)
primary_blue( ::Type{Gamut_NTSC}) = xyY(0.14, 0.08, 1)
whitepoint(::Type{Gamut_NTSC}) = WP_C

# SMPTE-C gamut
abstract type Gamut_SMPTE_C <: AbstractRGBGamut end
primary_red(  ::Type{Gamut_SMPTE_C}) = xyY(0.630, 0.340, 1)
primary_green(::Type{Gamut_SMPTE_C}) = xyY(0.310, 0.595, 1)
primary_blue( ::Type{Gamut_SMPTE_C}) = xyY(0.155, 0.070, 1)
whitepoint(::Type{Gamut_SMPTE_C}) = WP_D65

function mat_rgb_to_xyz(gamut::Type{<:AbstractRGBGamut},
                        wp::Union{XYZ, xyY}=whitepoint(gamut))
    pr, pg, pb = primary_red(gamut), primary_green(gamut), primary_blue(gamut)
    z(c::xyY) = 1 - c.x - c.y # Y == 1
    m_prim = [ pr.x  pg.x  pb.x
               pr.y  pg.y  pb.y
               z(pr) z(pg) z(pb) ]
    w = convert(XYZ, wp)
    sr, sg, sb = inv(m_prim) * [w.x, w.y, w.z] # diag.
    @inbounds [ m_prim[1,1]*sr m_prim[1,2]*sg m_prim[1,3]*sb
                m_prim[2,1]*sr m_prim[2,2]*sg m_prim[2,3]*sb
                m_prim[3,1]*sr m_prim[3,2]*sg m_prim[3,3]*sb ]
end

function mat_xyz_to_rgb(gamut::Type{<:AbstractRGBGamut},
                        wp::Union{XYZ, xyY}=whitepoint(gamut))
    inv(mat_rgb_to_xyz(gamut, wp))
end

function convert_gamut(c::T,
                       src::Type{<:AbstractRGBGamut},
                       dest::Type{<:AbstractRGBGamut}) where {T <: AbstractRGB}
    # awesome RGB -> XYZ -> RGB conversion 
end

[kimikage]

cf. #8

[kimikage]

I will change the gamut for YIQ to "SMPTE-C" instead of "sRGB" or (original) "NTSC".
We can use yiq_to_rgb_smptec/rgb_to_yiq_smptec as the compatible methods without gamut conversion.

However, the gamma and the narrow sense of the color gamut are independent concepts. There is already a lot of confusion. 😕

[kimikage]

I decided not to change this unless anyone else wanted the fix in Colors.jl.
cf. #412 (comment)

[kimikage]

FYI, the following is the conversion code from sRGB to "SMPTE 170M" YIQ.

function mat_srgb_to_xyz()
    m_prim = [640 300 150
              330 600  60
               30 100 790] ./ big"1000"
    w = [big"0.95047", big"1.0", big"1.08883"]
    sr, sg, sb = inv(m_prim) * w # diag.
    to_xyz = [ m_prim[1,1]*sr m_prim[1,2]*sg m_prim[1,3]*sb
               m_prim[2,1]*sr m_prim[2,2]*sg m_prim[2,3]*sb
               m_prim[3,1]*sr m_prim[3,2]*sg m_prim[3,3]*sb ]
    return to_xyz
end

function mat_xyz_to_smptergb()
    m_prim = [630 310 155
              340 595  70
               30  95 775] ./ big"1000"
    w = [big"0.95047", big"1.0", big"1.08883"]
    sr, sg, sb = inv(m_prim) * w # diag.
    to_xyz = [ m_prim[1,1]*sr m_prim[1,2]*sg m_prim[1,3]*sb
               m_prim[2,1]*sr m_prim[2,2]*sg m_prim[2,3]*sb
               m_prim[3,1]*sr m_prim[3,2]*sg m_prim[3,3]*sb ]
    return inv(to_xyz)
end

const M_SRGB_TO_XYZ = mat_srgb_to_xyz()
const M_XYZ_TO_SMPTERGB = mat_xyz_to_smptergb()
const M_SRGB_TO_SMPTERGB = M_XYZ_TO_SMPTERGB * M_SRGB_TO_XYZ

srgb_to_lin(v) = v <= big"0.04045" ? v / big"12.92" : ((v + big"0.055") / big"1.055")^big"2.4"
smpte_to_lin(v) = v < big"0.0812" ? v / big"4.5" : ((v + big"0.099") / big"1.099")^(100/big"45") # similar to Rec.709

function srgb_to_yiq(rgb::RGB)
    r0 = srgb_to_lin(rgb.r)
    g0 = srgb_to_lin(rgb.g)
    b0 = srgb_to_lin(rgb.b)
    rgb1 = M_SRGB_TO_SMPTERGB * [r0, g0, b0]
    r1 = rgb1[1]
    g1 = rgb1[2]
    b1 = rgb1[3]
    y = (587 * g1 + 114 * b1 + 299 * r1) / big"1000"
    kb = sqrt(big"209556997" / big"96146491") / big"3"
    kr = sqrt(big"221990474" / big"288439473")
    u = kb * (b1 - y)
    v = kr * (r1 - y)
    i = -u * sind(big"33") + v * cosd(big"33")
    q =  u * cosd(big"33") + v * sind(big"33")
    YIQ{Float64}(y, i, q)
end