Loading ICC profiles.

Profile usage.

Use of a profile is quite simple Let’s see the following code:

std::unique_ptr<profile_t> handler(profile_t::create("profile.icc"));
function_t<rgb_t, xyz_t> f = handler->pcs2device<rgb_t, xyz_t>();


Here we create a profile object and then we ask for a function giving a conversion from the profile connection space to the device space. If the profile contains a suitable transform (which may not be the case if we used, by mistake a profile for a CYMK device) then we get it in the function f. What if the PCS is lab instead of xyz? Well, in this case a xyz->lab conversion will be added automatically to the profile and added to the ‘normal’ content.

How is all this achieved? Well, I will keep all the boring details of the file parsing, the construction of the transforms, and so on out of the way. All I will tell you is that the profile_t object (actually a derived class) holds an algo_base_t object indexed by its signature (so there is an algo_base_t object for the BToA0Tag, for example). Then, given a rendering intent, the correct object is retrieved and turned to the type requested. But this step is not trivial, a dynamic cast it’s not the solution. Rather we need to implement some kind of type reflection and then, if required, add the necessary conversions.

This requires some metaprogramming, visitors and other stuff not for the faint of heart. And this requires a dedicated chapter in which I will explain my use of the visitor pattern in this context.

Multi dimensional lookup table.

Motivation

A multi dimensional lookup table is present in nearly every ICC profile. But they can be adopted also in other situations, for example where the color conversion is expensive to compute and may be useful to turn a generic function_t<Y, X>, build from the composition of several steps, into a single N- dimensional table with precomputed values on grid and then make use of interpolation for intermediate values. This library has the components to make this task easy, can be applied to dimensions bigger than three, and can use both multi-linear and tetrahedral interpolation. After all this advertising let’s start to give a look at how this is implemented.

The core class

This is the man class, I just stripped off a the trivial parts like the cloning function. The eval(const X &x) const function first of all finds the cube/square/element where x stays. This is used to fill in the delta vector and in this way the funcion has paid its dependency on input space by turning all into a vector. Then the result is handed to the data_ object that holds the grid points to make the interpolation. This part is a weighted sum of the grid points, and so the output space has to implement the mutiplication by a scalar and addition. That’s straightforward for floating point types, requires an extra trick when one has to deal with unsigned char

template 
class lut_t : public algo_t
{
public:
    enum {dimension = X::dimension};
    lut_t(const size_t *size, const X &step) : 
          data_(size), step_(step) {}
    virtual Y eval(const X &x) const override
    {
        delta_item_t deltas[dimension];
        // find the cube x belongs (and the remainders)
        input_traits_t::split(x, 
                              step_, 
                              deltas, 
                              data_.size1());
        // build the interpolation
        return data_.combine(deltas, I());
    }
 private:
    arrayn_t  data_; // TODO: consider
    X                  step_; // sharing this data..
};

Grid for 2D case

But now it’s time to see how the split is done. To make things simple to follow I will describe the 2D case, but the reasoning can be extended to any dimension. On the plane the multi-dimensional lookup table is a unit square with each side divided in N parts. The function has been evaluated at each corner of the grid. Given a generic point P(x,y) the first thing to do is to find out to which sub-square it falls in. That’s quit easy: just divide each coordinate by the side of the sub-square. The integer part will be the index (xi) for that dimension and the remainder (dx) will be related to the weights to interpolate the table values next to it.

Before carrying on let me stress again that this operation (finding the indexes and remainders) is dependent on the type of the input space. Is a single concept that needs to be implemented in a way if the input space is an N-dimensional vector, and in another if input space is e.g. rgb. The library defines a trait class for this operation, and provides (at the moment) two implementation: one for rgb and one for the vector template. They have been crafted carefully, so they work both for floating point types and for integer types like rgb_t<unsigned char>.

This is the code for the rgb_t case. The scalar_traits::one() function returns 1.0 for floating point types, and returns 255 for unsigned char.

