Merge pull request #1280 from jni/array_order

Describe array order terminology in 2D, 2D+c, 3D, 3D+c arrays

doc/source/user_guide/numpy_images.txt

@@ -1,4 +1,4 @@
-A crash course on Numpy for images
+A crash course on NumPy for images
 ----------------------------------
 
 Images manipulated by ``scikit-image`` are simply NumPy arrays. Hence, a
@@ -6,6 +6,8 @@ large fraction of operations on images will just consist in using NumPy::
 
     >>> from skimage import data
     >>> camera = data.camera()
+    >>> type(camera)
+    <type 'numpy.ndarray'>
 
 Retrieving the geometry of the image and the number of pixels: ::
 
@@ -21,22 +23,29 @@ Retrieving statistical information about gray values: ::
     >>> camera.mean()
     118.31400299072266
 
-Numpy arrays representing images can be of different integer of float
-numerical types. See :ref:`data_types` for more information about data
-types.
+NumPy arrays representing images can be of different integer of float
+numerical types. See :ref:`data_types` for more information about these
+types and how scikit-image treats them.
 
-Numpy indexing can be used both for looking at pixel values, and to
+
+NumPy indexing
+--------------
+
+NumPy indexing can be used both for looking at pixel values, and to
 modify pixel values: ::
 
-    >>> # Value of pixel on 10th line and 20th column
+    >>> # Get the value of the pixel on the 10th row and 20th column
     >>> camera[10, 20]
     153
-    >>> # Turn to black pixel on 3rd line and 10th column
+    >>> # Set to black the pixel on the 3rd row and 10th column
     >>> camera[3, 10] = 0
 
-Be careful that the first dimension (``camera.shape[0]``) corresponds to
-lines, while the second dimension (``camera.shape[1]``) stands for
-columns.
+Be careful: in NumPy indexing, the first dimension (``camera.shape[0]``)
+corresponds to rows, while the second (``camera.shape[1]``) corresponds
+to columns, with the origin (``camera[0, 0]``) on the top-left corner.
+This matches matrix/linear algebra notation, but is in contrast to
+Cartesian (x, y) coordinates. See `Coordinate conventions`_ below for
+more details.
 
 Beyond individual pixels, it is possible to access / modify values of
 whole sets of pixels, using the different indexing possibilities of
@@ -55,28 +64,193 @@ Masking (indexing with masks of booleans)::
 
 Fancy indexing (indexing with sets of indices) ::
 
-    >>> inds_x = np.arange(len(camera))
-    >>> inds_y = 4 * inds_x % len(camera)
-    >>> camera[inds_x, inds_y] = 0
+    >>> inds_r = np.arange(len(camera))
+    >>> inds_c = 4 * inds_r % len(camera)
+    >>> camera[inds_r, inds_c] = 0
 
 Using masks, especially, is very useful to select a set of pixels on
 which to perform further manipulations. The mask can be any boolean array
-of same shape as the image (or at least a shape broadcastable to the
-image shape). This can be useful to define a region of interest, as a
+of same shape as the image (or a shape broadcastable to the image shape).
+This can be useful to define a region of interest, such as a
 disk: ::
 
-    >>> l_x, l_y = camera.shape[0], camera.shape[1]
-    >>> X, Y = np.ogrid[:l_x, :l_y]
-    >>> outer_disk_mask = (X - l_x / 2)**2 + (Y - l_y / 2)**2 < (l_x / 2)**2
+    >>> nrows, ncols = camera.shape
+    >>> row, col = np.ogrid[:nrows, :ncols]
+    >>> cnt_row, cnt_col = nrows / 2, ncols / 2
+    >>> outer_disk_mask = ((row - cnt_row)**2 + (col - cnt_col)**2 <
+    ...                    (nrows / 2)**2)
     >>> camera[outer_disk_mask] = 0
 
 .. image:: ../../_images/plot_camera_numpy_1.png
     :width: 45%
     :target: ../auto_examples/plot_camera_numpy.html
 
-Boolean arithmetics can be used to define more complex masks: ::
+Boolean arithmetic can be used to define more complex masks: ::
 
-    >>> lower_half = X > l_x / 2
+    >>> lower_half = row > cnt_row
     >>> lower_half_disk = np.logical_and(lower_half, outer_disk_mask)
     >>> camera = data.camera()
     >>> camera[lower_half_disk] = 0 
