Writing your own Functions

Sooner or later, you will want to implement your own functions on the basis of VIGRA's functionality. Some people believe that this is very difficult because one needs to provide a lot of template magic and full genericity. However, this is not true: Your VIGRA functions need not be templated at all – function arguments can simply be hard-wired. In other cases, it makes sense to template on the pixel type, but leave averything else fixed. Full genericity should only be implemented step-by-step as needed.

As an example, consider again the image smoothing example program smooth_explicitly.cxx. It makes sense to encapsulate the smoothing algorithm into a function of its own. When we only need to support float images, the function is simply a verbatim copy of the algorithm. In contrast to the original version, we now allow an arbitrary window radius to be passed to the algorithm, so that the amount of smoothing can be controlled (this also nicely illustrates the use of vigra_precondition() for Error Reporting):

void smooth(MultiArrayView<2, float> input, MultiArrayView<2, float> result, int radius)
{
    vigra_precondition(radius >= 0, "smooth(): window radius must not be negative.");

    Shape2 current;
    for(current[1] = 0; current[1] < input.shape(1); ++current[1])
    {
        for (current[0] = 0; current[0] < input.shape(0); ++current[0])
        {
            Shape2 windowStart = max(Shape2(0),     current - Shape2(radius));
            Shape2 windowStop  = min(input.shape(), current + Shape2(radius+1));
            MultiArrayView<2, float> window = input.subarray(windowStart, windowStop);
            result[current] = window.sum<float>() / window.size();
        }
    }
}

If we don't need to support any higher dimension or other pixel type, we can just leave it at this – no templates are required then.

But suppose now that we want to generalize this code for arbitrary dimensional arrays. To do so, we specify the dimension N as a template parameter. Then we can no longer use Shape2 because this class only works for 2-dimensional arrays. Instead, we use the MultiArrayShape traits class to ask for the appropriate shape object. Moreover, we cannot iterate over the array with two explicitly nested loops because the number of loops must correspond to the (unknown) number of dimensions. We can solve this problems by means of a vigra::MultiCoordinateIterator from multi_iterator_coupled.hxx that iterates over all coordinates of an array, regardless of dimension. The current coordinate is returned by dereferencing the iterator:

#include <vigra/multi_iterator_coupled.hxx>

template <unsigned int N>
void smooth(MultiArrayView<N, float> input, MultiArrayView<N, float> result, int radius)
{
    vigra_precondition(radius >= 0, "smooth(): window radius must not be negative.");

    typedef typename MultiArrayShape<N>::type  Shape;

    typename MultiCoordinateIterator<N> current(input.shape()),
                                        end = current.getEndIterator();

    for(; current != end; ++current)
    {
        Shape windowStart = max(Shape(0),      *current - Shape(radius));
        Shape windowStop  = min(input.shape(), *current + Shape(radius+1));
        MultiArrayView<N, float> window = input.subarray(windowStart, windowStop);
        result[*current] = window.sum<float>() / window.size();
    }
}

Another useful generalization is in terms of the array's value_type. For one, we want to be able to smooth color images as well. Furthermore, most images are stored with pixel type unsigned char, and we don't want to force the user to convert them into float images before smoothing. We therefore specify the value_types as template parameters as well (notice that we allow input and result to have different types). In addition, we have to make the type of the sum in window.sum<...>() generic. However, there is a caveat: We cannot simply use the input value_type here, because this might lead to overflow. This is easily seen when the value_type is unsigned char: This type already overflows when the sum exceeds the value 255, which is very likely to happen even if the windows is only 3x3. In situations like this, a suitable temporary type for the sum can be obtained from the RealPromote type in VIGRA's NumericTraits class:

template <unsigned int N, class InputValue, class ResultValue>
void smooth(MultiArrayView<N, InputValue>  input,
            MultiArrayView<N, ResultValue> result,
            int radius)
{
    vigra_precondition(radius >= 0, "smooth(): window radius must not be negative.");

    typedef typename MultiArrayShape<N>::type                Shape;
    typedef typename NumericTraits<InputValue>::RealPromote  SumType;

    typename MultiCoordinateIterator<N> current(input.shape()),
                                        end = current.getEndIterator();

    for(; current != end; ++current)
    {
        Shape windowStart = max(Shape(0),      *current - Shape(radius));
        Shape windowStop  = min(input.shape(), *current + Shape(radius+1));
        MultiArrayView<N, InputValue> window = input.subarray(windowStart, windowStop);
        result[*current] = window.template sum<SumType>() / window.size();
    }
}

These simple tricks already get you a long way in the advanced use of VIGRA. You will notice, that many existing VIGRA functions are not implemented in temrs of vigra::MultiArrayView, but in terms of image iterators and hierarchical iterators. However, these iterators are more difficult to use, so the MultiArrayView approach is recommended for new code.