Switch to using C-order for image data (#9)

Within the realm of computing, there are two ways of storing multidimensional data:
* **C-order:** Gets its name from C where you typically create a multidimensional array like so `int data[z][y][x]`. This is backwards from traditional mathematics where you represent coordinates as (x, y, z). This format of storing data can also be thought of as slowest-varying to fastest-varying. The z-axis typically varies slowest, followed by y-axis, finished with the x axis.
* **Fortran-order:** Direct opposite of C-order and is how MATLAB and Fortran handle multidimensional array indexing. This method follows the traditional mathematical representation of (x, y, z). It makes sense that MATLAB and Fortran use this method because they are languages geared towards data scientists. This format of indexing data can be thought of as fastest-varying to slowest-varying. 

Previously, this package wasn't using C-order or Fortran-order but did a weird collage of both. This could be confusing to new users and I wanted to stick to a particular syntax.

# Why switch to C-order?

Python, Numpy, and many other libraries use this format for storing and indexing data. Numpy provides support for changing the order but it can be a pain to include for each call. Thus, I decided to stick with this syntax that Python follows.

# What changes were made?

Here is the general thought process I used. When returning coordinates, Fortran-order was used because that is intuitive for the average user who is familiar with coordinates from math courses. For example, it is much more intuitive to use (x, y) for a point rather than (y, x). The polar domain equivalent is (r, theta).

Images and image shapes or sizes are in C-order to stick with the overwhelming agreement in the Python community. 

Color channels are considered to be faster varying than the x axis.

For example, let's assume you have 50 RGB images with a width of 1024px and height of 2048px. The Numpy array of the images should be (50, 2048, 50, 3). This is what `polarTransform` expects when given an image. Let's say that you convert it to the polar domain, the resulting image will have a shape of (50, angular size, radial size, 3).

Color images with channels > 1 must be specified by using the parameter `hasColor`. This is handled by moving the color channel to the front of the dimensions, performing transformation, and then moving back after transformation is complete.

docs/source/user-guide.rst

@@ -10,6 +10,8 @@ As the names suggest, the two functions convert an image from the cartesian or p
 
 Since there are quite a few parameters that can be specified for the conversion functions, the class :class:`ImageTransform` is created and returned from the :class:`convertToPolarImage` or :class:`convertToCartesianImage` functions (along with the converted image) that contains the arguments specified. The benefit of this class is that if one wants to convert the image back to another domain or convert points on either image to/from the other domain, they can simply call the functions within the :class:`ImageTransform` class without specifying all of the arguments again.
 
+The examples below use images from the test suite. The code snippets should run without modification except for changing the paths to point to the correct image.
+
 Example 1
 --------------
 Let us take a B-mode echocardiogram and convert it to the polar domain. This is essentially reversing the scan conversion done internally by the ultrasound machine.
@@ -28,13 +30,19 @@ Here is the B-mode image:
     cartesianImage = imageio.imread('IMAGE_PATH_HERE')
 
     polarImage, ptSettings = polarTransform.convertToPolarImage(cartesianImage, center=[401, 365])
-    plt.imshow(polarImage, origin='lower')
+    plt.imshow(polarImage.T, origin='lower')
 
 Resulting polar domain image:
 
 .. image:: _static/shortAxisApexPolarImage.png
     :alt: Polar image of echocardiogram of short-axis apex view
 
+The example input image has a width of 800px and a height of 604px. Since many imaging libraries use C-order rather than Fortran order, the Numpy array containing the image data loaded from imageio has a shape of (604, 800). This is what polarTransform expects for an image where the first dimension is the slowest varying (y) and the last dimension is the fastest varying (x). Additional dimensions can be present before the y & x dimensions in which case polarTransform will transform each 2D slice individually.
+
+The center argument should be a list, tuple or Numpy array of length 2 with format (x, y). A common theme throughout this library is that points will be specified in Fortran order, i.e. (x, y) or (r, theta) whilst data and image sizes will be specified in C-order, i.e. (y, x) or (theta, r).
+
+The polar image returned in this example is in C-order. So this means that the data is (theta, r). When displaying an image using matplotlib, the first dimension is y and second is x. The image is transposed before displaying to flip it 90 degrees.
+
 Example 2
 --------------
 Input image:
@@ -52,12 +60,12 @@ Input image:
 
     polarImage, ptSettings = polarTransform.convertToPolarImage(verticalLinesImage, initialRadius=30,
                                                                 finalRadius=100, initialAngle=2 / 4 * np.pi,
-                                                                finalAngle=5 / 4 * np.pi)
+                                                                finalAngle=5 / 4 * np.pi, hasColor=True)
 
     cartesianImage = ptSettings.convertToCartesianImage(polarImage)
 
     plt.figure()
