Volumes concept philosophy

Volumes rely on the allocator system of somaio: see Allocators system.

First, a Volume is built using any allocator method. This means a DataSource and an AllocatorContext are provided to it. It can remain unallocated (no data buffer) if it is not used directly, but its data values should not be accessed until it is allocated.

Volumes and sub-volumes

A subvolume can be opened into an existing volume. This used to be a separate class, VolumeView, up to aims-4.3, but has been merged with the regular Volume class in aims-4.4. A view is used either to operate on smaller data blocks than the whole volume itself, or to simulate the "volume with borders" concept, or to "paste" a volume (buffer or file) into a larger one. Access to views is granted from one (view) volume into another one, using the Volume::refVolume(), and Volume::posInRefVolume() methods. Depending on the AllocatorContext and DataSource used for both the bigger volume and the view, different behaviours can be achieved. In the following we name "Volume" the "larger" volume, and "View" the view in the former.

• Volume allocated and linked to an existing DataSource (file or filled buffer), and View with no DataSource:
  This is the classical subvolume (region) case. All data for the whole volume are allocated (possibly using memory mapping), and the View is only useful to work on a restriced portion of the data, or to write / copy the View to another file or buffer.
  The AllocationContext of the View is not used in this case, the data block always belongs to the bigger Volume.
• Volume allocated with no DataSource, and View linked to an existing DataSource:
  This case is the "volume with borders" case: the buffer is allocated with borders included in the bigger Volume, and data are loaded directly in the View window. Work can be done in the View with no care of going in the borders (when borders are large enough).
  Another use of this is loading blocks of different files into a bigger volume without allocating nor copying a second volume.
  The AllocationContext of the View is not used in this case, the data block always belongs to the bigger Volume.
• Volume allocated and linked to an existing a DataSource, and View also linked to an existing DataSource:
  In this case the DataSource of the View will be used to load data into it and replace data from the original Volume (the Volume must be allocated in read/write mode). This case only differs from the previous one (borders) in that the "borders" are filled by the contents of another data source.
• Volume unallocated but linked to an existing DataSource, and View allocated with no DataSource:
  This is used to allocate and load only a part of a file: only the view is allocated (using its own AllocatorContext) and loaded from the bigger Volume Datasource: offsets in a large file correspond to the bigger Volume. If allocated read/write with memory mapping, changes will also affect the file. Otherwise, saving the whole modified volume is not possible at present.
  We may think about using the same trick for writing: writing from the View, but using the output DataSource of the Volume, would save only the viewed part, or copy the original file outside of it. But this will not be implemented at the beginning.
  

  Warning

  As buffers are actually held in the View, several views into the same unallocated Volume will not share data: this means that overlaps will be duplicated and independent in each view. This is not the case when using views on an allocated Volume.

• Volume unallocated with no DataSource:
  Whatever the SubVolume in it, I think this configuration is useless. The View will act as a full Volume, the bigger one won't be used.

Here I am speaking of Volume and View implemented as two separate classes, but they might be implemented in the same class: a Volume would be a view into another if it holds a reference to another (bigger) Volume.

In any case, a View is a complete Volume and has all functionalities of a Volume. This includes the capability to allow other views into it.

A View holds a reference-counting pointer to the Volume it is open on. So the original bigger Volume cannot be deleted while a View looking into it is still living.

Using cleverly the different combinations of allocators and DataSource on Volumes and Views allows to implement every memory-saving operations on very large files: memory mapping, block-by-block sequential reading, etc.

Attention

What will a View become if the Volume it views in resizes or is reallocated ?

Iterating over volumes

There are several ways to iterate over the voxels of a volume:

• Accessors: volume.at( position ). If the number of dimensions is fixed (typically 3, or 4) then the programmer may write the required number of nested loops:

  Volume<float> vol( 10, 10, 10 ); // fill volume
  std::vector<int> size = vol.getSize();
  float sum = 0.;
  for( int z=0; z<size[2]; ++z )
    for( int y=0; y<size[1]; ++y )
      for( int x=0; x<size[0]; ++x )
        sum += vol.at( x, y, z );

