File: tutorial.dxx

package info (click to toggle)
libvigraimpex 1.10.0%2Bgit20160211.167be93%2Bdfsg-2
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 48,024 kB
  • ctags: 25,337
  • sloc: cpp: 55,572; python: 8,495; ansic: 1,267; makefile: 150; sh: 38
file content (975 lines) | stat: -rw-r--r-- 56,966 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
/** \page Tutorial Tutorial
    
    This tutorial will help you to learn VIGRA's most important concepts by means of simple examples. The tutorial consists of the following parts:
    
    <ul style="list-style-image:url(documents/bullet.gif)">
    <li> \ref MultiDimensionalArrayTutorial
        <BR>&nbsp;&nbsp;&nbsp;<em>VIGRA's most important data structure</em>
        <ul style="list-style-image:url(documents/bullet.gif)">
        <li> \ref MultiArrayBasics
        <li> \ref MultiArrayIndexing
        <li> \ref MultiArrayScanOrder
        <li> \ref MultiArrayMethods
            <ul type="disc">
            <li> \ref MultiArrayViewBasics
            <li> \ref MultiArray_subarray
            <li> \ref MultiArray_bind
            <li> \ref MultiArray_vector_elements
            <li> \ref MultiArray_transpose
            <li> \ref MultiArray_unstrided
            </ul>
        </ul>
        
    <li> \ref MultiArrayArithmeticTutorial
        <BR>&nbsp;&nbsp;&nbsp;<em>mathematical operations on MultiArrays</em>
        <ul style="list-style-image:url(documents/bullet.gif)">
        <li> \ref MultiMathModule "Array Expressions"
        <li> \ref LinearAlgebraModule "Linear Algebra"
        <li> \ref MultiPointoperators "STL-style transformation algorithms"
        <li> \ref FeatureAccumulators
        </ul>
        
    <li> \ref ImageInputOutputTutorial
        <BR>&nbsp;&nbsp;&nbsp;<em>importing and exporting images and arbitrary-dimensional arrays</em>
        <ul style="list-style-image:url(documents/bullet.gif)">
        <li> \ref Impex2D
        <li> \ref ImpexND
        </ul>

    <li> \ref ImageProcessingTutorial
        <BR>&nbsp;&nbsp;&nbsp;<em>basic applications of VIGRA's functions</em>
        <ul style="list-style-image:url(documents/bullet.gif)">
        <li> \ref CallingConventions
        <li> \ref ImageInversion
        <li> \ref ImageBlending
        <li> \ref CompositeImage
        <li> \ref SmoothingTutorial
           <ul type="disc">
            <li> \ref Convolve2DTutorial
            <li> \ref SeparableConvolveTutorial
            </ul>
        </ul>

    <li> \ref OwnFunctionsTutorial
        <BR>&nbsp;&nbsp;&nbsp;<em>... without getting confused by templates</em>

    <li> \ref PythonBindingsTutorial
        <BR>&nbsp;&nbsp;&nbsp;<em>scripting with VIGRA made easy</em>
        
    </ul>
*/