-    plt.imshow(polarImage, origin='lower')
+    plt.imshow(polarImage.T, origin='lower')
 
     plt.figure()
     plt.imshow(cartesianImage, origin='lower')
@@ -70,4 +78,8 @@ Resulting polar domain image:
 Converting back to the cartesian image results in:
 
 .. image:: _static/verticalLinesCartesianImage_scaled.png
-    :alt: Cartesian image
+    :alt: Cartesian image
+
+Once again, when displaying polar images using matplotlib, the image is first transposed to rotate the image 90 degrees. This makes it easier to view the image because the theta dimension is longer than the radial dimension.
+
+The hasColor argument was set to True in this example because the image contains color images. The example RGB image has a width and height of 256px. The shape of the image loaded via imageio package is (256, 256, 3). By specified hasColor=True, the last dimension will be shifted to the front and then the polar transformation will occur on each channel separately. Before returning, the function will shift the channel dimension back to the end. If a RGBA image is loaded, it is advised to only transform the first 3 channels and then set the alpha channel to fully on.

polarTransform/convertToCartesianImage.py

@@ -5,7 +5,7 @@
 
 def convertToCartesianImage(image, center=None, initialRadius=None,
                             finalRadius=None, initialAngle=None,
-                            finalAngle=None, imageSize=None, order=3, border='constant',
+                            finalAngle=None, imageSize=None, hasColor=False, order=3, border='constant',
                             borderVal=0.0, useMultiThreading=False, settings=None):
     """Convert polar image to cartesian image.
 
@@ -15,16 +15,23 @@ def convertToCartesianImage(image, center=None, initialRadius=None,
 
     Parameters
     ----------
-    image : (N, M) or (N, M, P) :class:`numpy.ndarray`
+    image : N-dimensional :class:`numpy.ndarray`
         Polar image to convert to cartesian domain
 
-        .. note::
-            For a 3D array, polar transformation is applied separately across each 2D slice
+        Image should be structured in C-order, i.e. the axes should be ordered (..., z, theta, r, [ch]). The channel
+        axes should only be present if :obj:`hasColor` is :obj:`True`. This format is arbitrary but is selected to stay
+        consistent with the traditional C-order representation in the Cartesian domain.
+
+        In the mathematical domain, Cartesian coordinates are traditionally represented as (x, y, z) and as
+        (r, theta, z) in the polar domain. When storing Cartesian data in C-order, the axes are usually flipped and the
+        data is saved as (z, y, x). Thus, the polar domain coordinates are also flipped to stay consistent, hence the
+        format (z, theta, r).
 
         .. note::
-            If an alpha band (4th channel of image is present), then it will be converted. Typically, this is
-            unwanted, so the recommended solution is to transform the first 3 channels and set the 4th channel to
-            fully on.
+            For multi-dimensional images above 2D, the cartesian transformation is applied individually across each
+            2D slice. The last two dimensions should be the r & theta dimensions, unless :obj:`hasColor` is True in
+            which case the 2nd and 3rd to last dimensions should be. The multidimensional shape will be preserved
+            for the resulting cartesian image (besides the polar dimensions).
     center : :class:`str` or (2,) :class:`list`, :class:`tuple` or :class:`numpy.ndarray` of :class:`int`, optional
         Specifies the center in the cartesian image to use as the origin in polar domain. The center in the
         cartesian domain will be (0, 0) in the polar domain.
@@ -108,6 +115,19 @@ def convertToCartesianImage(image, center=None, initialRadius=None,
 
         If imageSize is not set, then it defaults to the size required to fit the entire polar image on a cartesian
         image.
+    hasColor : :class:`bool`, optional
+        Whether or not the polar image contains color channels
+
+        This means that the image is structured as (..., y, x, ch) or (..., theta, r, ch) for Cartesian or polar
+        images, respectively. If color channels are present, the last dimension (channel axes) will be shifted to
+        the front, converted and then shifted back to its original location.
+
+        Default is :obj:`False`
+
+        .. note::
+            If an alpha band (4th channel of image is present), then it will be converted. Typically, this is
+            unwanted, so the recommended solution is to transform the first 3 channels and set the 4th channel to
+            fully on.
     order : :class:`int` (0-5), optional
         The order of the spline interpolation, default is 3. The order has to be in the range 0-5.
 
@@ -162,15 +182,25 @@ def convertToCartesianImage(image, center=None, initialRadius=None,
 
     Returns
     -------
-    cartesianImage : (N, M) or (N, M, P) :class:`numpy.ndarray`
-        Cartesian image (3D cartesian image if 3D input image is given)
+    cartesianImage : N-dimensional :class:`numpy.ndarray`
+            Cartesian image
+
+            Resulting image is structured in C-order, i.e. the axes are ordered as (..., z, y, x, [ch]). This format is
+            the traditional method of storing image data in Python.
+
+            Resulting image shape will be the same as the input image except for the polar dimensions are
+            replaced with the Cartesian dimensions.
     settings : :class:`ImageTransform`
         Contains metadata for conversion between polar and cartesian image.
 
         Settings contains many of the arguments in :func:`convertToPolarImage` and :func:`convertToCartesianImage` and
         provides an easy way of passing these parameters along without having to specify them all again.
     """
 