• iterators: the iterator returned by volume.begin() will iterate over all dimensions. Depending on the implementation chosen at compilation time (blitz++ support or not), the behavior will be different: with blitz++ support, iterators will iterate as expected over each "actual" voxel of the volume, but with reduced efficiency. Witout blitz++ support, iterators will act like pointers: they will not take into account borders in the case of a view into a larger volume. Thus it will make too many iterations, and thus be wrong, in this case.

  Volume<float> vol( 10, 10, 10 ); // fill volume
  Volume<float>::iterator it, e = vol.end();
  float sum = 0.;
  for( it=vol.begin(); it!=e; ++it )
    sum += *it;

• N-D iterators: if the volume contains an arbitrary number of dimensions, then the nested loops cannot be hard-coded in the iterating code. The NDIterator and line_NDIterator classes in the CartoBase library provide convenient ways of walking through a Volume:

  Volume<float> vol( 10, 10, 10 ); // fill volume
  std::vector<size_t> sstrides = vol.getStrides();
  std::vector<int> strides;
  strides.insert( strides.end(), sstrides.begin(), sstrides.end() );
  float sum = 0.;
  NDIterator<float> it( &vol.at( 0 ), vol.getSize(), strides );
  for( ; !it.ended(); ++it )
    sum += *it;

  or, using line_NDIterator:

  Volume<float> vol( 10, 10, 10 ); // fill volume
  std::vector<size_t> sstrides = vol.getStrides();
  std::vector<int> strides;
  strides.insert( strides.end(), sstrides.begin(), sstrides.end() );
  float *p, *pp;
  float sum = 0.;
  line_NDIterator<float> it( &vol.at( 0 ), vol.getSize(), strides );
  for( ; !it.ended(); ++it )
  {
    p = &*it;
    for( pp=p + vol.getSizeX(); p!=pp; ++p )
      ++(*p);
  }

In terms of performance, the impact depends on the "weight" of the operation performed with each voxel, the lighter this operation, the more impacted it is by the iteration overhead. For a very simple operation (like incrementing each voxel value), if we take the Volume::iterator performance as reference 1., we can expect approximately:

Iteration type	Perf.	N-D support	view support
Volume::iterator with blitz++ support	1.	yes	yes
pointers	10.	yes	no
accessors with blitz++ support	5-10	no	yes
NDIterator	1.	yes	yes
line_NDIterator	10.	yes	yes

Thus in most cases the best compromise between performance and flexibility (support for the general N-dimension case, support for views inside larger volumes) is using line_NDIterator.

For comparison, here are some tests results on a Core i7-4600U 2.1GHz CPU (laptop) using different compilers on linux and trying different types of implementations, all running on the same machine. The reference (factor 1.) was actually 4.2e+08 voxels/s using iterators with gcc 4.9. The implementation choice was taken as the best between implementations (marked with a star), and depends on the compiler. It can be seen that there are striking behavior changes between compilers. Moreover it seems that there are also different behaviors using different machines / processors, so the choices we have made might not be the best ones for every machine. However for now they are hard-coded using preprocessor #if directives.

Iteration type	gcc 4.9	gcc 6	clang 3.8	clang 6
iterator	1.	1.	0.9	0.8
pointers	10.9	10.9	10.6	11.2
accessors (raw 4D)	4.6 *	4.6 *	5.4	5.4
accessors (blitz)	4.5	4.6	10.7 *	11. *
vector accessor (loop)	0.2	2.1 *	3.8 *	5.4 *
vector accessor (switch on dim)	3.5 *	1.2	0.5	0.5
vector accessor (blitz)	0.3	0.7	0.17	0.2
NDIterator	0.9	1.2	1.2	1.2
line_NDIterator	10.7	10.4	10.6	11.