/** \page MultiDimensionalArrayTutorial Multi-Dimensional Arrays

    <h2>Section Contents</h2>
    
    <ul style="list-style-image:url(documents/bullet.gif)">
    <li> \ref MultiArrayBasics
    <li> \ref MultiArrayIndexing
    <li> \ref MultiArrayScanOrder
    <li> \ref MultiArrayMethods
        <ul type="disc">
        <li> \ref MultiArrayViewBasics
        <li> \ref MultiArray_subarray
        <li> \ref MultiArray_bind
        <li> \ref MultiArray_vector_elements
        <li> \ref MultiArray_transpose
        <li> \ref MultiArray_unstrided
        </ul>
    </ul>
    
    \section MultiArrayBasics Basic MultiArray Usage
    
    \ref vigra::MultiArray is the most fundamental data structure in VIGRA. It holds a rectangular block of values in arbitrary many dimensions. Most VIGRA functions operate on top of MultiArray or the associated class MultiArrayView (see \ref MultiArrayViewBasics). 
    
    A 2D image can be interpreted as a matrix, i.e. a 2-dimensional array, where each element holds the information of a specific pixel. Internally, the data are stored in a single 1-dimensional piece of memory, and MultiArray encapsulates the entire mapping between our familiar 2-dimensional notation and the raw memory. Pixels in an image are identified by a coordinate pair (x,y), where indexing starts at 0. That is, the pixels in a 800x600 image are indexed by <tt>x = 0,...,799 and y = 0,...,599</tt>. The principle analoguously extends to higher dimensions.
    
    The structure of a multidimensional array is given by its <tt>shape</tt> vector, and the length of the shape vector is the array's <i>dimension</i>. The dimension must be fixed as a template parameter at compule time, while the shape is passed to the array's constructor. The second important template parameter is the pixel's <tt>value_type</tt>, as you know it form <tt>std::vector</tt>. 
    
    To represent the data of a gray scale image, we just need to store one value per pixel, so we 
    choose a 2-dimensional array, where each element has the <tt> unsigned char </tt> type 
    (in VIGRA, this type is also available as \ref vigra::UInt8). We instantiate a gray scale image object like this:

    \code
    #include <vigra/multi_array.hxx>
    
    using namespace vigra; // for brevity in the examples - don't do this in header files!
    
    int width = ...,  height = ...;
    MultiArray<2, UInt8> image(Shape2(width, height));
    \endcode
    
    By default, VIGRA arrays are <b>always zero-initialized</b>. Another initial value can be provided in the constructor, or later via the <tt>init()</tt> function or the assignment operator:
    
    \code
    MultiArray<2, UInt8> image(Shape2(width, height), 255); // init with value 255
    
    image.init(128);  // same effect, different initial value
    image = 100;      // yet another way
    \endcode
     
    The <tt>Shape2</tt> typedef also exists for higher dimensions up to five as <tt>Shape3</tt> etc. If you need even more dimensions, use <tt>MultiArrayShape<N>::type</tt> instead, were N is the number of dimensions: 

    \code
    // dimension 0 runs from 0, 1, ..., 299
    // dimension 1 runs from 0, 1, ..., 199
    // dimension 2 runs from 0, 1, ...,  99
    MultiArray<3, double> volume(Shape3(300, 200, 100));
    
    MultiArray<7, float> array7D(MultiArrayShape<7>::type(20, 10, ...));
    \endcode     

    When storing RGB images we obviously can't simply use the unsigned char type because every pixel contains 3 numbers: values for red, green and blue. Mathematically, you want to store a data vector for each pixel. To this end, VIGRA
    provides the <tt>vigra::RGBValue<ValueType></tt> class. So for RGB-images just use: </p>

    \code
    MultiArray<2, RGBValue<UInt8> > rgb_image(Shape2(256, 128));
    \endcode

    vigra::RGBValue<ValueType> is a specialised 3-dimensional vector containing ValueType elements. Arbitrary short vectors can be stored in the <tt>TinyVector<ValueType, SIZE></tt> class, which is the base class of RGBValue. It's length must be specified at compile time in the template parameter <tt>SIZE</tt>. Vectors whose length is known at compile time are very useful for the compiler to produce highly optimized code. Therefore, <tt>Shape2</tt> and it's higher-dimensional cousins are instances of <tt>TinyVector</tt> as well.

    Alternatively you can use a 3-dimensional array <tt>vigra::MultiArray<3, unsigned
    char></tt> to represent a color image. The third dimension has size 3 and contains the
    information for the red, green and blue channel. 
    
    \code
    MultiArray<3, UInt8> rgb_array(Shape3(256, 128, 3));
    \endcode

    However, we are not going to use this form in the tutorial because <tt>vigra::RGBValue<ValueType></tt> provides many helpful methods that are not available when color is just another array dimension.

    \section MultiArrayIndexing Array Indexing via Coordinates

    The easiest way to access the values of an array is via the coordinates. A tuple of coordinates can again be specified by the appropriate shape object, which must be passed to the array's indexing operator:

    \code
    int width = 300, height = 200;
    MultiArray<2, int> image(Shape2(width, height));

    // set all elements to 3
    image.init(3);

    // print pixel at x=1 and y=2
    int x=1, y=2;
    std::cout << image[Shape2(x,y)] << std::endl;  // output: 3
    \endcode

    <B>Important Remark:</B> Notice that VIGRA follows the mathematical convention of the index order: dimension 0 corresponds to the x (horizontal) coordinate, dimension 1 to the y (vertical) coordinate, and so on. Accordingly,  dimension 0 is changing fastest in memory: when we increase x by one, we get to the next memory location. In matrix jargon, this is also known as <i>Fortran order</i>. Many image processing libraries (e.g. <a href="http://www.imagemagick.org/">Image Magick</a>, <a href="http://opencv.willowgarage.com/">OpenCV</a>, and <a href="http://qt-project.org/">Qt</a>) use the same convention. However, others like Matlab and numpy, use the reverse order (so called <i>C order</i>). Don't be confused!
    
    Internally, shape objects are implemented in terms of the powerful \ref vigra::TinyVector class. This means that shape objects support the usual mathematical operations like addition, multiplication and scalar products. Coordinate computations can thus be performed on entire coordinate objects at once - there is no need to manipulate the individual coordinates separately.
    
    Nonetheless, in some circumstances it is more convenient to provide coordinates individually rather than in a shape object. This is possible with round brackets (x,y):

    \code
    // access via individual coordinates
    image(1,2) = 22;
    \endcode
	
    This kind of indexing is supported for dimensions up to five, and only if the array's dimension is known (this is not always the case: in a generic function where the dimension is a template parameter, you must use shape objects).
	In combination with the method <tt>shape(n)</tt>, that returns the length of the n-th dimension, 
    we can use the coordinates to set the element of an entire row or column:

    \code
    // set all elements of first row to 13
    for (int i = 0; i < image.shape(1); i++)
        image(1,i) = 13;
    \endcode

    On first sight, individual coordinates may seem to be necessary for iterating over the image or parts of it. But the following example shows that the same effect can be achieved with a shape object that is allocated outside the loop:
    3rd column of a 8x4-matrix (initialized with 5) to 7.

    \code
    // instantiate shape object (zero intialized by default)
    Shape2 p;
    // bind x=2
    p[0] = 2;
    // iterator over row 2
    for(p[1]=0; p[1]<image.shape(1); ++p[1])   
        image[p] = 7;
    \endcode

    We will discuss more powerful methods to access certain parts of an array in section \ref MultiArrayMethods.

    \section MultiArrayScanOrder One-dimensional Indexing and Scan-Order Iterator
	
    Regardless of the array's dimension, it is always possible to access elements with 1-dimensional index, its <i>scan-order index</i>, via the normal indexing operator. For example, <tt>array[1]</tt> refers to the index of the second array element. Defining a scan order is often called <i>flattening</i> of an array, because a high-dimensional data structure is accessed like a 1-dimensional vector. Notive that scan-order access in VIGRA does not require the data to be copied. 
    
	VIGRA defines scan-order by considering the dimensions from front to back. Thus, items are accessed such that only the x coordinate is incremented, while y (and possibly further coordinates) are held fixed at 0. When x is exhausted, y is incremented by one and the iteration starts again at x=0. To control iteration, the function <tt>array.size()</tt> returns the total number of elements: 

    \code
    MultiArray<2, int> intArray(Shape2(3,2));

    for(int k=0; k<intArray.size(); ++k=
        intArray[k] = k+1;
    
    // the array now contains the values
    //
    //   1 2 3
    //   4 5 6
    \endcode
	
    Alternatively, scan-order access can be achieved with an STL-compatible iterator pair obtained by calling <tt>array.begin()</tt> and <tt>array.end()</tt>. Continuing with the example above, we can write:

	\code
    // declare an alias for the iterator type
    typedef MultiArray<2, int>::iterator Iter;

    // iterate over intArray and print the elements in scan order
    for (Iter i = intArray.begin(); i != intArray.end(); ++i)
        std::cout << *i << " ";
    std::cout << std::endl;
    
    // output: 1 2 3 4 5 6
	\endcode

    The iterator is implemented by class <tt>StridedScanOrderIterator</tt> which encapsulates all the bookkeeping necessary to get the elements in the correct order, even when the array was transposed (see below). 

	Scan-order access is useful to implement pointwise operations, e.g. the addition of two matrices. The following code adds two matrices and stores the result in the first one: 

	\code
    MultiArray<2, int> matrix1(Shape2(3,3)),
                       matrix2(Shape2(3,3));
    ... // fill in data

    // use indexing
    for (int i=0; i < matrix1.size(); ++i)
        matrix1[i] += matrix2[i];
    
    // use iterators
    for (Iter i = matrix1.begin(), j = matrix2.begin(); i != matrix1.end(); ++i, ++j)
        *i += *j;
	\endcode

	This is convenient because the actual high-dimensional shapes of the arrays are of no significance for point-wise operations as long as the shapes match. Be careful: the arrays themselves have no way of checking this condition. So thefollowing code using two transposed shapes is perfectly valid C++, but has probably not the intended effect:

	\code
    MultiArray<2, int> matrix3(Shape2(3,2)),
                       matrix4(Shape2(2,3));  // transposed shape of matrix3
    ... // fill in data

    for (int i=0; i < matrix3.size(); ++i)
        matrix3[i] += matrix4[i];  // works, but may not have the intended effect
	\endcode

    By the way: VIGRA provides the += operator (and is cousins) to write this more concisely, and this operator throws an exception if the shapes don't match:

    \code
    matrix1 += matrix2;          // works fine!
    matrix3 += matrix4;          // error: shape mismatch!
    \endcode

    For more information on mathematical operations on arrays see the \ref MultiMathModule "multi_math" module.

    As mentioned, VIGRA's scan order is similar to the NumPy-method <tt>array.flatten()</tt>. You use it, 
    to copy a multi-dimensional array into an one-dimensional array, or to access elements in flattened order. The only
    difference is that NumPy uses "C-order" , i.e. the rightmost dimension takes priority, whereas
    VIGRA uses Fortran-order, i.e. the leftmost dimension takes priority. A method like flatten can be implemented in VIGRA as:

    \code
    MultiArray<2, int> intArray(Shape2(3,2));

    for(int k=0; k<intArray.size(); ++k=
        intArray[k] = k+1;
    
    // create 1D-array of appropriate size
    std::vector<int> flatArray(intArray.size());
    
    // copy 2D-array into 1D-array using the STL
    std::copy(intArray.begin(), intArray.end(), flatArray.begin());
    
    // print 1D-array on console
    // (same output as printing from the StridedScanOrderIterator directly) 
    for (std::vector<int>::iterator i = flatArray.begin(); i != flatArray.end(); ++i)
         std::cout << *iter << " ";
    std::cout << std::endl;
    \endcode

    To show the difference between VIGRA and NumPy we'll add the NumPy output, i.e. the result when we had
    used C-order in the code above:

    \verbatim
    flatArray - index     0       1       2       3       4       5 
    -----------------------------------------------------------------
    VIGRA-output:         1       2       3       4       5       6
    intArray - index    [0,0]   [1,0]   [2,0]   [0,1]   [1,1]   [2,1]
    -----------------------------------------------------------------
    NumPy-output:         1       4       2       5       3       6
    intArray - index    [0,0]   [0,1]   [1,0]   [1,1]   [2,0]   [2,1]
    \endverbatim

    To change the axis priorities of the StridedScanOrderIterator, look at the transpose-function
    in the next section.
    
    \section MultiArrayMethods Important MultiArray Methods

    This part of the tutorial explains important methods of MultiArray. However, before we proceed, we need to introduce the class \ref vigra::MultiArrayView.
    
    \subsection MultiArrayViewBasics The MultiArrayView Interface
    
    A \ref vigra::MultiArrayView has the same interface as a MultiArray (with the exception of <tt>reshape()</tt> and <tt>swap()</tt>), but it doesn't own its data. Instead, it provides a <i>view</i> onto the data of some other array. In contrast, a MultiArray owns its data and is responsible for freeing it in the destructor. MultiArrays are automatically converted into MultiArrayViews when needed.
    
    The point of this distinction is that MultiArrayViews can be used to access and manipulate the same data in many different ways <i>without any need for creating copies</i>. For example, we can work with a 2-dimensional slice of a volume dataset (i.e. a lower dimensional part of a 3D array) without first copying the slice into a 2D image. This is possible whenever the desired view can be realized by just manipulating the internal <i>mapping</i> from indices and shapes to memory locations, and not the memory layout itself.
    
    This possibility -- which is similarly implemented in other packages like Matlab and numpy -- is a key ingredient for efficient computations with multi-dimensional arrays. Thus, most VIGRA functions actually receive MultiArrayViews to maximize flexibility. This section describes the most important ways to create new MultiArrayViews from an existing array or view. The complete documentation is available in the \ref vigra::MultiArrayView reference.

    <hr>

    \subsection MultiArray_subarray subarray(p,q)
    
    This method creates a rectangular subarray of your array between the points p and q, where p (the starting point of the subregion) is included, q (the ending point) is not. <tt>subarray</tt> does not change the dimension of the array (this is the task of the various <tt>bind</tt>-methods).
    
    To give an example, we create a 4x4 array that consitst of a checkerboard with 2x2 squares:

    \code
    MultiArray<2, float> array_4x4(Shape2(4,4));  // zero (black) initialized
    
    // paint the upper left 2x2 square white
    array_4x4.subarray(Shape2(0,0), Shape2(2,2)) = 1.0;
    
    // likewise for the lower right 2x2 square, but this time we 
    // store the array view explicitly for illustration
    MultiArrayView<2, int> lower_right_square = array_4x4.subarray(Shape2(2,2), Shape2(4,4));
    lower_right_square = 1.0;
    
    // contents of array_4x4 now:
    //    1 1 0 0
    //    1 1 0 0
    //    0 0 1 1
    //    0 0 1 1
    \endcode

    The positions p and q are specified with the familiar <tt>Shape</tt> objects. In this example we simply overwrite parts of the array. The following larger example uses <tt>subarray</tt> to output a half-sized subimage around the center of the original image: <a href="subimage_tutorial_8cxx-example.html">subimage_tutorial.cxx</a> <br />
    The relevant part of this code is shown here (the functions <tt>importImage</tt> and <tt>exportImage</tt> are described in section \ref ImageInputOutputTutorial):

    \dontinclude subimage_tutorial.cxx
    \skip // read image
    \until exportImage

    After reading the (here: gray scale) image data to an array we need to calculate the 
    coordinates of our subimage. In this case we want to cut out the middle part of the image. 
    Afterwards we write the subimage into a new array. Look at the result:

    <Table cellspacing = "10">
    <TR valign="bottom">
    <TD> \image html lenna_small.gif "input file" </TD>
    <TD> \image html lenna_sub.gif "subimage output file" </TD>
    </TR>
    </Table>

    <hr>

    \subsection MultiArray_bind bind<M>(i) and bindAt(M, i)
    
    These methods bind axis M to the index i and thus reduce the dimension of the array by one. The only difference between the two forms is that the axis to be bound must be known at compile time in the first form, whereas it can be specified at runtime in the second. 
    
    Binding is useful when we want to access and/or manipulate a particular row or column of an image, or a single slice of a volume. In principle, the same can also be achieved by explicit loops, but use of <tt>bind</tt> often leads to more elegant and more generic code. Consider the following code to initialize the third column of an image with the constant 5:

    \code
    // initialize 200x100 image
    MultiArray<2, int> array2d(Shape2(200,100)); // zero initialized

    // initialize column 2 with value 5 using a loop
    for(int y=0; y<array2d.shape(1); ++y)
        array2d(2, y) = 5;
        
    // the same using bind
    array2d.bind<0>(2) = 5;
    \endcode

    NumPy-Users are familiar with the bind mechanism as "slicing". The example above written in numpy syntax becomes:
    \verbatim
    array2d[2, :] = 5      // NumPy-equivalent of array2d.bind<0>(2) = 5
    \endverbatim

    You can also initialize a lower-dimensional array with the bind-method:

    \code
    // initialize new 1D array with 3rd column of a 2D array
    MultiArray<1, int> array1d = array2d.bind<0>(2);
    \endcode

    The array <tt>array1d</tt> contains the elements the 3rd column of <tt>array2d</tt>. This bahavior nicely illustrates the difference between a copy and a view: <tt>array1d</tt> contains a copy of the 3rd column, whereas the <tt>bind</tt> function only creates a new view to the existing data in <tt>array2d</tt>. 
    
    At this point we have to distinguish between the classes <tt> MultiArray </tt> and 
    <tt> MultiArrayView </tt>. MultiArray inherits from MultiArrayView and contains the
    memory management of the array. With MultiArrayView we can view the data stored in a 
    MultiArray. The code above produces a copy of the 3rd column of intArray. If we change the
    elements of <tt>lowArray</tt> nothing happens to <tt> intArray </tt>.

    \code
         // initialize 200x100 image
    MultiArray<2, int> array2d(Shape2(200,100)); // zero initialized

    // initialize new 1D array with 3rd column of a 2D array
    MultiArray<1, int> array1d = array2d.bind<0>(2);
    
    // overwrite element [0] of array1d
    array1d[0] = 1;
    
    // this has no effect on the original array2d
    // output: 0 1
    std::cout << array2d(2, 0) << " " << array1d[0] << std::endl;
    
    // initialize a view and overwrite element [0]
    MultiArrayView<1, int> array_view = array2d.bind<0>(2);
    array_view[0] = 2;
    
    // now, the original array2d has changed as well
    // output: 2 2
    std::cout << array2d(2, 0) << " " << array_view[0] << std::endl;
    \endcode

    Moving on to image processing we'll give an example how you can flip an image by using 
    bind. We read a gray scale image into a 2-dimensional array called <tt> imageArray </tt>. 
    Then we initalize a new array <tt> newImageArray </tt> of the same dimension and size
    and set the first row of <tt> newImageArray </tt> to the values of the last row of 
    <tt> imageArray </tt>, the second row to the values of the second last row and so on.
    Hence, we flip the image top to bottom.

    \dontinclude mirror_tutorial.cxx
    \skip // mirror the image horizontally 
    \until }
    
    However, you don't need to implement a method like this yourself because VIGRA already provides the 
    function \ref reflectImage(). We use this function to flip the image left-to-right:
    
    \dontinclude mirror_tutorial.cxx
    \skip // mirror the image vertically
    \until reflectImage
   
    The complete example can be found in <a href="mirror_tutorial_8cxx-example.html">mirror_tutorial.cxx</a>.
    (This program needs an infile and an outfile as command-line arguments and contains additional I/O code 
    which will be explained in section \ref ImageInputOutputTutorial.) Here you can see what happens to an input file:

    <Table cellspacing = "10">
    <TR valign="bottom">
    <TD> \image html lenna_small.gif "input file" </TD>
    <TD> \image html lenna_mirror_horizontal.gif "mirrored top to bottom" <TD>
    <TD> \image html lenna_mirror_vertical.gif "mirrored left to right" <TD>
    </TR>
    </Table>

    For completeness, there are five additional versions of the bind()-method:
    
    <DL>
    <DT><b> bindInner(i) </b> with scalar or multi-dimensional index i:</DT>
        <DD> if i is an <tt> integer </tt>, the innermost dimension (axis 0) is fixed to i, <br>
             if i is <tt>MultiArrayShape<M>::type</tt> (a shape of size M), then the M innermost 
             dimensions (axes 0...M-1) are fixed to the values in the shape vector </DD>
    <DT><b> bindOuter(i) </b> with scalar or multi-dimensional index i:</DT>
        <DD> if i is an <tt> integer </tt>, the outmost dimension (axis N-1) is fixed to i, <br>
        if i is <tt>MultiArrayShape<M>::type</tt> (a shape of size M), then the M outmost dimensions 
        (axes N-M ... N-1) are fixed to the values in the shape vector </DD>
    <DT><b> diagonal() </b>:</DT>
        <DD> Create a 1-dimensional view to the diagonal elements of the original array 
        (i.e. <tt>view[i] == array(i,i,i)</tt> for a 3D original array). </DD>
    </DL>
    
    The opposite of binding - inserting a new axis - is also possible. However, since we cannot alter the internal memory layout and thus cannot insert additional data elements, a new axis must be singleton axis, i.e. an axis with length 1. The argument of <tt>insertSingletonDimension(k)</tt> determines the position of the new axis, with <tt>0 <= k <= N</tt> when the original array has <tt>N</tt> dimensions:
    
    \code
    MultiArray<2, int> array(20,10);
    std::cout << array.insertSingletonDimension(1).shape() << "\n"; // prints "(20, 1, 10)"
    \endcode
    
    <hr>

    \subsection MultiArray_vector_elements expandElements(k) and bindElementChannel(i)
    
    When the array elements are vectors (i.e. \ref vigra::TinyVector or \ref vigra::RGBValue), we can expand these elements into an addtional array dimension:
    \code
    MultiArray<2, TinyVector<int, 3> > vector_array(20, 10);
    std::cout << vector_array.shape() << "\n";  // prints "(20, 10)"
    
    MultiArrayView<3, int> expanded(vector_array.expandElements(2));
    std::cout << expanded.shape() << "\n";  // prints "(20, 10, 3)"
    \endcode
    
    The argument <tt>k</tt> of <tt>expandElements(k)</tt> determines the desired position of the channel axis, i.e. the index that refers to the vector elements. When the original vector array has <tt>N</tt> dimensions (not counting the channel axis), it is required that <tt>0 <= k <= N</tt>. 
    
    Often, we are only interested in a single channel of a vector-valued array. This can be achieved with the function <tt>bindElementChannel(i)</tt>. For example, we can extract the green channel (i.e. channel 1) from an RGB image like this:
    \code
    MultiArray<2, RGBValue<UInt8> > rgb_array(20, 10);
    MultiArrayView<2, UInt8> green_channel(rgb_array.bindElementChannel(1));
    \endcode
    This is simply an abbreviation for <tt>rgb_array.expandElements(0).bindInner(1)</tt>.

    <hr>

    \subsection MultiArray_transpose transpose()
    Everyone is familiar with the <tt>transpose()</tt> function of a matrix (i.e. a 2-dimensional array). Once again, this operation is possible without copying the data by just manipulating the internal access functions. The following example demonstrates the difference between a transposed copy and view:

    \dontinclude transpose.cxx
    \skip  create array
    \until set transarrayView to 5
    \until print

    The output is:

    \verbatim
    base_array:
    0  1  2  3
    0  1  2  3
    0  1  2  3
    0  1  2  3
    transarray:
    0  0  0  0
    1  1  1  1
    2  2  2  2
    3  3  3  3
    transArrayView:
    0  0  0  0
    1  1  1  1
    2  2  2  2
    3  3  3  3
    base_array after setting transarray to 5
    (no change, since transarray is a copy):
    0  1  2  3
    0  1  2  3
    0  1  2  3
    0  1  2  3
    base_array after setting transarrayView to 5
    (base_array changes because transarrayView is a view):
    5  5  5  5
    5  5  5  5
    5  5  5  5
    5  5  5  5
    \endverbatim

    The function <tt>MultiArrayView::transpose()</tt> generalizes transposition to arrays of arbitrary dimensions. Here, it just reverses the order of the axes: axis 0 becomes axis N-1, axis 1 becomes axis N-2 and so on. In the following example we transpose a 5D array and print out its shape.

    \dontinclude transpose.cxx
    \skip transposing a 5D array
    \until arrayview5D.shape

    The output is:

    \verbatim
    Shape of array5D: (1, 2, 3, 4, 5)
    Shape of array5D view after default transpose(): (5, 4, 3, 2, 1)
    \endverbatim
    
    Finally, <tt>MultiArrayView::transpose()</tt> can also be called with a shape object that specifies the desired permutation of the axes: When <tt>permutation[k] = j</tt>, axis <tt>j</tt> of the original array becomes axis <tt>k</tt> of the transposed array (remember, that VIGRA counts the axes from 0):
    
    \dontinclude transpose.cxx
    \skip transpose to an explicitly specified axis permutation
    \until applied permutation
    
    The permutation in the example is 2,1,3,4,0. Thus, original dimension 0 appears in the last position of the new view, original dimension 2 appears in the first position, and so on as demonstrated by the output of the example:
    
    \verbatim
    Shape of array5D view after user-defined transpose(): (3, 2, 4, 5, 1)
        (applied permutation 2 => 0, 1 => 1, 3 => 2, 4 => 3, 0 => 4 to the axes)
    \endverbatim

    When we transpose an image about the major diagonal, we can simply use the view created by <tt>MultiArrayView::transpose()</tt>. However, transposition about the minor diagonal requires a new image, which can be filled by \ref transposeImage() like this:

    \include transpose_image_tutorial.cxx

    The result is:

    <table cellspacing = "10">
    <TR valign="bottom">
    <TD> \image html lenna_small.gif "input file" </TD>
    <TD>\image html lenna_transposed_major.gif "transpose about major diagonal" </TD>
    <TD>\image html lenna_transposed_minor.gif "transpose about minor diagonal" </TD>
    </TR>
    </table>

    In VIGRA, image transposition is also implemented in function <tt>vigra::transposeImage(...)</tt>. The difference
    is that transposeImage() copies the image data, whereas MultiArrayView::transpose() just changes the internal mapping from indices to memory locations.

    <b>Important note:</b> Transposing an array also changes the direction of the StridedScanOrderIterator. Imagine a 3x4-
    matrix. Scan-order means that we iterate from left to right, row by row. Now, let's transpose the matrix to a 4x3 view. Than, scan-order in the new view is again left to right, row by row. However, in the original matrix this now corresponds to a transposed scan: from top to bottom, column by column. The same applies to the array's index operator with integer argument.
    
    <hr>

    \subsection MultiArray_unstrided isUnstrided(k)
    
    A MultiArray always accesses its elements in consecutive memory order, i.e. <tt>&array[i] == &array.data()[i]</tt> for all <tt>i</tt> in the range <tt>[0, array.size())</tt>. However, this does in general not hold for MultiArrayViews, because changing array access is the whole point of view creation. Sometimes, it is necessary to find out if a view still has consecutive, unstrided memory access, for example when you want to pass on the view's data to an external library that only accepts plain C arrays: When the view happens to be unstrided, you can avoid to create a copy of the data. You can determine this with the function <tt>isUnstrided(k)</tt> which returns <tt>true</tt> when the array is unstrided up to dimension <tt>k</tt> (<tt>k</tt> defaults to <tt>N-1</tt>, i.e. the entire array must be unstrided):
    \code 
    MultiArray<2, int> array(20,10);
    std::cout << array.isUnstrided() << " " << array.transpose().isUnstrided() << "\n";  // prints "true false"
    \endcode
*/


