An HDR PIxel Type

This is an implementation of a RGB/XYZ pixel format commonly used to store HDR lumiance maps. It falls somewhere between visualization and modeling & simulation in subject area, and so I am posting it here, in the hope of some commentary.

"""
    RealPixel{N}

HDR pixel format for radiance and luminance, as used in Radiance, and
widely used for HDR images.

This is an implementation of Greg Ward's "Real Pixels"[^rp] floating
point pixel compression scheme, intended for a compact representation
of absolute RGB or XYZ luminance values. It follows the *Radiance*
`color.c` code.  The code has been extended to arbitrarily many
spectral bands, supporting the new *Radiance* multispectral format.

[^rp]: Ward, Greg. “Real Pixels.” In *Graphics Gems II*, edited by
Arvo, James, 80–83. Graphics Gems Series. Boston: Academic Press,
1991.

The structure defines an n-byte mantissa, followed by a single byte
exponent.  Since physical values of luminance and radiance are always
positive and XYZ values are normalized to be positive the mantissas
are unsigned quantites.  The exponent is an unsigned representation of
a signed quantity; it is the signed quantity + 128.  In the important
special case of three spectral channels, RGB or XYZ, a RealPixel fits
in a 32-bit word.

"""
struct RealPixel{N} <: AbstractVector{Real}
    mantissa::NTuple{N,UInt8}
    exponent::UInt8
end

"""
    RealPixel(v::Vector{<:Real})

Construct a RealPixel from a vector of Reals, by compressing the Reals
into the RealPixel format.

"""
function RealPixel(v::Vector{<:Real})
    RealPixel(realpixelfromfloat(v))
end

"""
    RealPixel(v::Vector{UInt8})

Construct a RealPixel from raw UInt8s.  The last UInt8 is the
exponent.

"""
function RealPixel(v::Vector{UInt8})
    length(v) < 1 && error("RealPixel needs at least one UInt8")
    RealPixel(Tuple(v[begin:end-1]), v[end])
end

# AbstractVector interface
"""
    size(px::RealPixel)

The dimensions of the vector.

"""
size(px::RealPixel) = (length(px.mantissa),)

"""
    getindex(px::RealPixel,i::Int)::Float32

Subscript a RealPixel.

"""
function getindex(px::RealPixel,i::Int)::Float32
    1 <= i <= length(px) || throw(BoundsError(px, i))
    px.exponent == 0 && return 0.0
    return ldexp(px.mantissa[i] + 0.5, px.exponent - 136)
end

"""
    iterate(px::RealPixel, i=1)

Iterate through the elements of a RealPixel

"""
iterate(px::RealPixel, i=1) = (i <= length(px)) ? (px[i], i+1) : nothing

"""
    convert(T::Type{Vector}, px::RealPixel)::AbstractVector{Float32}

Convert a RealPixel to a vector of Float32 by decompression.

"""
convert(T::Type{Vector}, px::RealPixel)::AbstractVector{Float32} =
    floatfromrealpixel(vcat(collect(px.mantissa), px.exponent))


"""
    floatfromrealpixel(px::AbstractVector{UInt8})

Expand an 8(n+1)-bit real pixel to a 32n-bit pixel.

"""
function floatfromrealpixel(px::AbstractVector{UInt8})::Vector{Float32}
    exponent = px[end]
    mantissas = px[1:end-1]
    exponent == 0 && return zeros(Float32, length(mantissas))
    return [ldexp(v + 0.5, exponent - 136) for v in mantissas]
end

"""
    realpixelfromfloat(f::AbstractVector{<:Real})

Compress a pixel of some real type to an 8(n+1)-bit real pixel.
Returns a vector of UInt8, where the last value is the exponent.

"""
function realpixelfromfloat(f::AbstractVector{<:Real})
    maxchannel = maximum(f)
    maxchannel < 1e-32 && return zeros(UInt8, length(f)+1)
    x,e = frexp(maxchannel)
    d = x * 256.0 / maxchannel
    rp = [trunc(UInt8, v * d) for v in f]
    push!(rp, UInt8(e + 128))
    return rp
end
1 Like

Hi there :slight_smile: Looks already very promising. Reading through the code I noticed a couple of small things. Here are some comments in no particular order:

Docstrings

These docstrings are not very helpful. Of course iterate iterates through the elements of a RealPixel - that’s what iterate is supposed to do after all :smiley:
In Julia it is perfectly fine to not write docstrings for every method and instead have one applying to the whole function. You only should add docstring for methods if the method does something more specific than the what the function in general does. Here is the example of getindex:

help?> getindex
search: getindex setindex! lastindex firstindex checkindex Base.to_index

  getindex(A, inds...)

  Return a subset of array A as selected by the indices inds.

  Each index may be any supported index type, such as an Integer,
  CartesianIndex, range, or array of supported indices. A : may be used to
  select all elements along a specific dimension, and a boolean array (e.g. an
  Array{Bool} or a BitArray) may be used to filter for elements where the
  corresponding index is true.

  When inds selects multiple elements, this function returns a newly allocated
  array. To index multiple elements without making a copy, use view instead.

  See the manual section on array indexing for details.