+    # If there is a color channel present, move it to the front of axes
+    if settings.hasColor if settings is not None else hasColor:
+        image = np.moveaxis(image, -1, 0)
+
     # Create settings if none are given
     if settings is None:
         # Center is set to middle-middle, which means all four quadrants will be shown
@@ -187,7 +217,7 @@ def convertToCartesianImage(image, center=None, initialRadius=None,
         # In other words, what radius does the last row of polar image correspond to?
         # If not set, default is the largest radius from image
         if finalRadius is None:
-            finalRadius = image.shape[0]
+            finalRadius = image.shape[-1]
 
         # Initial angle of the source image
         # In other words, what angle does column 0 correspond to?
@@ -202,15 +232,15 @@ def convertToCartesianImage(image, center=None, initialRadius=None,
             finalAngle = 2 * np.pi
 
         # This is used to scale the result of the radius to get the appropriate Cartesian value
-        scaleRadius = image.shape[0] / (finalRadius - initialRadius)
+        scaleRadius = image.shape[-1] / (finalRadius - initialRadius)
 
         # This is used to scale the result of the angle to get the appropriate Cartesian value
-        scaleAngle = image.shape[1] / (finalAngle - initialAngle)
+        scaleAngle = image.shape[-2] / (finalAngle - initialAngle)
 
         if imageSize is None:
             # Obtain the image size by looping from initial to final source angle (every possible theta in the image
             # basically)
-            thetas = np.mod(np.linspace(0, (finalAngle - initialAngle), image.shape[1]) + initialAngle, 2 * np.pi)
+            thetas = np.mod(np.linspace(0, (finalAngle - initialAngle), image.shape[-2]) + initialAngle, 2 * np.pi)
             maxRadius = finalRadius * np.ones_like(thetas)
 
             # Then get the maximum radius of the image and compute the x/y coordinates for each option
@@ -289,13 +319,13 @@ def convertToCartesianImage(image, center=None, initialRadius=None,
         imageSize = tuple(imageSize)
 
         settings = ImageTransform(center, initialRadius, finalRadius, initialAngle, finalAngle, imageSize,
-                                  image.shape[0:2])
+                                  image.shape[-2:], hasColor)
     else:
         # This is used to scale the result of the radius to get the appropriate Cartesian value
-        scaleRadius = settings.polarImageSize[0] / (settings.finalRadius - settings.initialRadius)
+        scaleRadius = settings.polarImageSize[1] / (settings.finalRadius - settings.initialRadius)
 
         # This is used to scale the result of the angle to get the appropriate Cartesian value
-        scaleAngle = settings.polarImageSize[1] / (settings.finalAngle - settings.initialAngle)
+        scaleAngle = settings.polarImageSize[0] / (settings.finalAngle - settings.initialAngle)
 
     # Get list of cartesian x and y coordinate and create a 2D create of the coordinates using meshgrid
     xs = np.arange(0, settings.cartesianImageSize[1])
@@ -320,23 +350,25 @@ def convertToCartesianImage(image, center=None, initialRadius=None,
     theta = theta * scaleAngle
 
     # Flatten the desired x/y cartesian points into one 2xN array
-    desiredCoords = np.vstack((r.flatten(), theta.flatten()))
+    desiredCoords = np.vstack((theta.flatten(), r.flatten()))
 
     # Get the new shape which is the cartesian image shape plus any other dimensions
-    newShape = settings.cartesianImageSize + image.shape[2:]
+    # Get the new shape of the cartesian image which is the same shape of the polar image except the last two dimensions
+    # (r & theta) are replaced with the cartesian image size
+    newShape = image.shape[:-2] + settings.cartesianImageSize
 
     # Reshape the image to be 3D, flattens the array if > 3D otherwise it makes it 3D with the 3rd dimension a size of 1
-    image = image.reshape(image.shape[0:2] + (-1,))
+    image = image.reshape((-1,) + settings.polarImageSize)
 
     if border == 'constant':
         # Pad image by 3 pixels and then offset all of the desired coordinates by 3
-        image = np.pad(image, ((3, 3), (3, 3), (0, 0)), 'edge')
+        image = np.pad(image, ((0, 0), (3, 3), (3, 3)), 'edge')
         desiredCoords += 3
 
     if useMultiThreading:
         with concurrent.futures.ThreadPoolExecutor() as executor:
-            futures = [executor.submit(scipy.ndimage.map_coordinates, image[:, :, k], desiredCoords, mode=border,
-                                       cval=borderVal, order=order) for k in range(image.shape[2])]
+            futures = [executor.submit(scipy.ndimage.map_coordinates, slice, desiredCoords, mode=border, cval=borderVal,
+                                       order=order) for slice in image]
 
             concurrent.futures.wait(futures, return_when=concurrent.futures.ALL_COMPLETED)
 
@@ -345,13 +377,17 @@ def convertToCartesianImage(image, center=None, initialRadius=None,
         cartesianImages = []
 
         # Loop through the third dimension and map each 2D slice
-        for k in range(image.shape[2]):
-            imageSlice = scipy.ndimage.map_coordinates(image[:, :, k], desiredCoords, mode=border, cval=borderVal,
+        for slice in image:
+            imageSlice = scipy.ndimage.map_coordinates(slice, desiredCoords, mode=border, cval=borderVal,
                                                        order=order).reshape(x.shape)
             cartesianImages.append(imageSlice)
 
     # Stack all of the slices together and reshape it to what it should be
-    cartesianImage = np.dstack(cartesianImages).reshape(newShape)
+    cartesianImage = np.stack(cartesianImages, axis=0).reshape(newShape)
+
+    # If there is a color channel present, move it abck to the end of axes
+    if settings.hasColor:
+        cartesianImage = np.moveaxis(cartesianImage, 0, -1)
 
     return cartesianImage, settings