/** \example mirror_tutorial.cxx
    Mirror an image file (horizontal or vertical)
    <br>
    Usage: <TT>mirror_tutorial infile outfile</TT>
*/

/** \example transpose.cxx
    Demonstrate MultiArrayView::transpose()
    <br>
    Usage: <TT>transpose</TT>
*/

/** \example transpose_image_tutorial.cxx
    Transpose an image file about major or minor axis
    <br>
    Usage: <TT>transpose_image_tutorial infile outfile</TT>
*/

/** \example subimage_tutorial.cxx
    Create a half-sized subimage
    <br>
    Usage: <TT>subimage_tutorial infile outfile</TT>
*/


/** \page ImageInputOutputTutorial Image Input and Output 

    <h2>Section Contents</h2>
    <ul style="list-style-image:url(documents/bullet.gif)">
    <li> \ref Impex2D
    <li> \ref ImpexND
    </ul>

    \section Impex2D Two-Dimensional Images
    
    In this section we'll show you how to import and export an image with VIGRA. If you
    want to import an image from disk and enquire about its properties, you must use an 
    object of <tt>vigra::ImageImportInfo</tt> class. It reads the header of the image file. 
    The constructor expects the file name, the file type will be determined automatically. 
    
    The <tt>vigra::ImageImportInfo</tt> class currently recognizes the following file formats:

    <DL>
    <DT><b>BMP:</b></dt>
              <DD> Microsoft Windows bitmap image file.</DD>
    <DT><b>EXR:</b></dt>
              <DD> OpenEXR high dynamic range image format. (only available if libopenexr is installed)</DD>
    <DT><b>GIF:</b></dt>
              <DD> CompuServe graphics interchange format; 8-bit color.</DD>
    <DT><b>HDR:</b></dt>
              <DD> Radiance RGBE high dynamic range image format.</DD>
    <DT><b>JPEG:</b></dt>
              <DD> Joint Photographic Experts Group JFIF format - compressed 24-bit color (only available if libjpeg is installed).</DD>
    <DT><b>PNG:</b></dt>
              <DD> Portable Network Graphic (only available if libpng is installed).</DD>
    <DT><b>PBM:</b></dt>
              <DD> Portable bitmap format (black and white).</DD>
    <DT><b>PGM:</b></dt>
              <DD> Portable graymap format (gray scale).</DD>
    <DT><b>PNM:</b></dt>
              <DD> Portable anymap.</DD>
    <DT><b>PPM:</b></dt>
              <DD> Portable pixmap format (color).</DD>
    <DT><b>SUN:</b></dt>
              <DD> SUN Rasterfile.</DD>
    <DT><b>TIFF:</b></dt>
              <DD> Tagged Image File Format. (only available if libtiff is installed.)</DD>
    <DT><b>VIFF:</b></dt>
              <DD> Khoros Visualization image file.</DD>
    </DL>

    In the following example, the image file name is given in the first command line argument, and the most important image metadata are printed:

    \dontinclude imageImportInfo_tutorial.cxx
    \skip read image
    \until numBands
    

    As you can see, the <tt> ImageImportInfo </tt> object contains a lot of information,
    some of it is printed in the example. Using this image

    \image html lenna_small.gif "input file"

    we get the following output:

    \verbatim
    Image information:
      file format: GIF
      width:       154
      height:      145
      pixel type:  UINT8
      color image: no  (number of channels: 1)
    \endverbatim

    To process the image, we must load the actual image data into an array such as the one described in \ref MultiArrayBasics. To do so, we create a <tt>vigra::MultiArray</tt> with the appropriate shape and then call \ref vigra::importImage(). This function needs an <tt> ImageImportInfo </tt> object specifying the image to be loaded and a MultiArrayView object of apropriate size to copy the image data in. The code looks like this:

    \dontinclude imageIO_tutorial.cxx
    \skip read image
    \until importImage(imageInfo
    
    If you already know the type of data in the file, you can also just pass the filename and a MultiArray, which will automatically be resized as appropariate:

    \dontinclude imageIO_tutorial.cxx
    \skip if you don't need
    \until importImage(in_filename
    
    Writing the image data from an array to a file is quite similar. For this purpose, you use the function \ref vigra::exportImage(), which takes an 2D MultiArrayView and an \ref vigra::ImageExportInfo object or a string (the filename). The ImageExportInfo object also needs a file name, but gives you more control over how the image is written. The desired file format is guessed from the file name's extension (but can be overridden with the method <tt>ImageExportInfo::setFileType</tt>. Recognized extensions are: '.bmp', '.exr', '.gif', '.jpeg', '.jpg', '.p7', '.png', '.pbm', '.pgm', '.pnm', '.ppm', '.ras', '.tif', '.tiff', '.xv', '.hdr' (as for reading, '.exr' requires libopenexr, '.jpg' requires libjpeg, '.png' requires libpng and '.tif' requires libtiff). In the following example, we create and save a 160x160 pixels image, where the image is a checkerboard. The image is saved as "testimage.gif" in the same folder as the executed code.

    \include imageExportInfo_tutorial.cxx

    The resulting images are the following:

    <Table cellspacing = "10">
    <TR valign="bottom">
    <TD> \image html testimage.jpg "testimage.jpg" </TD>
    <TD> \image html testimage.gif "testimage.gif" </TD>
    </TR>
    </Table>

    Finally, we give a complete example of importing, editing and exporting an image.
    After importing, we set every other horizontal line to black. This can be done
    with the <tt>bind<M>(i)</tt> method explained in \ref MultiArray_bind.
    Input and output file names are specified via command line arguments.

    \include imageIO_tutorial.cxx

    The input image and the resulting output image are:
    
    <Table cellspacing = "10">
    <TR valign="bottom">
    <TD> \image html lenna_small.gif "input image" </TD>
    <TD> \image html lenna_stripes.gif "output image" </TD>
    </TR>
    </Table>

    The handling of <b>color images</b> is exactly the same, but instead of instantiating
    a <tt>vigra::MultiArray<2, UInt8></tt> you need a
    <tt> vigra::MultiArray<2, vigra::RGBValue<UInt8> > </tt> array as described in
    \ref MultiArrayBasics. Images with alpha channel are supported by \ref importImageAlpha() and \ref exportImageAlpha().

    Note that image processing often requires more complicated calculations than in these examples.
    In this case, it is better to import and convert the data into a <tt>float</tt> array (i.e. <tt>vigra::MultiArray<2, float></tt>) instead of the simple <tt>unsigned char</tt> type in order to minimize rounding errors. When a file is imported into such an array, the conversion is automatically performed by the importImage() function. When an array is to be exported, the handling of <tt>float</tt> depends on the file format: If the file format supports float (currently: TIFF and VIFF), the data are written verbatim (unless this is explicitly overridden, see below). Otherwise, the data are mapped to <tt>unsigned char</tt> via a linear transform of the orginal range, followed by rounding (use \ref vigra::linearRangeMapping() to override this behavior by an explicit user-defined mapping).

    The <tt>ImageExportInfo</tt> class provides a number of additional methods to <b>customize data export</b>, including:
    
    <DL>
    <DT><b>setCompression():</b></dt>
              <DD>Request compressed storage if the file format supports it.</DD>
    <DT><b>setFileType():</b></dt>
              <DD>Provide the file format explicitly instead of guessing from the filename extension.</DD>
    <DT><b>setPixelType():</b></dt>
              <DD>Request pixels to be stored as the given type (results in an exception if the type is unsupported by the file format).</DD>
    <DT><b> setXResolution(),  setYResolution:</b></dt>
              <DD>Store resolution information for the two axes (ignored if unsupported by the file format).</DD>
    </DL>
    
    See \ref vigra::ImageExportInfo for a complete list and more details.

    \section ImpexND Higher Dimensional Arrays
    
    The recommended file format for arrays of arbitrary dimension is <a href="http://www.hdfgroup.org/HDF5/">HDF5</a>. It supports all possible pixel types, arbitrary dimensions, on-the-fly compression, arbitrary many arrays per file, and flexible metadata storage along with arrays. See \ref VigraHDF5Impex for more information.

    The functions \ref importVolume() and \ref exportVolume() support three additional methods to read and write 3D volume data:
    <ul>
    <li> Multipage TIFF (read/write) and Andor SIF (read only) file formats.
    <li> Image stacks (consisting of one 2D image per slice - read/write).
    <li> Raw data if accompanied with an additional text file providing the necessary metadata (shape, pixel type etc.) - read only.
    </ul>
*/