[...]

  getindex(collection, key...)

  Retrieve the value(s) stored at the given key or index within a collection.
  The syntax a[i,j,...] is converted by the compiler to getindex(a, i, j, ...).
[...]

You see that the second docstring is the general one and the first one is specific to the context of AbstractArrays (actually this is not mentioned explicitly and that could perhaps be improved slightly but I hope you get the point).

Internal interface

I consider this interface where you store everything in a vector but the last element means something different entirely to be somewhat bad because it requires implicit knowledge which can cause mistakes easily. In Julia, it is perfectly fine to return multiple values from a function, so there is no need to pack the exponent in with the mantissas. This also makes your code more concise everywhere. Concrete examples (I cut the docstrings for brevity):

function realpixelfromfloat(f::AbstractVector{<:Real})
    maxchannel = maximum(f)
    maxchannel < 1e-32 && return zeros(UInt8, length(f)+1)
    x,e = frexp(maxchannel)
    d = x * 256.0 / maxchannel
    return [trunc(UInt8, v * d) for v in f], UInt8(e + 128) # returns 2 values
end

function RealPixel(v::Vector{<:Real})
    RealPixel(realpixelfromfloat(v)...) # needs slight adaption
end

function floatfromrealpixel(px::AbstractVector{UInt8}, exponent::UInt8)
    exponent == 0 && return zeros(Float32, length(mantissas))
    return [ldexp(v + 0.5, exponent - 136) for v in mantissas]
    # or
    # return @. ldexp(mantissas + 0.5, exponent - 136)
end

Also is there a reason you have realpixelfromfloat? This could just be a constructor I think. That would also be a bit more efficient since you can skip allocating an array and use a NTuple directly:

function RealPixel(f::AbstractVector{<:Real})
    maxchannel = maximum(f)
    maxchannel < 1e-32 && return zeros(UInt8, length(f)+1)
    x,e = frexp(maxchannel)
    d = x * 256.0 / maxchannel
    mantissas = Tuple(trunc(UInt8, v * d) for v in f)
    exponent = UInt8(e + 128)
    return RealPixel(mantissas, exponent)
end

Possible mistake?

I don’t know the underlying math/algorithm and also didn’t run your code but I noticed an asymmetry between floatfromrealpixel and realpixelfromfloat: In realpixelfromfloat you have UInt8(e + 128) as exponent while in floatfromrealpixel you do exponent - 136 and I am wondering a bit why these magic numbers are different.

Unnecessary implementation of iterate

Once you implement Base.size and Base.getindex you get iterate for free (for AbstractArrays not for general types).

A note on type stability

All constructors you defined are inherently type unstable because RealPixe{N} carries its length as type parameter but Vector{<:Real} which you use for construction does not. This may or may not be an issue for you. It is fine if:

  • You rarely convert vectors to RealPixel and/or do so only in the beginning of your programs such that the type instability is confined to a small portion of the code
  • You never have heterogeneous collections of RealPixels of different lengths

If you either convert often, do it at several points throughout your programs or simply don’t want to care you could consider changing RealPixel to:

struct RealPixel <: AbstractVector{Real}
    mantissa::Vector{UInt8}
    exponent::UInt8
end

Alternatively, if you absolutely need/crave the performance you could consider adding an optional argument to the constructors/conversion functions:

function RealPixel(f::AbstractVector{<:Real}, ::Val{N}=Val(length(f))) where N
    length(f) < N && error("Given to few datapoints for length $N")
    maxchannel = maximum(f)
    maxchannel < 1e-32 && return zeros(UInt8, length(f)+1)
    x,e = frexp(maxchannel)
    d = x * 256.0 / maxchannel
    mantissas = NTuple{N,UInt8}(trunc(UInt8, v * d) for v in f)
    exponent = UInt8(e + 128)
    return RealPixel(mantissas, exponent)
end

Then if you know the length beforehand you can construct a RealPixel like RealPixel(data, Val(3)).

4 Likes

Since getindex(::RealPixel, ::Int) returns a Float32, this should be

struct RealPixel{N} <: AbstractVector{Float32}
...

That will make eltype(px) return Float32.

That said, it may be worth considering whether you really want to make these iterable over color channels. In the ColorTypes/Colors ecosystem, for better or worse we decided not to. In my experience there are advantages either way: broadcasting is often easier if they are scalars (non-iterable), but people also need to be able address individual color channels. You can get the latter with a set of utility functions. But there are also advantages in iterability.

2 Likes

Thank you; I’ll do that. They are actually sort-of Float8 and I debated returning Float16 instead, but I suspect that would cost speed for a little bit of clarity.

I did consider iterability; multispectral pixels tipped the balance in favor of it. In studies, Greg Ward found that 12-18 bands was appropriate for visible light; he implemented 24 in Radiance to accommodate IR and UV bands. It seems to me that that many spectral bands made iterability desirable.

Thanks for the detailed thoughtful comments; I’m busy with other things right now, but I will come back to this & incorporate changes.

To be continued…

My point was that eltype should predict the return type of getindex(container, index::Int) (when that method is defined):

help?> eltype
search: eltype fieldtype fieldtypes

  eltype(type)

  Determine the type of the elements generated by iterating a collection of the given type....