Important decisions with respect to color math (please comment)

[timholy]

RGB multiplication and implications for var (#119, #125)

Currently multiplication of two RGBs is undefined, but sometimes this causes problems (#119, "tinting" as in https://discourse.julialang.org/t/loaderror-dimensionmismatch/35507?u=tim.holy). In #119, we hit on the possibility of defining ⋅ (\cdot) as the "dot product" (with channels viewed as a vector), * operating channelwise, and ⊗ (\otimes) some day being the outer product (though we currently lack a type to express this result). I'll modify this proposal below, but for now let's run with it and explore the implications. An important route lies through abs2...

Currently, abs2(c) returns a scalar, which in this proposed framework would be equivalent to defining abs2(c) = c⋅c. This is similar to how abs2 behaves with respect to the real and imaginary components of a complex number. This results in the following:

julia> using ColorVectorSpace, Statistics

julia> r = rand(RGB{Float32}, 2)
2-element Array{RGB{Float32},1} with eltype RGB{Float32}:
 RGB{Float32}(0.6334499f0,0.34752238f0,0.17738128f0)
 RGB{Float32}(0.35009634f0,0.38727677f0,0.6068548f0)

julia> var(r)
0.13315858f0

However, if you load Images it changes the behavior of var to be channelwise, because Images specializes var for arrays of Colorants ☹️ . Rather than overloading var for certain element types, we could let * be channelwise multiplication and define abs2(c) = c*c, and it should even make it robust with respect to the dims keyword argument (#125, CC @wizofe).

However, since this package is called ColorVectorSpace it's also worth noting this comparison:

julia> using Statistics, StaticArrays

julia> A = rand(3, 5)
3×5 Array{Float64,2}:
 0.901707   0.763397  0.0841844  0.164515  0.542334
 0.247929   0.690001  0.492793   0.275766  0.175159
 0.0860073  0.120248  0.0454239  0.403493  0.529168

julia> var(A; dims=2)
3×1 Array{Float64,2}:
 0.129401636194079  
 0.04475575432189567
 0.04655327916103   

julia> a = reinterpret(SVector{3,Float64}, A)
1×5 reinterpret(SArray{Tuple{3},Float64,1,3}, ::Array{Float64,2}):
 [0.901707, 0.247929, 0.0860073]  …  [0.542334, 0.175159, 0.529168]

julia> var(a; dims=2)
ERROR: MethodError: no method matching abs2(::SArray{Tuple{3},Float64,1,3})
Closest candidates are:
  abs2(::Missing) at missing.jl:100
  abs2(::Bool) at bool.jl:84
  abs2(::Real) at number.jl:157
  ...
Stacktrace:
 [1] (::Statistics.var"#12#13")(::SArray{Tuple{3},Float64,1,3}) at /home/tim/src/julia-1/usr/share/julia/stdlib/v1.3/Statistics/src/Statistics.jl:300
...

so this is an issue that doesn't have an agreed-upon resolution elsewhere in the ecosystem. We might actually want to be able to pass an operator, and the thinking above inspires me to propose var(⋅, a), var(*, a), and var(⊗, a).

However, this conflicts with the very clear decision that u*v for two vectors u and v is deliberately undefined, and specifically does not mean elementwise multiplication (that's what broadcasting is for). Unfortunately, it's not possible to pass broadcasted-* (.*) as an argument to a function. Consequently we might consider choosing an operator to be a synonym for elementwise multiplication, and leaving * undefined for RGB. Options include ∗ (\ast), • (\bullet), and ⊙ (\odot). I worry that \ast looks too similar to * and would be confused for it; if we go this way I think my preference is for \odot but curious to see what others think.

By this token abs2 should be undefined for colors, var(a) should throw an error for color arrays, but var(⊙, a) should do what users want.

Because of the analogy with arrays-of-vectors I'll open a similar issue in Statistics.jl and see what folks outside of the JuliaImages world think. (Update: https://github.com/JuliaLang/Statistics.jl/issues/29)

Interaction between color math and conversions

Another factor that would argue in favor of changing the behavior of abs2 is the following set of logical-sounding conclusions:

• if x is a real number, then Gray(x) is basically the same thing (other than being a display hint)
• RGB(Gray(x)) should return RGB(x, x, x), so by the same token RGB(x, x, x) is essentially equivalent to x

Our current scheme has a big gotcha:

julia> x = 0.5
0.5

julia> abs2(x)
0.25

julia> abs2(RGB(x, x, x))
0.75

so we're in a situation where a and b can be equivalent but abs2(a) is very different from abs2(b). If we changed abs2 to behave channelwise, we'd restore consistency. An error would also fix it; we don't object to -1 == -1+0im but having sqrt return an error for just one of them.

Result types from arithmetic (#38)

I originally designed the rules to be inspired to unitful arithmetic. However, most others who have commented seem to expect that colors will be "poisoning" in much the same way that NaN is poisoning in arithmetic. (#38 (comment) as well as the OP in #38).

Given that we support Gray(0.1)^2 and there is no type encoding Gray^2, I now think the poisoning metaphor is overall the better choice.

I suspect there is one interesting exception to the poisoning rule, however, which comes from linear algebra:

julia> v = SVector((0.2, 0.1, 0.8))
3-element SArray{Tuple{3},Float64,1,3} with indices SOneTo(3):
 0.2
 0.1
 0.8

julia> u = SVector((0.1, 0.2, 0.8))
3-element SArray{Tuple{3},Float64,1,3} with indices SOneTo(3):
 0.1
 0.2
 0.8

julia> u \ v
0.9855072463768115

Note the returned value is a scalar, not an SVector{3}. (The value minimizes norm(α*u - v) for real-valued α.) Should we define \ for Gray and RGB? / cannot be defined for RGB (EDIT: in the denominator), but it can for Gray; should Gray(0.2) / Gray(0.1) likewise be "colorless"? That is a little bit more like the unitful design. Is that a good thing, or would that be annoying?

[timholy]

This does not seem to have generated much controversy, so I'll begin implementing these changes. I will use \odot for elementwise multiplication.

[kimikage]

I have no objection about the multiplication.The multiplication of AbstractVector{<:Real} and Gray (or RGB) is controversial.

Regarding var or abs2, it is better to keep the backward compatibility for now.
TBH, I hate the implementation of abs and Normed-->Float32 (not floattype()) conversions below, but that is off-topic.:ghost:

ColorVectorSpace.jl/src/ColorVectorSpace.jl

Lines 160 to 166 in f2baf1d

	abs(c::AbstractRGB) = abs(red(c))+abs(green(c))+abs(blue(c)) # should this have a different name?
	abs(c::AbstractRGB{T}) where {T<:Normed} = Float32(red(c))+Float32(green(c))+Float32(blue(c)) # should this have a different name?
	abs(c::TransparentRGB) = abs(red(c))+abs(green(c))+abs(blue(c))+abs(alpha(c)) # should this have a different name?
	abs(c::TransparentRGB{T}) where {T<:Normed} = Float32(red(c))+Float32(green(c))+Float32(blue(c))+Float32(alpha(c)) # should this have a different name?
	abs2(c::AbstractRGB) = red(c)^2+green(c)^2+blue(c)^2
	abs2(c::AbstractRGB{T}) where {T<:Normed} = Float32(red(c))^2+Float32(green(c))^2+Float32(blue(c))^2
	abs2(c::TransparentRGB) = (ret = abs2(color(c)); ret + convert(typeof(ret), alpha(c))^2)

[timholy]

it is better to keep the backward compatibility for now

I'm eager to get this right, and so the sooner we fix problems the better. But I'm all in favor of a staged migration with good deprecation warnings. I haven't yet begun to plan for how to do the migration; I'm more interested in finding out whether the considerations I outline in the OP represent where we want to end up. For abs and abs2, the key issue is not the precise numeric type, instead it's deciding whether it should behave channelwise or return a single real-valued scalar. I.e., are they effectively norm in the same as they work for complex numbers, or are they vector-valued?

[johnnychen94]

If I understand it right, the motivation you propose var(*, a) is that you want to make it less ambiguity, just like how we did in Colors and FixedPointNumbers to not guess color types and storage types.

If we're thinking RGB pixel as a point in a vector space, then we don't have definitions like abs2(RGB(x, x, x)), instead it should be norm(2, RGB(x, x, x)). Abusing the name var and abs2 doesn't look like a good solution to me.

It's of conceptual interest, but I doubt that var(::Colorant) has its usage in practice compared to var(::AbstractArray{<:Colorant}).

Although it would be less convenient, making them vector-valued is probably the clearest way without confusion. It would be better to leave this level of decision to downstream packages whether they need to sum these vectors up or use the channel-wise results.

During my time developing ImageDistances and ImageQualityIndexes, I was having similar pains. -- For example, my mind blew when I found out PSNR of RGB image should be a real number while PSNR of other color images should be a vector, and that Minkowski(2) != Euclidean 😖

P.S. This issue is quite too general for me to find a way out, and my words here might not make much sense. But I feel a need to comment on something before you start playing with this.

[kimikage]

BTW, when var for RGB (i.e. not perceptually uniform) images is channel-wise, var and cov should be a pair from the standpoint of multivariate analysis. However, in practice, cov is usually ignored.:confused:

[timholy]

With this proposal, we won't even need cov because it will be (in your notation)

var((x, μ)->(x - μ)⊗(x - μ), collection)

[kimikage]

Abusing the name var and abs2 doesn't look like a good solution to me.

I agree to some degree. However, abs2 and 2-norm are confused. In particular, norm is in LinearAlgebra and abs2 is in Base, so it seems that abs2 is sometimes used in the context where norm should be used. In fact, the current Statistics do not use norm, so var (with dims) for vectors is broken. As you mentioned, at the price of strictness, we lose convenience. (I also love strictness, though.)

WRT var, we don't have to define var(::Colorant), even though a Colorant is like a vector.

[kimikage]

BTW, I'm going to create a new package ColorBlendModes, which implements the blend modes specified in the Compositing and Blending Level 1 (cf. JuliaGraphics/ColorTypes.jl#167).
TBH, I don't really want the functionality. However, I believe the new package will clarify the roles of other packages. For example, since the blend modes includes multiply, we will not need to fully support the multiply in the context of "compositing" in ColorVectorSpace.
ColorBlendModes will not define any arithmetic operators.

Edit:
WIP: https://kimikage.github.io/ColorBlendModes.jl/dev/
Edit2:
ColorBlendModes.jl has been registered.
xref: https://discourse.julialang.org/t/ann-colorblendmodes-jl/37480