/** \example imageImportInfo_tutorial.cxx
    Read an image and print its header information
    <br>
    Usage: <TT>imageImportInfo_tutorial infile</TT>
*/

/** \example imageExportInfo_tutorial.cxx
    Create an image and export it to a file
    <br>
    Usage: <TT>imageExportInfo_tutorial</TT>
*/

/** \example imageIO_tutorial.cxx
    A complete image I/O example
    <br>
    Usage: <TT>imageIO_tutorial infile outfile</TT>
*/

/** \page MultiArrayArithmeticTutorial Mathematics with Multi-Dimensional Arrays
    
    VIGRA supports various way to perform mathematical operations (arithmetic and algebraic functions, linear alebra) on arrays. Most of these functions operate element-wise.
    
    <ul style="list-style-image:url(documents/bullet.gif)">
    <li> \ref MultiMathModule "Array Expressions"
        <BR>The \ref MultiMathModule "vigra::multi_math" module overloads the usual arithmetic operators and algebraic functions for array arguments, similar to Matlab and numpy. This leads to very efficient and readable code.
    <li> \ref LinearAlgebraModule "Linear Algebra"
        <BR>The \ref LinearAlgebraModule "vigra::linalg" module implements linear algebra for 2-dimensional arrays. The main difference to multi_math (besides a different internal implementation) is that the multiplication operator realizes matrix multiplication here. In addition, this modul implements linear system solvers, eigenvalue decomposition, and other standard matrix algorithms.
    <li> \ref MultiPointoperators "STL-style transformation algorithms"
        <BR>VIGRA also provides functions like <tt>transformMultiArray()</tt> that generalize the corresponding STL functions to multiple dimensions. The functors needed for these functions are most easily created with the \ref FunctorExpressions module, VIGRA's "lambda library". This approach offers more flexibility than the array expressions above.
    <li> \ref FeatureAccumulators
        <BR>The \ref FeatureAccumulators "vigra::acc" module provides powerful and efficient methods to compute statistics accross entire arrays or arbitrary subparts of them.
    </ul>    
*/

