Re: [boost] Boost.Numeric: multi_array requirements

From: Ullrich Koethe (u.koethe_at_[hidden])
Date: 2001-04-02 11:56:38

Thanks for answering. I'd really like it very much if the array library
would be as good for image processing as it is for linear algabra.

Details below.

Regards
Ulli

Jeremy Siek wrote:
> 
> On Thu, 29 Mar 2001, Ullrich Koethe wrote:
> koethe> Indexing:
> koethe> ---------
> koethe>
> koethe> Most proposals I've seen so far are using integers for indexing.
> koethe> However, in higher dimensions this causes all kinds of trouble - Shall
> koethe> we use (i,j,k) or [i][j][k]? Is the horizontal index first or last? Is
> koethe> the first index 0 or 1? Therefore, I've found *index objects* a very
> koethe> powerful alternative. For example:
> koethe>
> koethe> Index2D start = matrix.upperLeft(),
> koethe>         end   = matrix.lowerRight().
> koethe>         i;
> koethe>
> koethe> for(i.y = start.y; i.y < end.y; ++i.y)
> koethe> {
> koethe>     for(i.x = start.x; i.x < end.x; ++i.x)
> koethe>     {
> koethe>         matrix[i] = 42;
> koethe>     }
> koethe> }
> 
> Actually, the multi_array spec included an elt() function which takes the
> equivalent of one of these index objects. There are two small differences:
> 
>   - the index object is just sequence of indices (not data members
>      x and y)
>   - the elt() function took an iterator to the beginning of the sequence.
> 

This could  easily be made into an object

> 
> So, your example above would have looked like this:
> 
> for (i[0] = start[0]; i[0] < end[0]; ++i[0])
>   for (i[1] = start[1]; i[1] < end[1]; ++i[1])
>     array.elt(i.begin());
> 
> I think you'll agree that the "x" and "y" are very must image specific,
> and it would be better to go with a sequence of indices.
> 

Yes and no. One important reason I've choosen "x" and "y" is that there
never is any confusion as to what is the horizontal and what is the
vertical index. If you use "i" and "j" or "0" and "1", it is never clear
which is which. (In your example above, I would have intuitively choosen
"1" for the outer, not the inner loop.) I've found it very important to
eliminate this error source altogether, because these errors are
otherwise extremely hard to find.

Also, I've experimented with different syntax variants:

++i.x
++i.x()
++i[x]
++i(x)

I found the first most readable, but the third is certainly a good
alternative - but it requires a global enum or object "x" which easily
interferes with other x'es. Ideally, all variants would be aliases for
the same thing, but this doesn't seem possible with C++ (Efficiently,
that is). I'd really appreaciate a good suggestion here.

> So what's left to debate:
>   - the name of the function. I'm not all that crazy about elt().
>     My current favorate is operator() because it is not operator[]
>     (which people are use to using with a single integer) and
>     because it is shorter.
> 

Using operator[] is useful because you can define functions
independently of array dimension. example (over-simplified):

template <class Array, class Index>
typename Array::value_type sqr(Array & a, Index i)
{
        return a[i]*a[i]; // works if i is int or an index object
}

> 
>     Either way, it be cool to make the element access funtion
>     templated, so people could make their own index objects.
> 

Agreed.

> koethe> Iterators:
> koethe> ----------
> koethe>
> koethe> A property of images is that all coordinate directions are 'equal'.
> koethe> Therefore, a design like the following, where one direction is the
> koethe> 'outer' direction, and the other the 'inner', is very artificial:
> koethe>
> koethe> JeremySiek> 1. Use two iterators, the "outer iterator" iterates over
> koethe> rows
> koethe> JeremySiek> of the image, and the "inner iterator" iterates  down each
> koethe> row.
> koethe> JeremySiek>
> koethe> JeremySiek> PixelRowIterator row_iter = rows_begin(myPix);
> koethe> JeremySiek> for (; row_iter != rows_end(myPix); ++row_iter) {
> koethe> JeremySiek>   PixelIterator iter = begin(*row_iter);
> koethe> JeremySiek>   for (; iter != end(*row_iter); ++iter)
> koethe> JeremySiek>     ...
> koethe> JeremySiek> }
> koethe>
> koethe> In VIGRA, I'm preferring a different iterator concept that doesn't have
> koethe> this asymmetry (expet for cache performance):
> koethe>
> koethe> ImageIterator start = image.upperLeft(),
> koethe>               end = image.lowerRight(),
> koethe>               i;
> koethe>
> koethe> for(i.y = start.y; i.y < end.y; ++i.y)
> koethe> {
> koethe>     for(i.x = start.x; i.x < end.x; ++i.x)
> koethe>     {
> koethe>        *i = 42;
> koethe>     }
> koethe> }
> koethe>
> 
> This isn't really a conceptual difference, just a difference in packaging.
> It's interesting, but it produces iterators that are not std conformant.
> If I see iterators, I want to be able to use them with std::copy(), etc.

I'm not sure about this. After all, a 2D array is conceptually very
different from a sequence. This is reflected in the 2D iterators. In
addition to them, I've lots of iterator adapters that realize different
mappings from 2D to STL-conformant 1D sequences (e.g. ScanOrderIterator,
RowIterator, ColumnIterator, ContourFollower)

> Also, there might be some performance advantage in the "hierarchical"
> approach because the index*extent multiplications can be moved to outer
> loops, effectivly performing what compiler guys call strength reduction.
> With the above iterators, you might as well just use the index objects
> instead.
> 

I've explicitly specified in VIGRA that it is faster to use x in the
inner loop. But some algorithms must be written the other way around
(the simplest being image transposition), or with no directional
preference whatsoever (e.g. contour following).

> koethe> Different views on the pixel data:
> koethe> ----------------------------------
> koethe>
> koethe> This has proven more important than I expected. Consider, for example,
> koethe> an image that contains complex numbers (e.g. the Fourier transform of an
> koethe> image). Now, the magnitude of a complex number is a scalar that can be
> koethe> used in any algorithm that requires a scalar image (for example, image
> koethe> display). By using *data accessors* to create a magnitude view, I can
> koethe> calculate the magnitude on the fly, without creating an extra magnitude
> koethe> image.
> koethe>
> koethe> displayImage(complex_image, MagnitudeAccessor());
> 
> Yes, this is a very important capability.
> 

Which would be very much simplified if "x[i] = a" weren't a requirement.
How about this:

x.set(i, a) is required.
x[i] = a is optional

This means that

1. If a data structure can implement "x[i] = a" efficiently, it shall do
so
2. A truly generic algorithm should only use x.set()
3. If you know that all your data structures implement "x[i] = a", you
may use it everywhere.

-- 
 ________________________________________________________________
|                                                                |
| Ullrich Koethe  Universität Hamburg / University of Hamburg    |
|                 FB Informatik / Dept. of Computer Science      |
|                 AB Kognitive Systeme / Cognitive Systems Group |
|                                                                |
| Phone: [phone]                Vogt-Koelln-Str. 30 |
| Fax:   [phone]                D - 22527 Hamburg   |
| Email: u.koethe_at_[hidden]               Germany             |
|        koethe_at_[hidden]                        |
| WWW:   http://kogs-www.informatik.uni-hamburg.de/~koethe/      |
|________________________________________________________________|