+
+
+Color images
+------------
+
+All of the above is true of color images, too: a color image is a
+NumPy array, with an additional trailing dimension for the channels:
+
+    >>> cat = data.chelsea()
+    >>> type(cat)
+    <type 'numpy.ndarray'>
+    >>> cat.shape
+    (300, 451, 3)
+
+This shows that ``cat`` is a 300-by-451 pixel image with three
+channels (red, green, and blue).
+As before, we can get and set pixel values:
+
+    >>> cat[10, 20]
+    array([151, 129, 115], dtype=uint8)
+    >>> # set the pixel at row 50, column 60 to black
+    >>> cat[50, 60] = 0
+    >>> # set the pixel at row 50, column 61 to green
+    >>> cat[50, 61] = [0, 255, 0] # [red, green, blue]
+
+We can also use 2D boolean masks for a 2D color image, as we did with
+the grayscale image above:
+
+.. plot::
+
+    Using a 2D mask on a 2D color image
+
+    >>> reddish = cat[:, :, 0] > 160
+    >>> cat[reddish] = [0, 255, 0]
+    >>> plt.imshow(cat)
+
+
+Coordinate conventions
+----------------------
+
+Because we represent images with numpy arrays, our coordinates must
+match accordingly. Two-dimensional (2D) grayscale images (such as
+`camera` above) are indexed by row and columns (abbreviated to either
+``row, col`` or ``r, c``), with the lowest element (0, 0) at the top-
+-left corner. In various parts of the library, you will
+also see ``rr`` and ``cc`` refer to lists of row and column
+coordinates. We distinguish this from (x, y), which commonly denote
+standard Cartesian coordinates, where x is the horizontal coordinate,
+y the vertical, and the origin is on the bottom right. (Matplotlib, for
+example, uses this convention.)
+
+In the case of color (or multichannel) images, the last dimension
+contains the color information and is denoted ``channel`` or ``ch``.
+
+Finally, for 3D images, such as videos, magnetic resonance imaging
+(MRI) scans, or confocal microscopy, we refer to the leading dimension
+as ``plane``, abbreviated as ``pln`` or ``p``.
+
+These conventions are summarized below:
+
+.. table:: Dimension name and order conventions in scikit-image
+
+========================   ========================================
+Image type                 coordinates
+========================   ========================================
+2D grayscale               (row, col)
+2D multichannel (eg. RGB)  (row, col, ch)
+3D grayscale               (pln, row, col)
+3D multichannel            (pln, row, col, ch)
+========================   ========================================
+
+
+Many functions in scikit-image operate on 3D images directly:
+
+    >>> im3d = np.random.rand(100, 1000, 1000)
+    >>> from skimage import morphology
+    >>> from scipy import ndimage as nd
+    >>> seeds = nd.label(im3d < 0.1)[0]
+    >>> ws = morphology.watershed(im3d, seeds)
+
+In many cases,
+the third imaging dimension has lower resolution than the other two.
+Some scikit-image functions provide a ``spacing`` keyword argument
+to process these images:
+
+    >>> from skimage import segmentation
+    >>> slics = segmentation.slic(im3d, spacing=[5, 1, 1], multichannel=False)
+
+
+Other times, processing must be done plane-wise. When planes are the
+leading dimension, we can use the following syntax:
+
+    >>> from skimage import filters
+    >>> edges = np.zeros_like(im3d)
+    >>> for pln, image in enumerate(im3d):
+    ...     # iterate over the leading dimension (planes)
+    ...     edges[pln] = filters.sobel(image)
+
+
+Notes on array order
+--------------------
+
+Although the labeling of the axes seems arbitrary, it can have a
+significant effect on speed of operations. This is because modern
+processors never retrieve just one item from memory, but rather a
+whole chunk of adjacent items. (This is called prefetching.)
+Therefore, processing elements that are
+next to each other in memory is faster than processing them
+in a different order, even if the number of operations is the same:
+
+    >>> def in_order_multiply(arr, scalar):
+    ...     for plane in list(range(arr.shape[0])):
+    ...         arr[plane, :, :] *= scalar
+    ... 
+    >>> def out_of_order_multiply(arr, scalar):
+    ...     for plane in list(range(arr.shape[2])):
+    ...         arr[:, :, plane] *= scalar
+    ... 
+    >>> import time
+    >>> im3d = np.random.rand(100, 1024, 1024)
+    >>> t0 = time.time(); x = in_order_multiply(im3d, 5); t1 = time.time()
+    >>> print("%.2f seconds" % (t1 - t0))  # doctest: +SKIP
+    0.14 seconds
+    >>> im3d_t = np.transpose(im3d).copy() # place "planes" dimension at end
+    >>> im3d_t.shape
+    (1024, 1024, 100)
+    >>> s0 = time.time(); x = out_of_order_multiply(im3d, 5); s1 = time.time()
+    >>> print("%.2f seconds" % (s1 - s0)) doctest: +SKIP
+    1.18 seconds
+    >>> print("Speedup: %.1fx" % ((s1 - s0) / (t1 - t0)))  doctest: +SKIP
+    Speedup: 8.6x
+
+
+When the dimension you are iterating over is even larger, the
+speedup is even more dramatic. It is worth thinking about this
+*data locality* when writing algorithms. In particular, know that
+scikit-image uses C-contiguous arrays unless otherwise specified, so
+one should iterate along the last/rightmost dimension in the
+innermost loop of the computation.
+
+A note on time
+--------------
+
+Although scikit-image does not currently (0.11) provide functions to
+work specifically with time-varying 3D data, our compatibility with
+numpy arrays allows us to work quite naturally with a 5D array of the
+shape (t, pln, row, col, ch):
+
+    >>> for timepoint in image5d:  # doctest: +SKIP
+    ...     # each timepoint is a 3D multichannel image
+    ...     do_something_with(timepoint)
+
+
+We can then supplement the above table as follows:
+
+.. table:: Addendum to dimension names and orders in scikit-image
+
+========================   ========================================
+Image type                 coordinates
+========================   ========================================
+2D color video             (t, row, col, ch)
+3D multichannel video      (t, pln, row, col, ch)
+========================   ========================================