/** \page OwnFunctionsTutorial 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 <i>not</i> 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 <a href="smooth_explicitly_8cxx-example.html">smooth_explicitly.cxx</a>. It makes sense to encapsulate the smoothing algorithm into a function of its own. When we only need to support <tt>float</tt> 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 <tt>vigra_precondition()</tt> for \ref ErrorReporting): 
    
    \code
    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();
            }
        }
    }
    \endcode
    
    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 <tt>N</tt> as a template parameter. Then we can no longer use <tt>Shape2</tt> because this class only works for 2-dimensional arrays. Instead, we use the <tt>MultiArrayShape</tt> 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 \ref vigra::MultiCoordinateIterator from <tt>multi_iterator_coupled.hxx</tt> that iterates over all coordinates of an array, regardless of dimension. The current coordinate is returned by dereferencing the iterator:
    
    \code
    #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();
        }
    }
    \endcode
    
    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 <tt>unsigned char</tt>, and we don't want to force the user to convert them into <tt>float</tt> 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 <tt>window.sum<...>()</tt> 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 <tt>unsigned char</tt>: 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 <tt>RealPromote</tt> type in VIGRA's \ref NumericTraits "NumericTraits" class:
   
    \code
    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();
        }
    }
    \endcode
    
    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 \ref vigra::MultiArrayView, but in terms of \ref ImageIterators "image iterators" and \ref MultiIteratorGroup "hierarchical iterators". However, these iterators are more difficult to use, so the MultiArrayView approach is recommended for new code.
*/