template 
struct lut_input_traits_t<rgb_t>
{
    typedef T scalar_type;
    typedef typename scalar_traits_t::weight_type weight_type;
    static void split(const rgb_t &x,
                      const rgb_t &step,
                      delta_item_t *deltas,
                      const size_t *size1)
    {
        split_0(x.red, step.red, deltas, size1);
        split_0(x.green, step.green, deltas + 1, size1 + 1);
        split_0(x.blue, step.blue, deltas + 2, size1 + 2);
    }
    static void split_0(T x, T step, 
                        delta_item_t *deltas, 
                        const size_t *size1)
    {
        T div = x / step;
        deltas->offset = static_cast(div);
        deltas->offset -= deltas->offset / size1[0];
        // deltas->value = x - deltas->offset * step;
        deltas->value = x * scalar_traits_t::one() / step - 
                deltas->offset * scalar_traits_t::one();
    }

};

Triangle
At this point if we want to do a multi-linear interpolation is all quite ready. But if we want to do tetrahedral interpolation (triangular in 2D) then there is another problem: we must find out in which triangle our point falls in. In 2D this is quite simple. If remainder (dx) on x axis is bigger than dy then the point is in the lower triangle (like in the picture on the right), otherwise in the upper. This means that to find the point involved in the interpolation we can do the following: starting from the bottom right corner we increment first the dimension that has the bigger remainder. This is a reasoning that can be easily generalized: if we have a space of dimension N then we take the remainders (dx1, dx2, …dxn) and sort them in decreasing order. I do not have a formal proof of this fact, just the fact that… works (and gives non-negative weights for the values at the vertexes, which is definitely a must).

Below there is the code for this operation. The perm_ object does the sorting. This is an actual sort only for N > 5, otherwise there is another object that evaluates dx[i] < dy[j] for every i < j and builds a sort of 'signature' of the input. For N = 3 this means just 3 comparisons and, more important, this means that the whole operation can be don without a single if instruction. So the processor pipeline of the processor never gets reset and computation is done at the highest possible speed. For values bigger then N the lookup table required starts to grow too much, and I think that a standard sort is better.

template <class Y, int N>
class arrayn_t
{
public:
    template <class S> // S is a scalar
    Y combine(delta_item_t<S> *deltas, interp_tetra_t) const
    {
        // computes offset of first cube ans substitutes
        // the indexes with offset for each dimension
        const Y *data = data_ + 
            iteration_t<S, N>::substitute(deltas, offset_);
        // sort deltas to have them in decreasing order
        perm_.indexes(deltas);
        // for byte type one() is 256, hence the trait
        return data[0] * 
               (scalar_traits_t<S>::one() - deltas->value) +
                voxel_sum_t<Y, S, N>::tetrahedral(data, deltas);
    }
private:
    size_t                    size_[N];
    size_t                    size1_[N]; 
    size_t                    offset_[N + 1];
    Y                        *data_;
    permutation_t<N, (N > 5)> perm_;
};

Now we can compute the weighted sum of the points at the vertexes. This means that the type Y must have defined an operator* (multiplication by scalar) and an addition. Then there is a trick for integer types. The trick is that the operator * for rgb_t<unsigned char> gives actually a rgb_t<int> that is later on down-scaled again to rgb_t<unsigned char>. Yes, is a bit tricky, but in this way I have a single implementation for the interpolation, which I consider a Good Thing (TM). Let me know your opinion about that. If you want to see the implementation of this part then look at the soure code. voxel_sum_t is a simple template that unrolls the loop over N.

It has been hard, but at last we have enough tools to read a profile and build a transform. But this will be matter for another post…

First building blocks: functions.

The cornerstone at the building of the basement is a template function_t<Y, X>, where X is the space (type) of the domain and Y is the space (type) of the range of the function. This template defines a member function Y eval(const X &) that does the colour transform. This class has polymorphic behavior, but having also value semantic would be handy.

The solution adopted to have both benefits is an envelope-letter idiom (see Coplien’s Advanced C++ Programming Styles and Idioms). The role of the letter is played by another template class algo_t<Y, X>.

But it’s time to show some code: here’s the outline of the function_t template:


