Hi, On Mon, Aug 3, 2015 at 8:09 AM, Nathaniel Smith <njs@pobox.com> wrote:
On Aug 2, 2015 11:06 PM, "Kang Wang" <kwang24@wisc.edu> wrote:
This is very good discussion. Thank you all for replying.
I can see the fundamental difference is that I always think/talk/read/write a 3D image as I(x, y, z), not (plane, row, column) . I am coming from MRI (magnetic resonance imaging) research, and I can assure you that the entire MRI community is using (x, y, z), including books, journal papers, conference abstracts, presentations, everything. We even talk about what we called "logical x/y/z" and "physical x/y/z", and the rotation matrix that converts the two coordinate systems. The radiologists are also used to (x, y, z). For example, we always say "my image is 256 by 256 by 20 slices", and we never say "20 by 256 by 256".
So, basically, at least in MRI, we always talk about an image as I(x, y, z), and we always assume that "x" is the fastest changing index. That's why I prefer column-major (because it is more natural).
Of course, I can totally get my work done by using row-major, I just have to always remind myself "write last dimension index first" when coding. I actually have done this before, and I found it would be so much easier if just using column-major.
Why not just use I[x, y, z] like you're used to, and let the computer worry about the physical layout in memory? Sometimes this will be Fortran order and sometimes C order and sometimes something else, but you don't have to know or care; 99% of the time it doesn't matter. The worst case is that when you use a python wrapper to call into a library that can only handle Fortran order, then the wrapper will quietly have to convert the memory order around and it will be slightly slower than if you had used Fortran order in the first place. But in practice you'll barely ever notice this, and when you do, *then* you can tell numpy explicitly what memory layout you want in the situation where it matters.
Yes - if you are using numpy, you really have to look numpy in the eye and say: "I will let you worry about the array element order in memory, and in return, you promise to make indexing work as I would expect" Just for example, let's say you loaded an MRI image into memory: In [1]: import nibabel In [2]: img = nibabel.load('my_mri.nii') In [3]: data = img.get_data() Because NIfTI images are Fortran memory layout, this happens to be the memory layout you get for your array: In [4]: data.flags Out[4]: C_CONTIGUOUS : False F_CONTIGUOUS : True OWNDATA : False WRITEABLE : True ALIGNED : True UPDATEIFCOPY : False But now - in Python - all I care about is what data I have on the first, second, third axes. For example, I could do this: In [5]: data_copy = data.copy() This has exactly the same values as the original array, and at the same index positions: In [7]: import numpy as np In [8]: np.all(data == data) Out[8]: memmap(True, dtype=bool) but I now have a C memory layout array. In [9]: data_copy.flags Out[9]: C_CONTIGUOUS : True F_CONTIGUOUS : False OWNDATA : True WRITEABLE : True ALIGNED : True UPDATEIFCOPY : False Worse than that, if I slice my original data array, then I get an array that is neither C- or Fortran- compatible in memory: In [10]: data_view = data[:, :, ::2] In [11]: data_view.flags Out[11]: C_CONTIGUOUS : False F_CONTIGUOUS : False OWNDATA : False WRITEABLE : True ALIGNED : True UPDATEIFCOPY : False So - if you want every array to be Fortran-contiguous in memory, I would not start with numpy at all, I would write your own array library. The alternative - or "the numpy way" - is to give up on enforcing a particular layout in memory, until you need to pass an array to some C or C++ or Fortran code that needs some particular layout, in which case you get your extension code to copy the array into the required layout on entry. Of course this is what numpy itself has to do when interfacing with external libraries like BLAS or LAPACK. Cheers, Matthew