/** \page PythonBindingsTutorial VIGRA Python Bindings
    
    See also the full <a href="../vigranumpy/index.html">vigranumpy reference</a>!
    
    When you configure VIGRA with the option <tt>-DWITH_VIGRANUMPY=1</tt> while running cmake, a Python module <tt>vigra</tt> will be compiled and installed. It exposes most of VIGRA's functionality for easy scripting and prototyping in Python. Most importantly, VIGRA's Python bindings are fully integrated with the popular 'numpy' package so that you can call vigra functions directly with numpy <tt>ndarrays</tt>. No explicit or implicit conversion of data formats is required.
    
    The syntax of the Python version is usually very similar to the C++ syntax, with one important difference: You do not have to pass pre-allocated result arrays to the functions. That is, while the call to <tt>gaussianSmoothing()</tt> in C++ is written like this
    
    \code
    MultiArray<2, float> inputImage(Shape2(width, height)), 
                         resultImage(inputImage.shape());  // pre-allocate result with correct shape
    ... // fill inputImage
    
    // smooth image with Gaussian filter with sigma=1.5
    // (pre-allocated resultImage must be passed to the function)
    gaussianSmoothing(inputImage, resultImage, 1.5);
    \endcode
    
    the corresponding Python call is
    
    \code
    >>> import numpy, vigra
    
    >>> inputImage = numpy.zeros((width, height), dtype=numpy.float32)
    ...  # fill inputImage
    
    >>> resultImage = vigra.filters.gaussianSmoothing(inputImage, 1.5);
    \endcode
    
    The result image is automatically allocated and returned by the function. Nonetheless, it is still possible to pass a result array of appropriate shape explicitly by means of the <tt>out</tt> parameter:
    
    \code
    >>> resultImage = numpy.zeros(inputImage.shape, dtype=numpy.float32)
    >>> vigra.filters.gaussianSmoothing(inputImage, 1.5, out=resultImage)
    \endcode
    
    This is, for example, useful when the same result image should be reused in several calls of the same function to avoid the repeated creation of new result arrays. Another possible use is the application of a function to only a rectangular region-of-interest: When the full result array is already allocated, you can pass a view of the approriate subarray to the <tt>out</tt> parameter in order to fill just the desired ROI.
    
    When a C++ function provides options, they are exposed on the Python side as keyword arguments:
    
    \code
    >>> labeling, max_label = vigra.analysis.watersheds(inputImage, seeds=seedImage, method='UnionFind')
    \endcode
    
    In general, the correspondence between a Python function and its C++ counterpart is straightforward, and the Python documentation frequently refers to the C++ documentation for details. However, there is a crucial difference: the default axis interpretation is different in VIGRA's <tt>MultiArray</tt> (which interpretes axes as x, y, z, so called 'Fortran' order) and in numpy's <tt>ndarray</tt> (which interpretes them as z, y, x, so called 'C'-order). To help you deal with this difficulty, vigranumpy provides a subclass <a href="../vigranumpy/index.html#axistags-and-the-vigraarray-data-structure">VigraArray</a> of <tt>ndarray</tt> and the concept of <a href="../vigranumpy/index.html#more-on-the-motivation-and-use-of-axistags">axistags</a>. Please take the time to read this material in order to avoid surprises.
    
    The full <a href="../vigranumpy/index.html">vigranumpy reference</a> is available via HTML or can be obtained directly at the Python prompt by the <tt>help()</tt> command:
    
    \code
    >>> help(vigra.filters.gaussianSmoothing)
    \endcode
    
    Another important difference between C++ and Python is that vigranumpy exposes most functions only for a restricted set of pixel types. This restriction is necessary because support for all possible type combinations would result in a combinatorial explosion and unreasonably large Python modules. In general, all functions are implemented for <tt>float</tt> pixel types (called <tt>numpy.float32</tt> on the Python side), and some provide <tt>uint8</tt> and/or <tt>uint32</tt> versions in addition. If you call a function with an unsupported pixel type, an error message listing the supported types will be printed:
    
    \code
    >>> a = vigra.ScalarImage((20,20), dtype=numpy.float64)

    >>> vigra.gaussianSmoothing(a, 1) # doesn't support numpy.float64
    ArgumentError: Python argument types in
        vigra.filters.gaussianSmoothing(numpy.ndarray, int)
    did not match C++ signature:
        gaussianSmoothing(class vigra::NumpyArray<3,struct vigra::Multiband<float>,struct vigra::StridedArrayTag> array, class boost::python::api::object sigma, class vigra::NumpyArray<3,struct vigra::Multiband<float>,struct vigra::StridedArrayTag> out=None, class boost::python::api::object sigma_d=0.0, class boost::python::api::object step_size=1.0, double window_size=0.0, class boost::python::api::object roi=None)    

        gaussianSmoothing(class vigra::NumpyArray<4,struct vigra::Multiband<float>,struct vigra::StridedArrayTag> array, class boost::python::api::object sigma, class vigra::NumpyArray<4,struct vigra::Multiband<float>,struct vigra::StridedArrayTag> out=None, class boost::python::api::object sigma_d=0.0, class boost::python::api::object step_size=1.0, double window_size=0.0, class boost::python::api::object roi=None)
    \endcode
    
    The error message is automatically generated by boost::python and therefore rather technical. It says that <tt>%gaussianSmoothing()</tt> supports 3- and 4-dimensional arrays where the rightmost dimension is interpreted as a channel axis, and the pixel type must be <tt>float</tt> (these properties are indicated by the type specifications <tt>NumpyArray<3,struct vigra::Multiband<float></tt> and <tt>NumpyArray<4,struct vigra::Multiband<float></tt> respectively). Thus, the input array must be a <tt>float32</tt> image or volume with either no explicit channel axis (in which case a singleton channel axis will be inserted automatically) or with arbitrary many channels (e.g. RGB). 
    
    <a href="../vigranumpy/index.html#more-on-the-motivation-and-use-of-axistags">Axistags</a> allow vigranumpy to distinguish if a given 3-dimensional array is to be interpreted as a 2D image with multiple channels, or as a 3D volume with only a single channel. If no axistags are attached to the array, it is unspecified which version of an algorithm will be called. Axistags are automatically specified when arrays are created with one of the factory functions in the <tt>vigra</tt> module, for example:
    
    \code
    >>> a = vigra.ScalarImage((30, 20))
    >>> print("%s \n %r" % (a.shape, a.axistags))
    (30L, 20L)
     x y

    >>> a = vigra.RGBImage((30, 20))
    >>> print("%s \n %r" % (a.shape, a.axistags))
    (30L, 20L, 3L)
     x y c

    >>> a = vigra.ScalarVolume((30, 20, 10))
    >>> print("%s \n %r" % (a.shape, a.axistags))
    (30L, 20L, 10L)
     x y z

    >>> a = vigra.RGBVolume((30, 20, 10))
    >>> print("%s \n %r" % (a.shape, a.axistags))
    (30L, 20L, 10L, 3L)
     x y z c
    \endcode
    
    Axistags are encoded 'x', 'y', 'z' for the three spatial axes, 'c' for a channel axis, and 't' for a time axis. If the channel axis is missing, vigranumpy will assume that the array has only a single channel. That is, arrays with shape (30, 20, 1) and axistags 'x y c' are equivalent to arrays with shape (30, 20) and axistags 'x y'. Functions that change the order of the axes (such as <tt>%array.transpose()</tt>) or reduce the number of axes (e.g. <tt>array[:, 1, :]</tt>) also modify the axistags accordingly, so that you can always ask for the axis meaning by simply calling <tt>array.axistags</tt>.
*/