template <class Y, class X>
class function_t
{
public:
    typedef X domain_type;
    typedef Y range_type;
    // Construction from an algorithm
    function_t(algo_t<Y, X> *algo) : imp_(algo) {}
    // Copy construction. Clones the algorithm
    function_t(const function_t<Y, X> &rhs) : imp_(rhs.clone()) {}
    operator bool() { return imp_.get() != nullptr; }
    // Operations
    range_type eval(const domain_type &x) const ///< Function evaluation
    {
        return imp_->eval(x);
    }
    range_type operator()(const domain_type &x) const
    {
        return eval(x);
    }
    algo_t<Y, X> *clone(void) const
    {
        return imp_->clone();
    }
    template <class Z> ///< Function composition
    function_t<Y, Z> operator*(const function_t<X, Z> &rhs) const operation
    {
        return function_t<Y, Z>(new composite_t<Z>(imp_->clone(), 
                                                      rhs.clone()));
    }
private:
    std::unique_ptr<algo_t<Y, X>> imp_;
};

Function composition

Next to evaluation the function_t<Y, X> class defines another operation: functional composition. function_t<Y, X> can be composed with a function_t<X, Z> giving a function<Y, Z> The algorithm that does the composition is a class nested inside function_t (so in the following code X and Y are the same of the previous section).


template <class Z>
class composite_t : public algo_t<Y, Z>
{
public:
    composite_t(algo_t<Y, X> *lhs,
              algo_t<X, Z> *rhs) : lhs_(lhs), rhs_(rhs) {}
    // evaluation of functional composition
    virtual Y eval(const Z &z) const override
    {
        X x = rhs_->eval(z);
        return lhs_->eval(x);
    }
    virtual composite_t<Z> *clone(void) const override
    {
        return new composite_t<Z>(lhs_->clone(), rhs_->clone());
    }
private:
    std::unique_ptr<algo_t<Y, X>> lhs_;
    std::unique_ptr<algo_t<X, Z>> rhs_;
};

A few standard colour spaces (types) are defined: Lab, XYZ, RGB. RGB is actually a template, since a single channel can be implemented with a floating point type (with values ranging from 0.0 to 1.0) or with an integer type (usually unsigned char).

Additional features

The algo_t<Y, X> template is derived from two base classes, algo_domain_t<X> and algo_range_t<Y> that are defined to allow the use of the visitor pattern. This will be used to do some type inspection on a non template algo_base_t class. Later on I will explain how this feature is used.

First actual algorithms.

Standard color conversion algorithms are defined to convert between spaces. There is a single template named color_conversion_t<Y, X> provided for this task. Actually the generic implementation (see code below) does nothing, just defines a type used to state the fact that the conversion is not available.


struct domain_null_t{};
template <class Y, class X>
class color_conversion_t : public algo_t<Y, X>
{
public:
    typedef domain_null_t domain_t;
};

Then there is a plethora of specializations of this template for actual conversions: from algorithmic one like xyz->lab, to other that does just a scaling, like rgb_t<double> -> rgb_t<unsigned char>. To have just one template for all conversions will be handy later on, when the library will be able to insert a conversion automatically…but this is another story, which will be told in another post 🙂

Other algorithms.

To implement the functions defined in an ICC profile the library defines a 1D lookup-table, a matrix transform and a multi-dimensional lookup table. But this is a rather complex component, which deserves a post on its own…

Motivation – Abstract – whatever

The raison d’être of this blog is to discuss about libicc++ library. The aim of this library is provide a simple and flexible way do deal with ICC color profiles and color transforms in general. There are already free libraries in this field, the most notable being http://www.littlecms.com/. I have something different in mind, a bit difficult to express, so I started this project and this blog, to explain what is in it.

Here I want discuss about this library, show its design guidelines, and look for hints, points of view, comments, and whatever can be useful to improve it. The source code repository is available at https://github.com/marco-o/LIBICCPP

While reading keep in mind that I could, in a peak of self-estimation, claim to be a professional programmer, but surely not a writer. Furthermore I am not a native English speaker so, well, be patient with the approximation of the language and all the rest…

The net result will be discussions about C++ coding styles, idioms, patterns and color management algorithms. You will see a seasoned developer finding its way into C++11. Enojy!